# Constraining scalar leptoquarks using COHERENT data

Roberta Calabrese<EMAIL_ADDRESS>Dipartimento di Fisica “Ettore
Pancini”, Università degli studi di Napoli “Federico II”, Complesso Univ.
Monte S. Angelo, I-80126 Napoli, Italy  INFN - Sezione di Napoli, Complesso
Univ. Monte S. Angelo, I-80126 Napoli, Italy  Jacob Gunn
<EMAIL_ADDRESS>Dipartimento di Fisica “Ettore Pancini”, Università
degli studi di Napoli “Federico II”, Complesso Univ. Monte S. Angelo, I-80126
Napoli, Italy  INFN - Sezione di Napoli, Complesso Univ. Monte S. Angelo,
I-80126 Napoli, Italy  Gennaro Miele<EMAIL_ADDRESS>Dipartimento di
Fisica “Ettore Pancini”, Università degli studi di Napoli “Federico II”,
Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy  INFN - Sezione di
Napoli, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy  Scuola
Superiore Meridionale, Università degli studi di Napoli “Federico II”, Largo
San Marcellino 10, 80138 Napoli, Italy  Stefano Morisi<EMAIL_ADDRESS>Dipartimento di Fisica “Ettore Pancini”, Università degli studi di Napoli
“Federico II”, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy  INFN -
Sezione di Napoli, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy
Samiran Roy<EMAIL_ADDRESS>INFN - Sezione di Napoli, Complesso Univ. Monte
S. Angelo, I-80126 Napoli, Italy  Pietro Santorelli
<EMAIL_ADDRESS>Dipartimento di Fisica “Ettore Pancini”,
Università degli studi di Napoli “Federico II”, Complesso Univ. Monte S.
Angelo, I-80126 Napoli, Italy  INFN - Sezione di Napoli, Complesso Univ. Monte
S. Angelo, I-80126 Napoli, Italy

###### Abstract

Neutrino-nucleus coherent scattering measurements by the COHERENT
collaboration provide us with a unique capability to test various beyond the
standard model scenarios. In this work, we constrain scalar leptoquarks (LQs)
using the COHERENT data. LQs arise in many extensions of the Standard Model
(SM). Generally, the mass of the LQs is assumed to be very high to avoid the
bounds from proton decay. However, there are low-scale LQ models which
prohibit proton decay by construction. We consider two electroweak doublet
scalar LQ models with hypercharge Y=1/6, and Y=7/6 and provide the bounds in
the plane of the Yukawa coupling and the mass of LQ. We also compare the
bounds on LQs coming from various other experiments and find that the COHERENT
one covers a wide range of LQ masses from MeV to TeV and in certain regions
the constraints are competitive with the others.

##  I Introduction

The detection of coherent elastic neutrino-nucleus scattering (CE$\nu$NS) was
finally achieved by the COHERENT collaboration Akimov  _et al._ (2017) more
than forty years after it was predicted by Freedman Freedman (1974). The
cross-section of CE$\nu$NS events depends on the momentum transfer to the
nucleus ($q^{2}$), and scales approximately as the square of the total number
of protons (Z) and neutrons (N) in the target nucleus for low $q^{2}$. In the
Standard Model (SM), CE$\nu$NS is mediated by the Z-boson, which couples
preferentially to the neutron. Hence, the SM coherent cross section depends
mainly on the square of N. The signals observed by the COHERENT collaboration
show no significant deviation from the SM prediction, for CsI Akimov  _et al._
(2017) and Ar Akimov  _et al._ (2021) nuclei.

The COHERENT data presents an opportunity to verify the SM parameters at low
energy scales ($\sim$ 10 MeV); in this energy range, the weak mixing angles
can be evaluated Miranda  _et al._ (2020a); Cadeddu  _et al._ (2020a, 2021a)
as well as the root mean square radii of the neutron in Cesium, Iodide, and
Argon Cadeddu  _et al._ (2018, 2020a); Canas  _et al._ (2020); Miranda  _et
al._ (2020a). Furthermore, the COHERENT data is also sensitive to new physics
which couples to quarks and leptons, such as extra mediators Cadeddu  _et al._
(2021b); Flores  _et al._ (2020); Banerjee  _et al._ (2021); de la Vega  _et
al._ (2021); Atzori Corona  _et al._ (2022); Binh  _et al._ (2021); Bertuzzo
_et al._ (2022); Chakraborty  _et al._ (2022); De Romeri  _et al._ (2022),
non-standard neutrino interactions (NSI) Miranda  _et al._ (2020a); Barranco
_et al._ (2005); Scholberg (2006); Liao and Marfatia (2017); Giunti (2020);
Denton and Gehrlein (2021); Khan  _et al._ (2021); Coloma  _et al._ (2022);
Chatterjee  _et al._ (2022), generalized neutrino interactions Dutta  _et al._
(2016); Flores  _et al._ (2022), neutrino magnetic moments Miranda  _et al._
(2019a), non-unitarity of the leptonic mixing matrix Miranda  _et al._
(2020b), light sterile neutrinos Dutta  _et al._ (2016); Kosmas  _et al._
(2017); Miranda  _et al._ (2019b, 2020b), dipole portal Dasgupta  _et al._
(2021); Bolton  _et al._ (2022), electromagnetic properties of neutrinos
Cadeddu  _et al._ (2020b) etc. For a wider review the reader is referred to De
Romeri  _et al._ (2022); Cadeddu  _et al._ (2020b) . In this work, we focus on
leptoquark (LQ) models.

LQs are particles beyond the SM, carrying both lepton and baryon numbers, that
in general can be scalar or vector. Here we focus on the scalar LQs case which
is less constrained compared to the vector one Dey and Mohanty (2016);
Valencia and Willenbrock (1994); Smirnov (2007). Indeed vector LQs couplings
are fixed from the gauge structure and there is less freedom. They arise in a
large variety of beyond the SM scenarios, such as grand unified theories
Georgi and Glashow (1974); Fritzsch and Minkowski (1975), the Pati-Salam model
Pati and Salam (1974), extended technicolor models Farhi and Susskind (1981);
Gabriel  _et al._ (1994); Andersen  _et al._ (2011) etc. In general, the mass
of such LQs are expected to be close to the GUT scale to avoid proton decay.
There are however models with low scale LQs Buchmüller  _et al._ (1987);
Belyaev  _et al._ (2005); Dorsner and Fileviez Perez (2005); Fileviez Perez
and Wise (2013) which forbid rapid proton decays. We consider scalar LQs which
do not give rise to tree level interactions of the type ($qlqq$), so proton
decay does not occur.

In this paper, we utilize the COHERENT data of CsI to derive the exclusion
limits on scalar LQs via modification to the SM CE$\nu$NS event rate in the
presence of LQs. The bounds cover a wide mass range from MeV to TeV scales. We
also compare our results with existing constraints on LQs.

The structure of this paper is as follows: section II details our model and
notation and provides the relevant part of the lagrangian that mediates the
coherent process. In section III, we describe CE$\nu$NS SM and LQ signal
predictions. Section IV contains the $\chi^{2}$ analysis and the bounds on LQ
parameters from coherent and other measurements. Finally, section V summarizes
our conclusions.

##  II Models

| $SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$ | LQ
---|---|---
$(\overline{Q^{c}}_{L})^{i}L^{j}_{L}$ | ${\bf(3,3,-1/3)}$ | $S_{3}{\bf(\overline{3},3,1/3)}$
$(\overline{Q^{c}}_{L})^{i}L^{j}_{L}$ | ${\bf(3,1,-1/3)}$ | $S_{1}{\bf(\overline{3},1,1/3)}$
$(\overline{u^{c}}_{R})^{i}\ell_{R}^{j}$ | ${\bf(3,1,-1/3)}$ | $S_{1}{\bf(\overline{3},1,1/3)}$
$\overline{u}_{R}^{i}L^{j}_{L}$ | ${\bf(\overline{3},2,-7/6)}$ | $R_{2}{\bf(3,2,7/6)}$
$\overline{\ell}_{R}^{i}Q_{L}^{j}$ | ${\bf(\overline{3},2,-7/6)}$ | $R_{2}{\bf(3,2,7/6)}$
$\overline{d}_{R}^{i}L^{j}_{L}$ | ${\bf(\overline{3},2,-1/6)}$ | $\tilde{R}_{2}{\bf(3,2,1/6)}$
$(\overline{d^{c}}_{R})^{i}\ell_{R}^{j}$ | ${\bf(3,1,-4/3)}$ | $\tilde{S}_{1}{\bf(\overline{3},1,4/3)}$
Table 1:  First column corresponds to the list of $q\ell$ interaction terms.
The second column shows the quantum number of $q\ell$ terms. In the third
column, we report the LQ symbol and the corresponding quantum number.

LQs couple to quarks and leptons. Therefore, LQ quantum numbers are not
arbitrary and must be fixed accordingly. In Tab 1, we give the
$SU(2)_{L}\times U(1)_{Y}$ quantum numbers of all possible terms involving
quarks, leptons and a scalar LQ. However, the above table does not
exhaustively report the possible LQ interactions which also include “diquark”
type interactions which occur for $S_{3}$, $S_{1}$, and $\tilde{S}_{1}$. Two
scalar LQs with quantum numbers (3, 2, 1/6) and (3, 2, 7/6) i.e.
$\tilde{R}_{2}$ and $R_{2}$ do not contain any “diquark” type interactions
which cause the proton decay Doršner  _et al._ (2016). But in principle,
$\tilde{R}$ can generate the renormalizable operator
$H\tilde{R}\tilde{R}\tilde{R}$ via the SM Higgs (H). This operator leads to
proton decay ($p\rightarrow\pi^{+}+\pi^{-}+e^{-}+\nu\nu$) Arnold  _et al._
(2013). We can either neglect or forbid this operator by (some) symmetry. For
simplicity, we consider each LQ separately and do not analyze the scenarios
with two LQs simultaneously. The SM Lagrangian is minimally extended to
include the Yukawa interaction of the LQs with the quark and lepton fields.
The relevant part of the Lagrangian that contributes in the neutrino-nucleon
coherent process is given by

$\begin{split}\mathcal{L}^{\Delta_{1}}&\supset-
y_{ij}^{(1)}\bar{d}^{i}_{R}\tilde{R}_{2}L^{j}_{L}+h.c.\\\
&=-y_{ij}^{(1)}\bar{d}^{i}P_{L}\ell^{j}\Delta^{2/3}_{1}-y_{ij}^{(1)}\bar{d}^{i}P_{L}\nu^{j}\Delta^{-1/3}_{1}+h.c.\end{split}$
(1)

and

$\begin{split}\mathcal{L}^{\Delta_{2}}&\supset-
y_{ij}^{(2)}\bar{u}^{i}_{R}R_{2}L^{j}_{L}+h.c.\\\
&=-y_{ij}^{(2)}\bar{u}^{i}P_{L}\ell^{j}\Delta^{5/3}_{2}-y_{ij}^{(2)}\bar{u}^{i}P_{L}\nu^{j}\Delta^{2/3}_{2}+h.c.\,,\end{split}$
(2)

where $\tilde{R}_{2}=(\Delta^{2/3}_{1},\Delta^{-1/3}_{1})^{T}$ and
$R_{2}=(\Delta^{5/3}_{2},\Delta^{2/3}_{2})^{T}$ under $SU(2)_{L}$, $i,j$ are
generation indices, $d^{i}$ and $u^{i}$ are the down-type and up-type quark
fields, $l^{j}_{L}$ are the lepton doublets
$L^{j}_{L}=\begin{pmatrix}\nu^{j}&\ell^{j}\end{pmatrix}^{T}$, $P_{L,R}$, are
the usual left and right chiral projection operators and $y^{\prime}$s are the
Yukawa coupling matrices. The masses of LQs in the doublet structure could be
different. We can use this feature to avoid many constraints on the Yukawa
couplings. Since CE$\nu$NS involves quark mass eigenstates and neutrino
interaction eigenstates, the coupling matrices we are interested in are
$\tilde{y}^{(1)}=V_{d}^{\dagger}\,y^{(1)}$ and
$\tilde{y}^{(2)}=V_{u}^{\dagger}\,y^{(2)}$ where $V_{d}$ ($V_{u}$) are
arbitrary, unitary matrices which rotate the quark fields into their mass
basis. Therefore, since $y^{(1)}$ ($y^{(2)}$) are arbitrary matrices,
$\tilde{y}^{(1)}$ ($\tilde{y}^{(2)}$) are also completely arbitrary a priori.
To avoid having Flavour Changing Neutral Current problems at tree level, we
consider only couplings with the first generation of quarks. This choice does
not affect our prediction for the CE$\nu$NS since it only involves the valence
quarks, hence only the first generation.

In principle, one should also include terms like

$\mathcal{L}\supset
y_{ij}^{(3)}\overline{\ell}^{i}_{R}(R_{2}^{*})^{a}(Q_{L}^{j})^{a}+h.c.\,.$ (3)

Such terms do not contribute to CE$\nu$NS but can modify the pion decay. There
are strong constraint on LQ coming from pion decay if we consider
simultaneously both the Eq. 2 and Eq.3 type operators Davidson  _et al._
(1994). We assume the coupling in Eq.3 to be very small and explore the
COHERENT data to constrain the couplings in Eqns. 1 and 2.

##  III Signal Prediction in COHERENT

In this section, we present the signal prediction at COHERENT for both the SM
and LQs scenarios.

At the COHERENT experiment, the neutrino fluxes come from the Spallation
Neutron Source (SNS). Prompt pion decay
($\pi^{+}\rightarrow\mu^{+}+\nu_{\mu}$) at rest and the subsequently delayed
decay of muons ($\mu^{+}\rightarrow e^{+}+\nu_{e}+\bar{\nu}_{\mu}$) produce
three different neutrino fluxes. The $\nu_{\mu}$ reach the detector within
$\sim 1.5~{}\mu$s after protons-on-target while the $\bar{\nu}_{\mu}$ and
$\nu_{e}$ fluxes arrive in a comparatively longer time interval ($\sim
10~{}\mu$s). The differential fluxes for neutrinos are as follows :

$\displaystyle\dfrac{dN_{\nu_{\mu}}}{dE}=\eta\,\delta\,\Big{(}E-\dfrac{m^{2}_{\pi}-m^{2}_{\mu}}{2m_{\pi}}\Big{)}$
(4)
$\displaystyle\dfrac{dN_{\bar{\nu}_{\mu}}}{dE}=\eta\,\dfrac{64E^{2}}{m^{3}_{\mu}}\,\Big{(}\dfrac{3}{4}-\dfrac{E}{m_{\mu}}\Big{)}$
(5)
$\displaystyle\dfrac{dN_{\nu_{e}}}{dE}=\eta\,\dfrac{192E^{2}}{m^{3}_{\mu}}\,\Big{(}\dfrac{1}{2}-\dfrac{E}{m_{\mu}}\Big{)},$
(6)

where $m_{\mu}$ and $m_{\pi}$ are the muon and pion masses respectively,
$\eta=r\,\rm{N_{POT}}/4\pi L^{2}$ is a normalization factor, $\rm{N_{POT}}$ is
the total number of protons on target (POT), $r$ corresponds to neutrinos per
flavour produced for each POT, and L measures the distance between the source
and the detector. For the CsI detector r = 0.08, $N_{\rm{POT}}=17.6\times
10^{22}$ and L = 19.3 m.

The SM contribution to the CE$\nu$NS cross-section for a neutrino ($\nu_{l}$)
as a function of the nuclear recoil kinetic energy ($T_{nr}$) is given by

$\displaystyle\dfrac{d\sigma_{\nu_{l}-N}}{dT_{nr}}(E,T_{nr})=\dfrac{G^{2}_{F}M}{\pi}\Big{(}1-\dfrac{MT_{nr}}{2E^{2}}\Big{)}Q^{2}_{l,SM}.$
(7)

In the above, $l$ is the neutrino flavour index, $M$ is the mass of the
detector nucleus with $Z$ protons and $N$ neutrons, and $E$ is the energy of
incoming neutrino, $G_{F}$ is the Fermi constant, and
$Q^{2}_{l,SM}=[g^{p}_{V}(\nu_{l})ZF_{Z}(|\vec{q}|^{2})+g^{n}_{V}NF_{N}(|\vec{q}|^{2})]^{2}$.
Considering the radiative corrections in the minimal subtraction,
$\rm{\overline{MS}}$ scheme, the values of the $g^{p}_{V}$ and $g^{n}_{V}$
couplings are Erler and Su (2013); Cadeddu  _et al._ (2020b)

$\displaystyle g^{p}_{V}(\nu_{e})$ $\displaystyle=$ $\displaystyle 0.0401,$
(8) $\displaystyle g^{p}_{V}(\nu_{\mu})$ $\displaystyle=$ $\displaystyle
0.0318,$ (9) $\displaystyle g^{n}_{V}$ $\displaystyle=$
$\displaystyle-0.5094.$ (10)

The proton and neutron distributions in the nucleus are represented by the
form factors ($F_{Z}(|\vec{q}|^{2})$ and $F_{N}(|\vec{q}|^{2})$) which depend
on the three-momentum transfer $|\vec{q}|\simeq\sqrt{2MT_{nr}}$. We use the
Helm parameterization for the form factors Helm (1956)

$\displaystyle
F(|\vec{q}|^{2})=3\dfrac{j_{1}(|\vec{q}|R_{0})}{|\vec{q}|R_{0}}e^{-|\vec{q}|^{2}s^{2}/2},$
(11)

where $j_{1}(x)=\sin(x)/x^{2}-\cos(x)/x$ is the spherical Bessel function and
$s=0.9$ fm, and the rms radius $R$ is related to $R_{0}$ by
$R^{2}_{p,n}=3R^{2}_{0}/5+3s^{2}$. Helm’s form factors are effectively
identical to the Fermi parameterization for our purposes Giunti (2020). We
take $R_{p}(\rm{Cs})=4.804\,\rm{fm}$ and $R_{p}(\rm{I})=4.749\,\rm{fm}$ as the
rms radius of the proton Fricke  _et al._ (1995); Angeli and Marinova (2013).
The neutron rms radius is not known with good accuracy for CsI. In our
analysis, we take the values $R_{n}(Cs)=5.01$ fm and $R_{n}(I)=4.94$ fm Bender
_et al._ (1999). If the transfer of momentum is such that
$|\vec{q}|R_{p}\gtrsim 1$ and $|\vec{q}|R_{n}\gtrsim 1$, the scattering loses
coherence.

The number of events in the i-th bin of nuclear-recoil energy is given by

$\displaystyle
N_{i}^{\mathrm{CE}\nu\mathrm{NS}}=N(\mathcal{N})\int^{T_{\mathrm{nr}}^{i+1}}_{T_{\mathrm{nr}}^{i}}dT_{\mathrm{nr}}\,A(T_{\mathrm{nr}})\int^{E_{\mathrm{max}}}_{E_{\mathrm{min}}}dE$
$\displaystyle\sum_{\nu=\nu_{e},\nu_{\mu},\overline{\nu}_{\mu}}\frac{dN_{\nu}}{dE}\dfrac{d\sigma_{\nu\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}}(E,T_{\mathrm{nr}})\,,$
(12)

where $E_{\mathrm{min}}=\sqrt{MT_{\mathrm{nr}}/2}$,
$E_{\mathrm{max}}=m_{\mu}/2\sim 52.8$ MeV, and $A(T_{\mathrm{nr}})$ is the
energy-dependent reconstruction efficiency. The nuclear recoil energy
$T_{\mathrm{nr}}\,[\mathrm{keV}]$ is converted into the corresponding electron
recoil energy $T_{ee}\,[\mathrm{keV}]$

$T_{ee}=(1-\alpha_{3})f_{Q}(T_{\mathrm{nr}})T_{\mathrm{nr}}\,,$ (13)

where $f_{Q}$ is the quenching factor and $\alpha_{3}$ is a nuisance parameter
quantifying the uncertainty on $f_{Q}$. We take the CsI quenching factor from
Collar  _et al._ (2019). $N(\mathcal{N})$ is the total number of atoms present
in the active detector mass ($M_{\rm det}$) and
$N(\mathcal{N})=N_{\mathrm{A}}\,M_{\mathrm{det}}/M_{\mathcal{N}}$, where
$N_{\mathrm{A}}$ is Avogadro’s number and $M_{\mathcal{N}}$ is the molar mass
of the detector material. For CsI we take $M_{\mathrm{det}}=14.6$ kg and
$M_{\mathrm{CsI}}=259.8$ g/mol.

In the presence of LQs the standard CE$\nu$NS cross section will be modified.
Processes contributing to the CE$\nu$NS are shown Fig. 1. In this figure we
can see that we have vertices with a neutrino, LQ, and a quark. As the typical
momentum transfer to the nucleus is the order of a few tens of MeV, these
diagrams can be represented by the effective four fermion interaction for
$m_{\Delta}\gtrapprox 10$ MeV like

$\displaystyle\mathcal{L}^{\Delta}_{\rm{eff}}=\dfrac{y^{2}}{m^{2}_{\Delta}}\,(\bar{\psi}_{N}P_{L}\nu)\,(\bar{\nu}P_{R}\psi_{N})\,.$
(14)

where $\psi_{N}$ is either a $u$ or $d$ quark field. However, we need to
factorize this diagram into a neutrino current and a hadronic current to
compute the CE$\nu$NS events. This can be achieved through Fierz
transformations giving the following effective interaction

$\displaystyle\mathcal{L}^{\Delta}_{\rm{eff}}\sim-\dfrac{y^{2}}{2m^{2}_{\Delta}}\,(\bar{\psi}_{N}\gamma^{\mu}P_{R}\psi_{N})\,(\bar{\nu}\gamma_{\mu}P_{L}\nu).$
(15)

Figure 1:  Tree level processes involving the Leptoquark which modifies the SM
prediction for the CE$\nu$NS cross-section. The final state neutral lepton
flavor is not experimentally distinguishable.

Writing the effective interaction in this way, with a $V-A$ axial structure,
allows us to write down the CE$\nu$NS cross section including the contribution
from LQs

$\frac{d\sigma_{\nu_{i}-N}}{dT_{nr}}(E,T_{nr})=\dfrac{G^{2}_{F}M}{\pi}\left(1-\dfrac{MT_{nr}}{2E^{2}}\right)Q^{2}_{i,\,k}$
(16)
$Q^{2}_{i,\,k}=\left(\left(Q_{i,\,SM}+Q_{ii,\,\Delta_{k}}\right)^{2}+\sum_{j\neq
i}Q^{2}_{ij,\,\Delta_{k}}\right)$ (17)

where

$\begin{split}&Q_{ij,\,\Delta_{1}}=\frac{\tilde{y}_{1i}^{(1)}\tilde{y}_{1j}^{(1)}}{4\sqrt{2}\,G_{F}}\,\frac{ZF_{Z}(|\vec{q}|^{2})+2\,NF_{N}(|\vec{q}|^{2})}{|\vec{q}|^{2}+m_{\Delta_{1}^{-1/3}}^{2}},\\\
&Q_{ij,\,\Delta_{2}}=\frac{\tilde{y}_{1i}^{(2)}\tilde{y}_{1j}^{(2)}}{4\sqrt{2}\,G_{F}}\,\frac{2\,ZF_{Z}(|\vec{q}|^{2})+NF_{N}(|\vec{q}|^{2})}{|\vec{q}|^{2}+m_{\Delta_{2}^{2/3}}^{2}}.\end{split}$
(18)

Figure 2:  The plot shows the number of events vs the photo-electron number.
The experimental CE$\nu$NS events minus the background are presneted in cyan.
The dashed black line represents the SM event prediction, the purple solid
line corresponds to the prediction including $\Delta_{1}^{-1/3}$, the orange
solid line is the prediction for $\Delta_{2}^{2/3}$. The theoretical
predictions for the case A (see Tab. 2) are obtained assuming
$m_{\Delta}=1$GeV and $g_{\Delta}=4\cdot 10^{-3}$.

Note that typical $U^{\prime}(1)$ based models Cadeddu  _et al._ (2021b)
consider the processes like $N+\nu_{i}\to N+\nu_{i}$ for CE$\nu$NS. On the
contrary here we have $N+\nu_{i}\to N+\nu_{j}$ where $i$ can be different from
$j$ as can be seen from Fig. 1.

LQ Case | A | B | C
---|---|---|---
$\Delta_{1}^{-1/3}$ | $\tilde{y}^{(1)}=\left(\begin{array}[]{ccc}g_{\Delta_{1}}&g_{\Delta_{1}}&0\\\ 0&0&0\\\ 0&0&0\end{array}\right)$ | $\tilde{y}^{(1)}=\left(\begin{array}[]{ccc}0&g_{\Delta_{1}}&0\\\ 0&0&0\\\ 0&0&0\end{array}\right)$ | $\tilde{y}^{(1)}=\left(\begin{array}[]{ccc}g_{\Delta_{1}}&0&0\\\ 0&0&0\\\ 0&0&0\end{array}\right)$
$\Delta_{2}^{2/3}$ | $\tilde{y}^{(2)}=\left(\begin{array}[]{ccc}g_{\Delta_{2}}&g_{\Delta_{2}}&0\\\ 0&0&0\\\ 0&0&0\end{array}\right)$ | $\tilde{y}^{(2)}=\left(\begin{array}[]{ccc}0&g_{\Delta_{2}}&0\\\ 0&0&0\\\ 0&0&0\end{array}\right)$ | $\tilde{y}^{(2)}=\left(\begin{array}[]{ccc}g_{\Delta_{2}}&0&0\\\ 0&0&0\\\ 0&0&0\end{array}\right)$
Table 2:  In the table we report the benchmark cases studied in this paper.

##  IV Results

Figure 3:  Constraints obtained in the plane $g_{\Delta}-m_{\Delta}$. In the
left panel are shown the constraints on $\Delta^{-1/3}_{1}$, while on the
right we have the constraints on $\Delta^{2/3}_{2}$. In solid colored lines
are reported the constraints at 90% C.L. obtained with our analysis. In dashed
colored lines are reported the constraints valid in our scenario from CMS (in
cyan) Tumasyan  _et al._ (2021) and IceCube (in violet) Dey and Mohanty
(2016). The dashed grey lines are constraints reported for comparison but are
not valid in our scenarios. Such constraints are taken from Ref. Taxil  _et
al._ (2000).

Using Eqns. (7, 16, 17, 18), we compute the theoretical prediction for the
event number of CE$\nu$NS assuming only the SM or assuming the existence of
either $\Delta_{1}$ or $\Delta_{2}$. As mentioned above, we are considering
the LQs that couple only to the first generation of quarks since CE$\nu$NS is
sensitive only to valence quarks. Moreover, this prevents us from having FCNC
at tree level mediated by LQs. For each LQ, we consider three possible
scenarios. $\Delta_{k}$ couples to both $\nu_{e}$ and $\nu_{\mu}$ (case A) or
$\Delta_{k}$ couples only to one of $\nu_{\mu}$ and $\nu_{e}$ (Case B and C)
as summarized in Tab. 2. It is possible to consider other possibilities, but
they are less conservative.

To compare our predictions with the measured event rate, we need to write our
prediction as a function of the photo-electron number instead of the nuclear
recoil energy. In Ref. Akimov  _et al._ (2018), it is shown that

$N_{PE}=Y_{T}T_{ee}\,,$ (19)

where $Y_{T}=13.35\,N_{PE}/keV$ is the light yield of the phototubes. Eq. III
can be rewritten in terms of $N_{PE}$ as

$\displaystyle
N_{i}^{\mathrm{CE}\nu\mathrm{NS}}=N(\mathcal{N})\int^{N_{PE}^{i+1}}_{N_{PE}^{i}}dN_{PE}\,f(N_{PE})\int^{E_{\mathrm{max}}}_{E_{\mathrm{min}}}dE$
$\displaystyle\sum_{\nu=\nu_{e},\nu_{\mu},\overline{\nu}_{\mu}}\frac{dN_{\nu}}{dE}\dfrac{d\sigma_{\nu\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}}(E,T_{\mathrm{nr}})\,,$
(20)

where $f(N_{PE})$ is the detector efficiency in terms of the photo-electron
content defined as Akimov  _et al._ (2018)

$f(x)=\frac{a}{1+\exp(-k(x-x_{0}))}\Theta(x)\,,$ (21)

where

$\begin{split}&a=0.6655^{+0.0212}_{-0.0384}\,,\\\
&k=0.4942^{+0.0335}_{-0.0131}\,,\\\ &x_{0}=10.8507^{+0.1838}_{-0.3995}\,.\\\
\end{split}$ (22)

The function $\Theta(x)$ is defined as

$\Theta(x)=\left\\{\begin{array}[]{cc}0&x<5\\\ 0.5&5\leq x\leq 6\\\
1&x>6\end{array}\right.\,.$ (23)

In Fig. 2, we show the events for CsI in COHERENT. The cyan band corresponds
to the experimental minus the background CE$\nu$NS events in the range
$N_{PE}\in[6-30]$ and the SM predictions are depicted by a dashed black line.
Assuming Case A (see Tab. 2), we report the events for $\Delta_{1}$ and
$\Delta_{2}$ with a solid purple line and with a solid orange line
respectively.

To constrain the models, we perform a $\chi^{2}$ analysis analogous to that in
Ref. Cadeddu  _et al._ (2021b)

$\begin{split}\chi^{2}=\min_{\boldsymbol{\alpha}}\sum_{i=4}^{15}&\left[\left(\frac{N_{i}^{\rm
exp}-(1+\alpha_{1})N_{i}^{CE\nu
NS}-(1+\alpha_{2})B_{i}}{\sigma_{i}}\right)^{2}\right.\\\
&\left.+\sum_{j}\left(\frac{\alpha_{j}}{\sigma_{\alpha_{j}}}\right)^{2}\right]\,.\end{split}$
(24)

For each bin $i$, we indicate with $N_{i}^{\rm exp}$ the measured event number
reported in Ref. Akimov  _et al._ (2018), $N_{i}^{CE\nu NS}$ is the expected
event rate, $B_{i}$ is the background event rate, and $\sigma_{i}$ is the
statistical uncertainty taken to be equal to the square root of the event
number. As done in Refs. Cadeddu  _et al._ (2020a, 2021b), we consider only
the energy bins from $i=4$ to $i=15$ because they span the recoil kinetic
energy range covered by the most recent quenching factor measurement performed
at Chicago-3 Collar  _et al._ (2019). We are considering 3 nuisance parameters
$\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},\alpha_{3})$ which quantify the
systematic uncertainty of the signal rate, background rate, and quenching
factor respectively. The corresponding standard deviations are
$\sigma_{\boldsymbol{\alpha}}=(0.112,0.25,0.051)$ Akimov  _et al._ (2017). The
obtained constraints on $\Delta_{1}^{-1/3}$ and on $\Delta_{2}^{2/3}$ are
reported on the left and right side of Fig. 3 as continuous lines. In the
lower side of the same figure, we place a zoomed version of the upper figure
in the range $m_{\Delta_{k}}\in[10^{2}-10^{4}]$ GeV for clarity. In each plot,
we report existing constraints on LQs for comparison. We report as dashed
lines the constraints that are valid in our case, and with dotted lines the
ones that are not directly applicable to our scenarios. The cyan dashed line
corresponds to constraints from CMS taken from Ref. Tumasyan  _et al._ (2021).
These constraints are obtained assuming that $BR(\Delta\to q+\nu)=1$ and
assume that $\Delta$ is pair produced via gluon fusion. The violet dashed line
refers to constraints obtained using IceCube data and were taken from Ref. Dey
and Mohanty (2016). This is not the only constraint coming from IceCube (see
Refs. Dey  _et al._ (2018); Bečirević  _et al._ (2018)). The dark and light
dotted lines refer to LEP and Deep Inelastic Scattering constraints taken from
Ref. Taxil  _et al._ (2000). These constraints do not apply to our scenario
since they constrain the element of the LQ doublet which couples to the
charged leptons.

The choices we made for the Yukawa couplings prevent us from having
constraints from FCNC and B-physics. Moreover, we can overcome the constraints
from LFV considering a mass splitting between the member of the doublets.

##  V Conclusions

In this work, we utilize the neutrino-nucleus coherent elastic scattering data
to put bounds on scalar LQs over a wide mass range from MeV to TeV scales. We
consider scalar LQs, in particular, with hypercharge Y=1/6 and Y=7/6. These
LQs do not have any “diquark” couplings which contribute to the proton decay.
Hence effectively the LQ masses could be well below the GUT scale. We use the
Fierz transformation to recast the effective four-fermion interaction vertex
into neutrino-neutrino and quark-quark currents. This allows us to calculate
the hadronic matrix elements for the coherent process. Only the valence quarks
(u and d) of the nuclei participate in neutrino-nucleus coherent scattering.
Therefore, we can constrain the specific structure of the Yukawa matrix that
couples to the first generation of quarks, and neutrinos. These choices of
couplings evade bounds from FCNCs and LFV. Hence, COHERENT measurements
provide distinctive quality to probe this particular Yukawa structure. We also
show the constraints coming from LHC, IceCube, LEP, and DIS. The constraint is
comparable with other experiments for some mass regions. However, putting
constraints on LQs over a wide mass range is the unique capability of coherent
neutrino-nucleus scattering measurements.

Acknowledgments

We thank Carlo Giunti, Ujjal Kumar Dey, and Marco Chianese for the useful
discussion.

This work was partially supported by the research grant number 2017W4HA7S
“NAT-NET: Neutrino and Astroparticle Theory Network” under the program PRIN
2017 funded by the Italian Ministero dell’Università e della Ricerca (MUR).

## References

  * Akimov  _et al._ (2017) D. Akimov _et al._ (COHERENT), Science 357, 1123 (2017), arXiv:1708.01294 [nucl-ex] .
  * Freedman (1974) D. Z. Freedman, Phys. Rev. D 9, 1389 (1974).
  * Akimov  _et al._ (2021) D. Akimov _et al._ (COHERENT), Phys. Rev. Lett. 126, 012002 (2021), arXiv:2003.10630 [nucl-ex] .
  * Miranda  _et al._ (2020a) O. G. Miranda, D. K. Papoulias, G. Sanchez Garcia, O. Sanders, M. Tórtola, and J. W. F. Valle, JHEP 05, 130 (2020a), [Erratum: JHEP 01, 067 (2021)], arXiv:2003.12050 [hep-ph] .
  * Cadeddu  _et al._ (2020a) M. Cadeddu, F. Dordei, C. Giunti, Y. F. Li, and Y. Y. Zhang, Phys. Rev. D 101, 033004 (2020a), arXiv:1908.06045 [hep-ph] .
  * Cadeddu  _et al._ (2021a) M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, E. Picciau, C. A. Ternes, and Y. Y. Zhang, Phys. Rev. C 104, 065502 (2021a), arXiv:2102.06153 [hep-ph] .
  * Cadeddu  _et al._ (2018) M. Cadeddu, C. Giunti, Y. F. Li, and Y. Y. Zhang, Phys. Rev. Lett. 120, 072501 (2018), arXiv:1710.02730 [hep-ph] .
  * Canas  _et al._ (2020) B. C. Canas, E. A. Garces, O. G. Miranda, A. Parada, and G. Sanchez Garcia, Phys. Rev. D 101, 035012 (2020), arXiv:1911.09831 [hep-ph] .
  * Cadeddu  _et al._ (2021b) M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, E. Picciau, and Y. Y. Zhang, JHEP 01, 116 (2021b), arXiv:2008.05022 [hep-ph] .
  * Flores  _et al._ (2020) L. J. Flores, N. Nath, and E. Peinado, JHEP 06, 045 (2020), arXiv:2002.12342 [hep-ph] .
  * Banerjee  _et al._ (2021) H. Banerjee, B. Dutta, and S. Roy, Phys. Rev. D 104, 015015 (2021), arXiv:2103.10196 [hep-ph] .
  * de la Vega  _et al._ (2021) L. M. G. de la Vega, L. J. Flores, N. Nath, and E. Peinado, JHEP 09, 146 (2021), arXiv:2107.04037 [hep-ph] .
  * Atzori Corona  _et al._ (2022) M. Atzori Corona, M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, E. Picciau, C. A. Ternes, and Y. Y. Zhang, JHEP 05, 109 (2022), arXiv:2202.11002 [hep-ph] .
  * Binh  _et al._ (2021) D. T. Binh, L. T. Hue, V. H. Binh, D. V. Soa, and H. N. Long, (2021), arXiv:2109.08118 [hep-ph] .
  * Bertuzzo  _et al._ (2022) E. Bertuzzo, G. Grilli di Cortona, and L. M. D. Ramos, JHEP 06, 075 (2022), arXiv:2112.04020 [hep-ph] .
  * Chakraborty  _et al._ (2022) K. Chakraborty, A. Das, S. Goswami, and S. Roy, JHEP 04, 008 (2022), arXiv:2111.08767 [hep-ph] .
  * De Romeri  _et al._ (2022) V. De Romeri, O. G. Miranda, D. K. Papoulias, G. Sanchez Garcia, M. Tórtola, and J. W. F. Valle, (2022), arXiv:2211.11905 [hep-ph] .
  * Barranco  _et al._ (2005) J. Barranco, O. G. Miranda, and T. I. Rashba, JHEP 12, 021 (2005), arXiv:hep-ph/0508299 .
  * Scholberg (2006) K. Scholberg, Phys. Rev. D 73, 033005 (2006), arXiv:hep-ex/0511042 .
  * Liao and Marfatia (2017) J. Liao and D. Marfatia, Phys. Lett. B 775, 54 (2017), arXiv:1708.04255 [hep-ph] .
  * Giunti (2020) C. Giunti, Phys. Rev. D 101, 035039 (2020), arXiv:1909.00466 [hep-ph] .
  * Denton and Gehrlein (2021) P. B. Denton and J. Gehrlein, JHEP 04, 266 (2021), arXiv:2008.06062 [hep-ph] .
  * Khan  _et al._ (2021) A. N. Khan, D. W. McKay, and W. Rodejohann, Phys. Rev. D 104, 015019 (2021), arXiv:2104.00425 [hep-ph] .
  * Coloma  _et al._ (2022) P. Coloma, I. Esteban, M. C. Gonzalez-Garcia, L. Larizgoitia, F. Monrabal, and S. Palomares-Ruiz, JHEP 05, 037 (2022), arXiv:2202.10829 [hep-ph] .
  * Chatterjee  _et al._ (2022) S. S. Chatterjee, S. Lavignac, O. G. Miranda, and G. Sanchez Garcia, (2022), arXiv:2208.11771 [hep-ph] .
  * Dutta  _et al._ (2016) B. Dutta, Y. Gao, R. Mahapatra, N. Mirabolfathi, L. E. Strigari, and J. W. Walker, Phys. Rev. D 94, 093002 (2016), arXiv:1511.02834 [hep-ph] .
  * Flores  _et al._ (2022) L. J. Flores, N. Nath, and E. Peinado, Phys. Rev. D 105, 055010 (2022), arXiv:2112.05103 [hep-ph] .
  * Miranda  _et al._ (2019a) O. G. Miranda, D. K. Papoulias, M. Tórtola, and J. W. F. Valle, JHEP 07, 103 (2019a), arXiv:1905.03750 [hep-ph] .
  * Miranda  _et al._ (2020b) O. G. Miranda, D. K. Papoulias, O. Sanders, M. Tórtola, and J. W. F. Valle, Phys. Rev. D 102, 113014 (2020b), arXiv:2008.02759 [hep-ph] .
  * Kosmas  _et al._ (2017) T. S. Kosmas, D. K. Papoulias, M. Tortola, and J. W. F. Valle, Phys. Rev. D 96, 063013 (2017), arXiv:1703.00054 [hep-ph] .
  * Miranda  _et al._ (2019b) O. G. Miranda, G. Sanchez Garcia, and O. Sanders, Adv. High Energy Phys. 2019, 3902819 (2019b), arXiv:1902.09036 [hep-ph] .
  * Dasgupta  _et al._ (2021) A. Dasgupta, S. K. Kang, and J. E. Kim, JHEP 11, 120 (2021), arXiv:2108.12998 [hep-ph] .
  * Bolton  _et al._ (2022) P. D. Bolton, F. F. Deppisch, K. Fridell, J. Harz, C. Hati, and S. Kulkarni, Phys. Rev. D 106, 035036 (2022), arXiv:2110.02233 [hep-ph] .
  * Cadeddu  _et al._ (2020b) M. Cadeddu, F. Dordei, C. Giunti, Y. F. Li, E. Picciau, and Y. Y. Zhang, Phys. Rev. D 102, 015030 (2020b), arXiv:2005.01645 [hep-ph] .
  * Dey and Mohanty (2016) U. K. Dey and S. Mohanty, JHEP 04, 187 (2016), arXiv:1505.01037 [hep-ph] .
  * Valencia and Willenbrock (1994) G. Valencia and S. Willenbrock, Phys. Rev. D 50, 6843 (1994), arXiv:hep-ph/9409201 .
  * Smirnov (2007) A. D. Smirnov, Mod. Phys. Lett. A 22, 2353 (2007), arXiv:0705.0308 [hep-ph] .
  * Georgi and Glashow (1974) H. Georgi and S. L. Glashow, Phys. Rev. Lett. 32, 438 (1974).
  * Fritzsch and Minkowski (1975) H. Fritzsch and P. Minkowski, Annals of Physics 93, 193 (1975).
  * Pati and Salam (1974) J. C. Pati and A. Salam, Phys. Rev. D 10, 275 (1974).
  * Farhi and Susskind (1981) E. Farhi and L. Susskind, Physics Reports 74, 277 (1981).
  * Gabriel  _et al._ (1994) T. Gabriel, D. Groom, P. Job, N. Mokhov, and G. Stevenson, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 338, 336 (1994).
  * Andersen  _et al._ (2011) J. R. Andersen _et al._ , Eur. Phys. J. Plus 126, 81 (2011), arXiv:1104.1255 [hep-ph] .
  * Buchmüller  _et al._ (1987) W. Buchmüller, R. Rückl, and D. Wyler, Physics Letters B 191, 442 (1987).
  * Belyaev  _et al._ (2005) A. Belyaev, C. Leroy, R. Mehdiyev, and A. Pukhov, JHEP 09, 005 (2005), arXiv:hep-ph/0502067 .
  * Dorsner and Fileviez Perez (2005) I. Dorsner and P. Fileviez Perez, Nucl. Phys. B 723, 53 (2005), arXiv:hep-ph/0504276 .
  * Fileviez Perez and Wise (2013) P. Fileviez Perez and M. B. Wise, Phys. Rev. D 88, 057703 (2013), arXiv:1307.6213 [hep-ph] .
  * Doršner  _et al._ (2016) I. Doršner, S. Fajfer, A. Greljo, J. F. Kamenik, and N. Košnik, Phys. Rept. 641, 1 (2016), arXiv:1603.04993 [hep-ph] .
  * Arnold  _et al._ (2013) J. M. Arnold, B. Fornal, and M. B. Wise, Phys. Rev. D 87, 075004 (2013), arXiv:1212.4556 [hep-ph] .
  * Davidson  _et al._ (1994) S. Davidson, D. C. Bailey, and B. A. Campbell, Z. Phys. C 61, 613 (1994), arXiv:hep-ph/9309310 .
  * Erler and Su (2013) J. Erler and S. Su, Prog. Part. Nucl. Phys. 71, 119 (2013), arXiv:1303.5522 [hep-ph] .
  * Helm (1956) R. H. Helm, Phys. Rev. 104, 1466 (1956).
  * Fricke  _et al._ (1995) G. Fricke, C. Bernhardt, K. Heilig, L. A. Schaller, L. Schellenberg, E. B. Shera, and C. W. de Jager, Atom. Data Nucl. Data Tabl. 60, 177 (1995).
  * Angeli and Marinova (2013) I. Angeli and K. P. Marinova, Atom. Data Nucl. Data Tabl. 99, 69 (2013).
  * Bender  _et al._ (1999) M. Bender, K. Rutz, P. G. Reinhard, J. A. Maruhn, and W. Greiner, Phys. Rev. C 60, 034304 (1999), arXiv:nucl-th/9906030 .
  * Collar  _et al._ (2019) J. I. Collar, A. R. L. Kavner, and C. M. Lewis, Phys. Rev. D 100, 033003 (2019), arXiv:1907.04828 [nucl-ex] .
  * Tumasyan  _et al._ (2021) A. Tumasyan _et al._ (CMS), JHEP 11, 153 (2021), arXiv:2107.13021 [hep-ex] .
  * Taxil  _et al._ (2000) P. Taxil, E. Tugcu, and J. M. Virey, Eur. Phys. J. C 14, 165 (2000), arXiv:hep-ph/9912272 .
  * Akimov  _et al._ (2018) D. Akimov _et al._ (COHERENT),  (2018), 10.5281/zenodo.1228631, arXiv:1804.09459 [nucl-ex] .
  * Dey  _et al._ (2018) U. K. Dey, D. Kar, M. Mitra, M. Spannowsky, and A. C. Vincent, Phys. Rev. D 98, 035014 (2018), arXiv:1709.02009 [hep-ph] .
  * Bečirević  _et al._ (2018) D. Bečirević, B. Panes, O. Sumensari, and R. Zukanovich Funchal, JHEP 06, 032 (2018), arXiv:1803.10112 [hep-ph] .


#  $0^{+}$ fully-charmed tetraquark states

Jian-Rong Zhang  Department of Physics, College of Liberal Arts and Sciences,
National University of Defense Technology, Changsha 410073, Hunan, People’s
Republic of China

###### Abstract

Motivated by the LHCb’s new observation of structures in the $J/\psi$-pair
invariant mass spectrum, for which could be classified as possible
$cc\bar{c}\bar{c}$ tetraquark candidates, we systematically study $0^{+}$
fully-charmed tetraquark states through QCD sum rules. Making the development
of calculation techniques to fourfold heavy hadronic systems, four different
configuration currents with $0^{+}$ are considered and vacuum condensates up
to dimension $6$ are included in the operator product expansion (OPE).
Finally, mass values acquired for $0^{+}$ $cc\bar{c}\bar{c}$ tetraquark states
agree well with the experimental data of the broad structure, which supports
that it could be a $0^{+}$ fully-charmed tetraquark state.

###### pacs:

11.55.Hx, 12.38.Lg, 12.39.Mk

##  I Introduction

By far, the topic of fully-charmed tetraquark state has attracted much
attention. For example, a variety of phenomenological models were employed to
predict the existence of some states merely made up of four heavy quarks
Predict-1 ; Predict-2 ; Predict-3 ; Predict-4 ; Predict-5 ; Predict-6 ;
Predict-7 ; Predict-8 ; Predict-9 ; Predict-10 ; Predict-11 ; Predict-12 ;
Predict-13 ; Predict-14 ; Predict-15 ; Predict-16 ; Predict-16-a ; Predict-17
; Predict-18 ; Predict-19 ; Predict-20 ; Predict-21 ; Predict-22 ; Predict-23
; Predict-24 ; Predict-25 ; Predict-26 ; Predict-27 . Particularly, without
any light quark contamination, fully-charmed tetraquark states are ideal
prototypes to refine one’s understanding on heavy quark dynamics.

Recently, the invariant mass spectrum of double-$J/\psi$ was researched using
proton-proton collision data recorded by the LHCb experiment, which shows a
broad structure just above twice the $J/\psi$ mass ranging from $6.2$ to
$6.8~{}\mbox{GeV}$ and a narrower structure around $6.9~{}\mbox{GeV}$,
referred to as $X(6900)$ LHCb . Soon after the LHCb’s new results, various
investigations were presented to explain them via different approaches Th-1 ;
Th-2 ; Th-2-a ; Th-2-aa ; Th-3 ; Th-4 ; Th-5 ; Th-6 ; Th-7 ; Th-8 ; Th-9 ;
Th-10 ; Th-11 ; Th-12 ; Th-13 ; Th-14 ; Th-15 ; Th-16 ; Th-17 ; Th-18 ; Th-19
; Th-20 ; Th-21 ; Th-22 ; Th-23 ; Th-24 ; Th-25 ; Th-26 . To probe a real
hadron, one inevitably has to face the sophisticated nonperturbative QCD
problem. As one trustable way for evaluating nonperturbative effects, the QCD
sum rule svzsum  is firmly established on the basic theory, and has been
widely applied to hadronic systems (for reviews see Refs. overview ; overview1
; overview2 ; overview3  and references therein). In particular on $0^{+}$
$cc\bar{c}\bar{c}$ tetraquark states, there are some existing works Predict-11
; Predict-12 ; Th-2 ; Th-10 ; Th-19  in different versions of QCD sum rules.

By comparison, Ref. Predict-11  used a moment QCD sum rule method augmented by
fundamental inequalities to explore the doubly hidden-charm/bottom tetraquark
states; Working with the Finite Energy version of the QCD Inverse Laplace sum
rules, Ref. Th-10  investigated doubly-hidden scalar heavy molecules and
tetraquarks states; Besides, there have some other QCD sum rule analysis of
fully-heavy tetraquark states involving condensate contributions up to
dimension $4$ in the OPE, specially choosing the axial vector-axial vector
configuration current in Ref. Predict-12 , paying attention to the scalar’s
first radial excited states Th-2 , or introducing a relative P-wave to the
diquark operator of the tetraquark current Th-19 . Consumingly motivated by
the so exciting and significant structures observed in the di-$J/\psi$ mass
spectrum, we would follow the previous QCD sum rule studies Nielsen ; Nielsen1
; Zhang  on hadrons containing one or two heavy quarks and devote to
developing the corresponding calculation techniques to fourfold heavy hadronic
systems. And then, we intend to systematically study $0^{+}$ fully-charmed
tetraquark states with QCD sum rules, by taking into account four possible
configuration currents and calculating condensates up to dimension $6$.

The paper is organized as follows. After the Introduction, the QCD sum rule is
derived for $0^{+}$ fully-charmed tetraquark states in Sec. II, along with
numerical analysis and discussions in Sec. III. The last part contains a brief
summary.

##  II $0^{+}$ fully-charmed tetraquark state QCD sum rules

Considering a tetraquark state, its interpolating current can ordinarily be
represented by a diquark-antidiquark configuration. Thus, following forms of
currents could be constructed for $0^{+}$ $cc\bar{c}\bar{c}$ tetraquark
states, with

$\displaystyle j$ $\displaystyle=$
$\displaystyle(Q_{a}^{T}C\gamma_{5}Q_{b})(\bar{Q}_{a}\gamma_{5}C\bar{Q}_{b}^{T})$

for the scalar-scalar configuration,

$\displaystyle j$ $\displaystyle=$
$\displaystyle(Q_{a}^{T}CQ_{b})(\bar{Q}_{a}C\bar{Q}_{b}^{T})$

for the pseudoscalar-pseudoscalar one,

$\displaystyle j$ $\displaystyle=$
$\displaystyle(Q_{a}^{T}C\gamma_{\mu}Q_{b})(\bar{Q}_{a}\gamma^{\mu}C\bar{Q}_{b}^{T})$

for the axial vector-axial vector (shortened to axial-axial) one, and

$\displaystyle j$ $\displaystyle=$
$\displaystyle(Q_{a}^{T}C\gamma_{5}\gamma_{\mu}Q_{b})(\bar{Q}_{a}\gamma^{\mu}\gamma_{5}C\bar{Q}_{b}^{T})$

for the vector-vector one. Here the index $T$ indicates matrix transposition,
$C$ means the charge conjugation matrix, $Q$ is the heavy charm quark, as well
as $a$ and $b$ are color indices. It is needed to state that these overall
scalar forms of currents are constructed as existing works mainly taking into
consideration that all of their Lorentz indices being contracted. One should
note that the situation for axial-axial and vector-vector could be more
complicated while characterizing a scalar state. In a general way, the
combination of two objects with spin $1$ may have total spin equal to $0$, $1$
and $2$, and the currents for axial-axial and vector-vector may represent a
mixture of spins. To select the spin $0$ part of these two currents, one could
try to project their correlation functions into the spin $0$ state with the
help of the appropriate projection operators, for which have been introduced
in Refs. remark ; remark1 .

To derive QCD sum rules, one can start with the two-point correlator

$\displaystyle\Pi(q^{2})=i\int d^{4}x\mbox{e}^{iq.x}\langle
0|T[j(x)j^{{\dagger}}(0)]|0\rangle.$ (1)

In phenomenology, it can be expressed as

$\displaystyle\Pi(q^{2})=\frac{\lambda_{H}^{2}}{M_{H}^{2}-q^{2}}+\frac{1}{\pi}\int_{s_{0}}^{\infty}\frac{\mbox{Im}\big{[}\Pi^{\mbox{phen}}(s)\big{]}}{s-q^{2}}ds,$
(2)

in which $M_{H}$ is the hadron’s mass, $s_{0}$ denotes the continuum
threshold, and $\lambda_{H}$ displays the coupling of the current to the
hadron $\langle 0|j|H\rangle=\lambda_{H}$. Moreover, it can theoretically be
rewritten as

$\displaystyle\Pi(q^{2})=\int_{(4m_{Q})^{2}}^{\infty}\frac{\rho}{s-q^{2}}ds,$
(3)

where $m_{Q}$ is the heavy charm mass, and the spectral density
$\rho=\frac{1}{\pi}\mbox{Im}\big{[}\Pi(s)\big{]}$. After equating Eqs. (2) and
(3), adopting quark-hadron duality, and applying a Borel transform, it yields

$\displaystyle\lambda_{H}^{2}e^{-M_{H}^{2}/M^{2}}$ $\displaystyle=$
$\displaystyle\int_{(4m_{Q})^{2}}^{s_{0}}\rho e^{-s/M^{2}}ds,$ (4)

with $M^{2}$ the Borel parameter. Taking the derivative of Eq. (4) with
respect to $-\frac{1}{M^{2}}$ and then dividing the result by Eq. (4) itself,
one could achieve the mass sum rule

$\displaystyle M_{H}$ $\displaystyle=$
$\displaystyle\sqrt{\int_{(4m_{Q})^{2}}^{s_{0}}\rho
se^{-s/M^{2}}ds/\int_{(4m_{Q})^{2}}^{s_{0}}\rho e^{-s/M^{2}}ds}.$ (5)

In the OPE calculation, one could work at the momentum-space making use of the
heavy-quark propagator reinders , and then the result is dimensionally
regularized at $D=4$, by making an extension of the calculation techniques
Nielsen ; Nielsen1 ; Zhang  to fully heavy tetraquark systems. The spectral
density is concretely expressed as $\rho=\rho^{\mbox{pert}}+\rho^{\langle
g^{2}G^{2}\rangle}+\rho^{\langle g^{3}G^{3}\rangle}$, with

$\displaystyle\rho^{\mbox{pert}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2^{10}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{1}{\textbf{H}^{6}}\bigg{[}m_{Q}^{2}-\frac{(1-\alpha-\beta-\gamma)s}{\textbf{H}}\bigg{]}^{2}$
$\displaystyle\times$
$\displaystyle\bigg{\\{}\Big{[}-3(1-\alpha-\beta-\gamma)+4\Big{(}\gamma(1-\alpha-\beta-\gamma)+\alpha\beta\Big{)}\textbf{H}-6\alpha\beta\gamma\textbf{H}^{2}\Big{]}\textbf{H}^{2}m_{Q}^{4}$
$\displaystyle+$
$\displaystyle(1-\alpha-\beta-\gamma)\Big{[}18(1-\alpha-\beta-\gamma)-10\Big{(}\gamma(1-\alpha-\beta-\gamma)+\alpha\beta\Big{)}\textbf{H}\Big{]}\textbf{H}m_{Q}^{2}s-21(1-\alpha-\beta-\gamma)^{3}s^{2}\bigg{\\}},$

$\displaystyle\rho^{\langle g^{2}G^{2}\rangle}$ $\displaystyle=$
$\displaystyle-\frac{m_{Q}^{2}\langle g^{2}G^{2}\rangle}{3\cdot
2^{9}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{(1-\alpha-\beta-\gamma)^{2}}{\textbf{H}^{5}}$
$\displaystyle\times$
$\displaystyle\bigg{\\{}\Big{[}-6(1-\alpha-\beta-\gamma)^{2}+(1-\alpha-\beta-\gamma)\Big{(}6\gamma+2\gamma(1-\alpha-\beta-\gamma)+2\alpha\beta\Big{)}\textbf{H}-3\alpha\beta\gamma\textbf{H}^{2}\Big{]}\textbf{H}m_{Q}^{2}$
$\displaystyle+$ $\displaystyle
3(1-\alpha-\beta-\gamma)^{2}\Big{[}4(1-\alpha-\beta-\gamma)-3\gamma\textbf{H}\Big{]}s\bigg{\\}},$

$\displaystyle\rho^{\langle g^{3}G^{3}\rangle}$ $\displaystyle=$
$\displaystyle-\frac{\langle g^{3}G^{3}\rangle}{3\cdot
2^{10}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{(1-\alpha-\beta-\gamma)^{3}}{\textbf{H}^{5}}\bigg{\\{}\Big{[}-3(1-\alpha-\beta-\gamma)$
$\displaystyle-$ $\displaystyle
6(1-\alpha-\beta-\gamma)^{2}+6\gamma(1-\alpha-\beta-\gamma)\textbf{H}+\alpha\beta\textbf{H}\Big{]}\textbf{H}m_{Q}^{2}+6(1-\alpha-\beta-\gamma)^{2}s\bigg{\\}},$

for the scalar-scalar current,

$\displaystyle\rho^{\mbox{pert}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2^{10}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{1}{\textbf{H}^{6}}\bigg{[}m_{Q}^{2}-\frac{(1-\alpha-\beta-\gamma)s}{\textbf{H}}\bigg{]}^{2}$
$\displaystyle\times$
$\displaystyle\bigg{\\{}\Big{[}-3(1-\alpha-\beta-\gamma)-4\Big{(}\gamma(1-\alpha-\beta-\gamma)+\alpha\beta\Big{)}\textbf{H}-6\alpha\beta\gamma\textbf{H}^{2}\Big{]}\textbf{H}^{2}m_{Q}^{4}$
$\displaystyle+$
$\displaystyle(1-\alpha-\beta-\gamma)\Big{[}18(1-\alpha-\beta-\gamma)+10\Big{(}\gamma(1-\alpha-\beta-\gamma)+\alpha\beta\Big{)}\textbf{H}\Big{]}\textbf{H}m_{Q}^{2}s-21(1-\alpha-\beta-\gamma)^{3}s^{2}\bigg{\\}},$

$\displaystyle\rho^{\langle g^{2}G^{2}\rangle}$ $\displaystyle=$
$\displaystyle-\frac{m_{Q}^{2}\langle g^{2}G^{2}\rangle}{3\cdot
2^{9}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{(1-\alpha-\beta-\gamma)^{2}}{\textbf{H}^{5}}$
$\displaystyle\times$
$\displaystyle\bigg{\\{}\Big{[}-6(1-\alpha-\beta-\gamma)^{2}-(1-\alpha-\beta-\gamma)\Big{(}6\gamma+2\gamma(1-\alpha-\beta-\gamma)+2\alpha\beta\Big{)}\textbf{H}-3\alpha\beta\gamma\textbf{H}^{2}\Big{]}\textbf{H}m_{Q}^{2}$
$\displaystyle+$ $\displaystyle
3(1-\alpha-\beta-\gamma)^{2}\Big{[}4(1-\alpha-\beta-\gamma)+3\gamma\textbf{H}\Big{]}s\bigg{\\}},$

$\displaystyle\rho^{\langle g^{3}G^{3}\rangle}$ $\displaystyle=$
$\displaystyle-\frac{\langle g^{3}G^{3}\rangle}{3\cdot
2^{10}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{(1-\alpha-\beta-\gamma)^{3}}{\textbf{H}^{5}}\bigg{\\{}\Big{[}-3(1-\alpha-\beta-\gamma)$
$\displaystyle-$ $\displaystyle
6(1-\alpha-\beta-\gamma)^{2}-6\gamma(1-\alpha-\beta-\gamma)\textbf{H}-\alpha\beta\textbf{H}\Big{]}\textbf{H}m_{Q}^{2}+6(1-\alpha-\beta-\gamma)^{2}s\bigg{\\}},$

for the pseudoscalar-pseudoscalar current,

$\displaystyle\rho^{\mbox{pert}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2^{8}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{1}{\textbf{H}^{6}}\bigg{[}m_{Q}^{2}-\frac{(1-\alpha-\beta-\gamma)s}{\textbf{H}}\bigg{]}^{2}$
$\displaystyle\times$
$\displaystyle\bigg{\\{}\Big{[}-3(1-\alpha-\beta-\gamma)+2\Big{(}\gamma(1-\alpha-\beta-\gamma)+\alpha\beta\Big{)}\textbf{H}-6\alpha\beta\gamma\textbf{H}^{2}\Big{]}\textbf{H}^{2}m_{Q}^{4}$
$\displaystyle+$
$\displaystyle(1-\alpha-\beta-\gamma)\Big{[}18(1-\alpha-\beta-\gamma)-5\Big{(}\gamma(1-\alpha-\beta-\gamma)+\alpha\beta\Big{)}\textbf{H}\Big{]}\textbf{H}m_{Q}^{2}s-21(1-\alpha-\beta-\gamma)^{3}s^{2}\bigg{\\}},$

$\displaystyle\rho^{\langle g^{2}G^{2}\rangle}$ $\displaystyle=$
$\displaystyle-\frac{m_{Q}^{2}\langle g^{2}G^{2}\rangle}{3\cdot
2^{8}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{(1-\alpha-\beta-\gamma)^{2}}{\textbf{H}^{5}}$
$\displaystyle\times$
$\displaystyle\bigg{\\{}\Big{[}-12(1-\alpha-\beta-\gamma)^{2}+(1-\alpha-\beta-\gamma)\Big{(}6\gamma+2\gamma(1-\alpha-\beta-\gamma)+2\alpha\beta\Big{)}\textbf{H}-6\alpha\beta\gamma\textbf{H}^{2}\Big{]}\textbf{H}m_{Q}^{2}$
$\displaystyle+$ $\displaystyle
3(1-\alpha-\beta-\gamma)^{2}\Big{[}8(1-\alpha-\beta-\gamma)-3\gamma\textbf{H}\Big{]}s\bigg{\\}},$

$\displaystyle\rho^{\langle g^{3}G^{3}\rangle}$ $\displaystyle=$
$\displaystyle-\frac{\langle g^{3}G^{3}\rangle}{3\cdot
2^{9}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{(1-\alpha-\beta-\gamma)^{3}}{\textbf{H}^{5}}\bigg{\\{}\Big{[}-6(1-\alpha-\beta-\gamma)$
$\displaystyle-$ $\displaystyle
12(1-\alpha-\beta-\gamma)^{2}+6\gamma(1-\alpha-\beta-\gamma)\textbf{H}+\alpha\beta\textbf{H}\Big{]}\textbf{H}m_{Q}^{2}+12(1-\alpha-\beta-\gamma)^{2}s\bigg{\\}},$

for the axial-axial current, and

$\displaystyle\rho^{\mbox{pert}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2^{8}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{1}{\textbf{H}^{6}}\bigg{[}m_{Q}^{2}-\frac{(1-\alpha-\beta-\gamma)s}{\textbf{H}}\bigg{]}^{2}$
$\displaystyle\times$
$\displaystyle\bigg{\\{}\Big{[}-3(1-\alpha-\beta-\gamma)-2\Big{(}\gamma(1-\alpha-\beta-\gamma)+\alpha\beta\Big{)}\textbf{H}-6\alpha\beta\gamma\textbf{H}^{2}\Big{]}\textbf{H}^{2}m_{Q}^{4}$
$\displaystyle+$
$\displaystyle(1-\alpha-\beta-\gamma)\Big{[}18(1-\alpha-\beta-\gamma)+5\Big{(}\gamma(1-\alpha-\beta-\gamma)+\alpha\beta\Big{)}\textbf{H}\Big{]}\textbf{H}m_{Q}^{2}s-21(1-\alpha-\beta-\gamma)^{3}s^{2}\bigg{\\}},$

$\displaystyle\rho^{\langle g^{2}G^{2}\rangle}$ $\displaystyle=$
$\displaystyle-\frac{m_{Q}^{2}\langle g^{2}G^{2}\rangle}{3\cdot
2^{8}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{(1-\alpha-\beta-\gamma)^{2}}{\textbf{H}^{5}}$
$\displaystyle\times$
$\displaystyle\bigg{\\{}\Big{[}-12(1-\alpha-\beta-\gamma)^{2}-(1-\alpha-\beta-\gamma)\Big{(}6\gamma+2\gamma(1-\alpha-\beta-\gamma)+2\alpha\beta\Big{)}\textbf{H}-6\alpha\beta\gamma\textbf{H}^{2}\Big{]}\textbf{H}m_{Q}^{2}$
$\displaystyle+$ $\displaystyle
3(1-\alpha-\beta-\gamma)^{2}\Big{[}8(1-\alpha-\beta-\gamma)+3\gamma\textbf{H}\Big{]}s\bigg{\\}},$

$\displaystyle\rho^{\langle g^{3}G^{3}\rangle}$ $\displaystyle=$
$\displaystyle-\frac{\langle g^{3}G^{3}\rangle}{3\cdot
2^{9}\pi^{6}}\int_{\alpha_{min}}^{\alpha_{max}}\frac{d\alpha}{\alpha^{3}}\int_{\beta_{min}}^{\beta_{max}}\frac{d\beta}{\beta^{3}}\int_{\gamma_{min}}^{\gamma_{max}}\frac{d\gamma}{\gamma^{3}}\frac{(1-\alpha-\beta-\gamma)^{3}}{\textbf{H}^{5}}\bigg{\\{}\Big{[}-6(1-\alpha-\beta-\gamma)$
$\displaystyle-$ $\displaystyle
12(1-\alpha-\beta-\gamma)^{2}-6\gamma(1-\alpha-\beta-\gamma)\textbf{H}-\alpha\beta\textbf{H}\Big{]}\textbf{H}m_{Q}^{2}+12(1-\alpha-\beta-\gamma)^{2}s\bigg{\\}},$

for the vector-vector current. It is defined as

H $\displaystyle=$ $\displaystyle
1+(1-\alpha-\beta-\gamma)\Big{(}\frac{1}{\alpha}+\frac{1}{\beta}+\frac{1}{\gamma}\Big{)},$
$\displaystyle\alpha_{min}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigg{[}\bigg{(}1-\frac{8m_{Q}^{2}}{s}\bigg{)}-\sqrt{\bigg{(}1-\frac{8m_{Q}^{2}}{s}\bigg{)}^{2}-\frac{4m_{Q}^{2}}{s}}\Bigg{]},$
$\displaystyle\alpha_{max}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigg{[}\bigg{(}1-\frac{8m_{Q}^{2}}{s}\bigg{)}+\sqrt{\bigg{(}1-\frac{8m_{Q}^{2}}{s}\bigg{)}^{2}-\frac{4m_{Q}^{2}}{s}}\Bigg{]},$
$\displaystyle\beta_{min}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigg{[}\Big{(}1+2\alpha-\frac{3\alpha^{2}s}{\alpha
s-m_{Q}^{2}}\Big{)}-\sqrt{\frac{[\alpha(1-\alpha)s-m_{Q}^{2}][\alpha(1-\alpha)s-(1+8\alpha)m_{Q}^{2}]}{(\alpha
s-m_{Q}^{2})^{2}}}\Bigg{]},$ $\displaystyle\beta_{max}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigg{[}\Big{(}1+2\alpha-\frac{3\alpha^{2}s}{\alpha
s-m_{Q}^{2}}\Big{)}+\sqrt{\frac{[\alpha(1-\alpha)s-m_{Q}^{2}][\alpha(1-\alpha)s-(1+8\alpha)m_{Q}^{2}]}{(\alpha
s-m_{Q}^{2})^{2}}}\Bigg{]},$ $\displaystyle\gamma_{min}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigg{[}(1-\alpha-\beta)-\sqrt{\frac{(1-\alpha-\beta)\Big{\\{}(1-\alpha-\beta)[\alpha\beta
s-(\alpha+\beta)m_{Q}^{2}]-4\alpha\beta m_{Q}^{2}\Big{\\}}}{\alpha\beta
s-(\alpha+\beta)m_{Q}^{2}}}\Bigg{]},$

and

$\displaystyle\gamma_{max}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigg{[}(1-\alpha-\beta)+\sqrt{\frac{(1-\alpha-\beta)\Big{\\{}(1-\alpha-\beta)[\alpha\beta
s-(\alpha+\beta)m_{Q}^{2}]-4\alpha\beta m_{Q}^{2}\Big{\\}}}{\alpha\beta
s-(\alpha+\beta)m_{Q}^{2}}}\Bigg{]}.$

##  III Numerical analysis and discussions

In this part, the heavy $m_{Q}$ is taken as the running charm mass
$m_{c}=1.27\pm 0.02~{}\mbox{GeV}$ PDG , and other input parameters are
$\langle g^{2}G^{2}\rangle=0.88\pm 0.25~{}\mbox{GeV}^{4}$ and $\langle
g^{3}G^{3}\rangle=0.58\pm 0.18~{}\mbox{GeV}^{6}$ svzsum ; overview3 ; Narison
. Following the standard criterion of sum rule analysis, one could find proper
work windows for the threshold parameter $\sqrt{s_{0}}$ and the Borel
parameter $M^{2}$. The lower bound of $M^{2}$ could be gained in view of the
OPE convergence, and the upper one is obtained from the pole dominance.
Meanwhile, the threshold $\sqrt{s_{0}}$ is around $0.4\sim 0.6~{}\mbox{GeV}$
higher than the extracted $M_{H}$ empirically, for which describes the
beginning of continuum state.

Taking the scalar-scalar case an example, the inputs are kept at their central
values at the start. To find the lower bound of $M^{2}$, its OPE convergence
is shown by comparing the relative contributions of various condensates from
sum rule (4) for $\sqrt{s_{0}}=7.0~{}\mbox{GeV}$ in FIG. 1, and one could note
that the relative contributions of two-gluon condensate $\langle
g^{2}G^{2}\rangle$ and three-gluon condensate $\langle g^{3}G^{3}\rangle$ are
very small. Numerically, it is taken as $M^{2}\geq 2.5~{}\mbox{GeV}^{2}$ with
an eye to the OPE convergence analysis. On the other hand, the upper one of
$M^{2}$ is obtained from the pole contribution dominance phenomenologically.
FIG. 2 makes the comparison between pole and continuum contribution from sum
rule (4) for $\sqrt{s_{0}}=7.0~{}\mbox{GeV}$. The relative pole contribution
is about $50\%$ at $M^{2}=3.0~{}\mbox{GeV}^{2}$ and decreasing with $M^{2}$.
Thereby, it could satisfy the pole dominance requirement while $M^{2}\leq
3.0~{}\mbox{GeV}^{2}$, and the Borel window is fixed as $M^{2}=2.5\sim
3.0~{}\mbox{GeV}^{2}$ for $\sqrt{s_{0}}=7.0~{}\mbox{GeV}$. Analogously, they
are taken as $M^{2}=2.5\sim 2.9~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=6.9~{}\mbox{GeV}$, and $M^{2}=2.5\sim 3.2~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.1~{}\mbox{GeV}$, respectively. The mass $M_{H}$ dependence on
$M^{2}$ is shown in FIG. 3 for the scalar-scalar case, and it is computed to
be $6.44\pm 0.13~{}\mbox{GeV}$ in the chosen windows. And then varying all the
input values, the attained mass is $6.44\pm 0.13^{+0.02}_{-0.03}~{}\mbox{GeV}$
(the first error from variation of $s_{0}$ and $M^{2}$, and the second from
the uncertainty of QCD parameters) or $6.44^{+0.15}_{-0.16}~{}\mbox{GeV}$ in a
short form.

In the very similar analyzing processes, proper work windows for other three
cases could also be found and their corresponding Borel curves are
respectively given in FIG. 4–6. After considering the uncertainty from both
work windows and variation of input parameters, mass values of $0^{+}$
$cc\bar{c}\bar{c}$ tetraquark states are gained as
$6.45^{+0.14}_{-0.16}~{}\mbox{GeV}$ for the pseudoscalar-pseudoscalar
configuration, $6.46^{+0.13}_{-0.17}~{}\mbox{GeV}$ for the axial-axial one,
and $6.47^{+0.12}_{-0.18}~{}\mbox{GeV}$ for the vector-vector one,
respectively.

Figure 1: The OPE convergence for the $0^{+}$ fully-charmed tetraquark state
with a scalar-scalar configuration is shown by comparing the relative
contributions of perturbative, two-gluon condensate $\langle
g^{2}G^{2}\rangle$, and three-gluon condensate $\langle g^{3}G^{3}\rangle$
from sum rule (4) for $\sqrt{s_{0}}=7.0~{}\mbox{GeV}$.

Figure 2: The phenomenological contribution in sum rule (4) for
$\sqrt{s_{0}}=7.0~{}\mbox{GeV}$ for the $0^{+}$ fully-charmed tetraquark state
with a scalar-scalar configuration. The solid line is the relative pole
contribution as a function of $M^{2}$ and the dashed line is the relative
continuum contribution.

Figure 3: The mass $M_{H}$ dependence on $M^{2}$ for the $0^{+}$ fully-charmed
tetraquark state with a scalar-scalar configuration from sum rule (5) is
shown. The Borel windows of $M^{2}$ are $2.5\sim 2.9~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=6.9~{}\mbox{GeV}$, $2.5\sim 3.0~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.0~{}\mbox{GeV}$, and $2.5\sim 3.2~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.1~{}\mbox{GeV}$, respectively.

Figure 4: The mass $M_{H}$ dependence on $M^{2}$ for the $0^{+}$ fully-charmed
tetraquark state with a pseudoscalar-pseudoscalar configuration from sum rule
(5) is shown. The Borel windows of $M^{2}$ are $2.5\sim 2.9~{}\mbox{GeV}^{2}$
for $\sqrt{s_{0}}=6.9~{}\mbox{GeV}$, $2.5\sim 3.0~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.0~{}\mbox{GeV}$, and $2.5\sim 3.2~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.1~{}\mbox{GeV}$, respectively.

Figure 5: The mass $M_{H}$ dependence on $M^{2}$ for the $0^{+}$ fully-charmed
tetraquark state with a axial-axial configuration from sum rule (5) is shown.
The Borel windows of $M^{2}$ are $2.5\sim 2.9~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=6.9~{}\mbox{GeV}$, $2.5\sim 3.0~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.0~{}\mbox{GeV}$, and $2.5\sim 3.2~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.1~{}\mbox{GeV}$, respectively.

Figure 6: The mass $M_{H}$ dependence on $M^{2}$ for the $0^{+}$ fully-charmed
tetraquark state with a vector-vector configuration from sum rule (5) is
shown. The Borel windows of $M^{2}$ are $2.5\sim 2.9~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=6.9~{}\mbox{GeV}$, $2.5\sim 3.0~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.0~{}\mbox{GeV}$, and $2.5\sim 3.2~{}\mbox{GeV}^{2}$ for
$\sqrt{s_{0}}=7.1~{}\mbox{GeV}$, respectively.

##  IV Summary

Focusing on the LHCb’s new observation in the di-$J/\psi$ mass spectrum, we
systematically investigate $0^{+}$ fully-charmed tetraquark states in the
framework of QCD sum rules. By developing related calculation techniques to
fourfold heavy tetraquark states, four types of currents with different
configurations are taken into consideration and condensates up to dimension
$6$ are involved in the OPE side. At last, mass values of $0^{+}$
$cc\bar{c}\bar{c}$ tetraquark states are calculated to be
$6.44^{+0.15}_{-0.16}~{}\mbox{GeV}$ for the scalar-scalar case,
$6.45^{+0.14}_{-0.16}~{}\mbox{GeV}$ for the pseudoscalar-pseudoscalar case,
$6.46^{+0.13}_{-0.17}~{}\mbox{GeV}$ for the axial-axial case, and
$6.47^{+0.12}_{-0.18}~{}\mbox{GeV}$ for the vector-vector case, respectively.
All these results are numerically consistent with the LHCb’s experimental data
$6.2\sim 6.8~{}\mbox{GeV}$ for the broad structure, which could support its
internal structure as a $0^{+}$ $cc\bar{c}\bar{c}$ tetraquark state. For the
future, it is expected that further experimental and theoretical efforts may
reveal more on the nature of the exotic states.

###### Acknowledgements.

The author is very grateful to Xiang Liu for first bringing the information
“Latest results on exotic hadrons at LHCb” to his attention, and also thanks
Zhi-Gang Wang for recent communication and discussion. This work was supported
by the National Natural Science Foundation of China under Contract Nos.
11475258 and 11675263, and by the project for excellent youth talents in NUDT.

## References

  * (1) Y. Iwasaki, Phys. Rev. Lett. 36, 1266 (1976).
  * (2) K. T. Chao, Z. Phys. C 7, 317 (1981).
  * (3) J. P. Ader, J. M. Richard, and P. Taxil, Phys. Rev. D 25, 2370 (1982).
  * (4) A. M. Badalian, B. L. Ioffe, and A. V. Smilga, Nucl. Phys. B281, 85 (1987).
  * (5) R. J. Lloyd and J. P. Vary, Phys. Rev. D 70, 014009 (2004).
  * (6) N. Barnea, J. Vijande, and A. Valcarce, Phys. Rev. D 73, 054004 (2006).
  * (7) A. V. Berezhnoy, A. V. Luchinsky, and A. A. Novoselov, Phys. Rev. D 86, 034004 (2012).
  * (8) A. V. Berezhnoy, A. K. Likhoded, and A. A. Novoselov, Phys. Rev. D 87, 054023 (2013).
  * (9) M. Karliner, S. Nussinov, and J. L. Rosner, Phys. Rev. D 95, 034011 (2017).
  * (10) J. M. Richard, A. Valcarce, and J. Vijande, Phys. Rev. D 95, 054019 (2017).
  * (11) W. Chen, H. X. Chen, X. Liu, T. G. Steele, and S. L. Zhu, Phys. Lett. B 773, 247 (2017).
  * (12) Z. G. Wang, Eur. Phys. J. C 77, 432 (2017).
  * (13) Y. Chen and R. Vega-Morales, arXiv:1710.02738 [hep-ph].
  * (14) J. Wu, Y. R. Liu, K. Chen, X. Liu, and S. L. Zhu, Phys. Rev. D 97, 094015 (2018).
  * (15) M. N. Anwar, J. Ferretti, F. K. Guo, E. Santopinto, and B. S. Zou, Eur. Phys. J. C 78, 647 (2018).
  * (16) A. Esposito and A. D. Polosa, Eur. Phys. J. C 78, 782 (2018).
  * (17) C. Hughes, E. Eichten, and C. T. H. Davies, Phys. Rev. D 97, 054505 (2018).
  * (18) V. R. Debastiani and F. S. Navarra, Chin. Phys. C 43, 013105 (2019).
  * (19) M. S. Liu, Q. F. Lü, X. H. Zhong, and Q. Zhao, Phys. Rev. D 100, 016006 (2019).
  * (20) G. J. Wang, L. Meng, and S. L. Zhu, Phys. Rev. D 100, 096013 (2019).
  * (21) M. A. Bedolla, J. Ferretti, C. D. Roberts, and E. Santopinto, arXiv:1911.00960 [hep-ph].
  * (22) Z. G. Wang and Z. Y. Di, Acta Phys. Polon. B 50, 1335 (2019).
  * (23) X. Chen, Phys. Rev. D 100, 094009 (2019).
  * (24) Y. R. Liu, H. X. Chen, W. Chen, X. Liu, and S. L. Zhu, Prog. Part. Nucl. Phys. 107, 237 (2019).
  * (25) Y. Bai, S. Lu, and J. Osborne, Phys. Lett. B 798, 134930 (2019).
  * (26) C. Becchi, A. Giachino, L. Maiani, and E. Santopinto, Phys. Lett. B 806, 135495 (2020).
  * (27) X. Chen, arXiv:2001.06755 [hep-ph].
  * (28) C. Deng, H. Chen, and J. Ping, arXiv:2003.05154 [hep-ph].
  * (29) R. Aaij et al., (LHCb collaboration), arXiv:2006.16957 [hep-ex].
  * (30) M. S. Liu, F. X. Liu, X. H. Zhong, and Q. Zhao, arXiv:2006.11952 [hep-ph].
  * (31) Z. G. Wang, arXiv:2006.13028 [hep-ph].
  * (32) X. Jin, Y. Xue, H. Huang, and J. Ping, arXiv:2006.13745 [hep-ph].
  * (33) G. Yang, J. Ping, L. He, and Q. Wang, arXiv:2006.13756 [hep-ph].
  * (34) C. Becchi, J. Ferretti, A. Giachino, L. Maiani, and E. Santopinto, arXiv:2006.14388 [hep-ph].
  * (35) Q. F. Lü, D. Y. Chen, Y. B. Dong, Eur. Phys. J. C 80, 871 (2020).
  * (36) H. X. Chen, W. Chen, X. Liu, and S. L. Zhu, arXiv:2006.16027 [hep-ph].
  * (37) M. Mikhasenko, L. An, and R. McNulty, arXiv:2007.05501 [hep-ph].
  * (38) X. Y. Wang, Q. Y. Lin, H. Xu, Y. P. Xie, Y. Huang, and X. Chen, arXiv:2007.09697 [hep-ph].
  * (39) M. Albaladejo et al., (Joint Physics Analysis Center), arXiv:2008.01001 [hep-ph].
  * (40) J. Sonnenschein and D. Weissman, arXiv:2008.01095 [hep-ph].
  * (41) R. M. Albuquerque, S. Narison, A. Rabemananjara, D. Rabetiarivony, and G. Randriamanatrika, arXiv:2008.01569 [hep-ph].
  * (42) J. F. Giron and R. F. Lebed, arXiv:2008.01631 [hep-ph].
  * (43) L. Maiani, arXiv:2008.01637 [hep-ph].
  * (44) J. M. Richard, arXiv:2008.01962 [hep-ph].
  * (45) J. Z. Wang, D. Y. Chen, X. Liu, and T. Matsuki, arXiv:2008.07430 [hep-ph].
  * (46) K. T. Chao and S. L. Zhu, arXiv:2008.07670 [hep-ph].
  * (47) G. Yang, J. Ping, and J. Segovia, arXiv:2009.00238 [hep-ph].
  * (48) R. Maciuła, W. Schäfer, and A. Szczurek, arXiv:2009.02100 [hep-ph].
  * (49) M. Karliner and J. L. Rosner, arXiv:2009.04429 [hep-ph].
  * (50) Z. G. Wang, arXiv:2009.05371 [hep-ph].
  * (51) X. K. Dong, V. Baru, F. K. Guo, C. Hanhart, and A. Nefediev, arXiv:2009.07795 [hep-ph].
  * (52) Y. Q. Ma and H. F. Zhang, arXiv:2009.08376 [hep-ph].
  * (53) F. Feng, Y. Huang, Y. Jia, W. L. Sang, X. Xiong, and J. Y. Zhang, arXiv:2009.08450 [hep-ph].
  * (54) J. Zhao, S. Shi, and P. Zhuang, arXiv:2009.10319 [hep-ph].
  * (55) M. C. Gordillo, F. De Soto, and J. Segovia, arXiv:2009.11889 [hep-ph].
  * (56) R. N. Faustov, V. O. Galkin, and E. M. Savchenko, arXiv:2009.13237 [hep-ph].
  * (57) X. Z. Weng, X. L. Chen, W. Z. Deng, and S. L. Zhu, arXiv:2010.05163 [hep-ph].
  * (58) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B147, 385 (1979); B147, 448 (1979); V. A. Novikov, M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Fortschr. Phys. 32, 585 (1984).
  * (59) B. L. Ioffe, in _The Spin Structure of The Nucleon_ , edited by B. Frois, V. W. Hughes, and N. de Groot (World Scientific, Singapore, 1997).
  * (60) S. Narison, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 17, 1 (2002), arXiv:hep-ph/0205006.
  * (61) P. Colangelo and A. Khodjamirian, in _At the Frontier of Particle Physics: Handbook of QCD_ , edited by M. Shifman, Boris Ioffe Festschrift Vol. 3 (World Scientific, Singapore, 2001), pp. 1495-1576.
  * (62) M. Nielsen, F. S. Navarra, and S. H. Lee, Phys. Rep. 497, 41 (2010).
  * (63) H. Kim and Y. Oh, Phys. Rev. D 72, 074012 (2005); M. E. Bracco, A. Lozea, R. D. Matheus, F. S. Navarra, and M. Nielsen, Phys. Lett. B 624, 217 (2005); R. D. Matheus, S. Narison, M. Nielsen, and J. M. Richard, Phys. Rev. D 75, 014005 (2007).
  * (64) F. S. Navarra, M. Nielsen, and S. H. Lee, Phys. Lett. B 649, 166 (2007); S. H. Lee, A. Mihara, F. S. Navarra, and M. Nielsen, Phys. Lett. B 661, 28 (2008).
  * (65) J. R. Zhang and M. Q. Huang, JHEP 1011, 057 (2010); J. R. Zhang and M. Q. Huang, Phys. Rev. D 83, 036005 (2011); J. R. Zhang, J. L. Zou, and J. Y. Wu, Chin. Phys. C 42, 043101 (2018); J. R. Zhang, Phys. Lett. B 789, 432 (2019); J. R. Zhang, Phys. Rev. D 102, 054006 (2020); J. R. Zhang, Phys. Rev. D 102, 074013 (2020).
  * (66) A. Martínez Torres, K. P. Khemchandani, M. Nielsen, F. S. Navarra, and E. Oset, Phys. Rev. D 88, 074033 (2013).
  * (67) K. P. Khemchandani, A. Martínez Torres, M. Nielsen, and F. S. Navarra, Phys. Rev. D 89, 014029 (2014).
  * (68) L. J. Reinders, H. R. Rubinstein, and S. Yazaki, Phys. Rep. 127, 1 (1985).
  * (69) M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018) and 2019 update.
  * (70) S. Narison, Phys. Rep. 84, 263 (1982); G. Launer, S. Narison, and R. Tarrach, Z. Phys. C 26, 433 (1984); S. Narison, Phys. Lett. B 673, 30 (2009).


# Multi-Task Graph Autoencoders

Phi Vu Tran
Strategic Innovation Group
Booz $|$ Allen $|$ Hamilton
San Diego, CA USA
<EMAIL_ADDRESS>

##  1 Autoencoder Architecture for Link Prediction and Node Classification

As the world is becoming increasingly interconnected, relational data are also
growing in ubiquity. In this work, we examine the task of learning to make
predictions on graphs for a broad range of real-world applications.
Specifically, we study two canonical subtasks associated with relational,
graph-structured datasets: link prediction and node classification (LPNC). A
graph is a partially observed set of edges and nodes (or vertices), and the
learning task is to predict the labels for edges and nodes. In real-world
applications, the input graph is a network with nodes representing unique
entities and edges representing relationships (or links) between entities.
Further, the labels of nodes and edges in a graph are often correlated,
exhibiting complex relational structures that violate the general assumption
of independent and identical distribution fundamental in traditional machine
learning [5]. Therefore, models capable of exploiting topological structures
of graphs have been shown to achieve superior predictive performances on many
LPNC tasks [10].

We introduce the Multi-Task Graph Autoencoder (MTGAE) architecture,
schematically depicted in Figure 1, capable of learning a shared
representation of latent node embeddings from local graph topology and
available explicit node features for LPNC. Our simple, yet effective and
versatile model is efficiently trained end-to-end for the joint, simultaneous
multi-task learning (MTL) of unsupervised link prediction and semi-supervised
node classification in a single stage, whereas previous related deep graph
embedding methods require multiple training steps that are difficult to
optimize. We present an empirical evaluation of the MTGAE model on five
challenging benchmark graph-structured datasets and show significant
improvement in predictive performance over three strong baselines designed
specifically for LPNC. Reference code and data are available at
https://github.com/vuptran/graph-representation-learning.

Figure 1: Schematic depiction of the Multi-Task Graph Autoencoder (MTGAE)
architecture. _Left_ : A partially observed graph with positive links (solid
lines) and negative links (dashed lines) between two nodes; pairs of nodes not
yet connected have unknown status links. _Middle_ : A symmetrical, densely
connected autoencoder with parameter sharing is trained end-to-end to learn
node embeddings from the adjacency vector. _Right_ : Exemplar multi-task
output for link prediction and node classification.

Problem Formulation and Notation The input to the MTGAE model is a graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ of $N=|\mathcal{V}|$ nodes. Graph
$\mathcal{G}$ is represented by its adjacency matrix
$\mathbf{A}\in\mathbb{R}^{N\times N}$ paired with a unique ordering of
vertices. For a partially observed graph,
$\mathbf{A}\in\\{1,0,\textsc{unk}\\}^{N\times N}$, where $1$ denotes a known
positive edge, $0$ denotes a known negative edge, and unk denotes an unknown
status (missing or unobserved) edge. In general, the input to the model can be
directed or undirected, weighted or unweighted, and/or bipartite graphs.

Optionally, we are given a matrix of available node features, i.e. side
information $\mathbf{X}\in\mathbb{R}^{N\times F}$. The aim of the MTGAE model
$h(\mathbf{A,X})$ is to learn a set of low-dimensional latent variables for
the nodes $\mathbf{Z}\in\mathbb{R}^{N\times D}$ that can produce an
approximate reconstruction output $\mathbf{\hat{A}}$ such that the empirical
error between $\mathbf{A}$ and $\mathbf{\hat{A}}$ is minimized, thereby
preserving the global graph structure. In this paper, we use capital variables
(e.g., $\mathbf{A}$) to denote matrices and lower-case variables (e.g.,
$\mathbf{a}$) to denote row vectors. For example, we use $\mathbf{a}_{i}$ to
mean the $i$th row of the matrix $\mathbf{A}$.

Unsupervised Link Prediction Let $\mathbf{a}_{i}\in\mathbb{R}^{N}$ be an
_adjacency vector_ of $\mathbf{A}$ that contains the local neighborhood of the
$i$th node. Our proposed MTGAE architecture comprises a set of non-linear
transformations on $\mathbf{a}_{i}$ summarized in two component parts: encoder
$g(\mathbf{a}_{i})\colon\mathbb{R}^{N}\to\mathbb{R}^{D}$ and decoder
$f\left(g\left(\mathbf{a}_{i}\right)\right)\colon\mathbb{R}^{D}\to\mathbb{R}^{N}$.
We stack two layers of the encoder part to derive $D$-dimensional latent
feature representation of the $i$th node $\mathbf{z}_{i}\in\mathbb{R}^{D}$,
and then stack two layers of the decoder part to obtain an approximate
reconstruction output $\mathbf{\hat{a}}_{i}\in\mathbb{R}^{N}$, resulting in a
four-layer autoencoder architecture. Note that $\mathbf{a}_{i}$ is highly
sparse, with up to 80 percent of the edges missing at random in some of our
experiments, and the dense reconstructed output $\mathbf{\hat{a}}_{i}$
contains the predictions for the missing edges. The hidden representations for
the encoder and decoder parts are computed as follows:

$\displaystyle\text{Encoder}\qquad\mathbf{z}_{i}$
$\displaystyle=g\left(\mathbf{a}_{i}\right)=\text{ReLU}\left(\mathbf{W}\boldsymbol{\cdot}\text{ReLU}\left(\mathbf{V}\mathbf{a}_{i}+\mathbf{b}^{(1)}\right)+\mathbf{b}^{(2)}\right).$
$\displaystyle\text{Decoder}\qquad\mathbf{\hat{a}}_{i}$
$\displaystyle=f\left(\mathbf{z}_{i}\right)=\mathbf{V}^{\text{T}}\boldsymbol{\cdot}\text{ReLU}\left(\mathbf{W}^{\text{T}}\mathbf{z}_{i}+\mathbf{b}^{(3)}\right)+\mathbf{b}^{(4)}.$
$\displaystyle\text{Autoencoder}\qquad\mathbf{\hat{a}}_{i}$
$\displaystyle=h\left(\mathbf{a}_{i}\right)=f\left(g\left(\mathbf{a}_{i}\right)\right).$

The choice of non-linear, element-wise activation function is the rectified
linear unit $\text{ReLU}(\mathbf{x})=\text{max}(0,\mathbf{x})$. The last
decoder layer computes a linear transformation to score the missing links as
part of the reconstruction. We constrain the MTGAE architecture to be
symmetrical with shared parameters for $\\{\mathbf{W},\mathbf{V}\\}$ between
the encoder and decoder parts, resulting in almost $2\times$ fewer parameters
than an unconstrained architecture. Parameter sharing is a powerful form of
regularization that helps improve learning and generalization, and is also the
main motivation behind MTL [2, 16]. Notice the bias units $\mathbf{b}$ do not
share parameters, and
$\left\\{\mathbf{W}^{\text{T}},\mathbf{V}^{\text{T}}\right\\}$ are transposed
copies of $\\{\mathbf{W},\mathbf{V}\\}$. For brevity of notation, we summarize
the parameters to be learned in
$\theta=\left\\{\mathbf{W},\mathbf{V},\mathbf{b}^{(k)}\right\\},k=1,...,4$.

Optionally, if a matrix of node features $\mathbf{X}\in\mathbb{R}^{N\times F}$
is available, then we concatenate $(\mathbf{A},\mathbf{X})$ to obtain an
_augmented_ adjacency matrix $\mathbf{\bar{A}}\in\mathbb{R}^{N\times(N+F)}$
and perform the above encoder-decoder transformations on
$\mathbf{\bar{a}}_{i}$ for unsupervised link prediction. The intuition behind
the concatenation of node features is to enable a shared representation of
both graph and node features throughout the autoencoding transformations by
way of the tied parameters $\\{\mathbf{W},\mathbf{V}\\}$.

During the forward pass, or inference, the model takes as input an adjacency
vector $\mathbf{a}_{i}$ and computes its reconstructed output
$\mathbf{\hat{a}}_{i}=h(\mathbf{a}_{i})$ for unsupervised link prediction.
During the backward pass, we learn $\theta$ by minimizing the Masked Balanced
Cross-Entropy (MBCE) loss, which allows only the contributions of those
parameters associated with observed edges, as in [9, 11]. We handle class
imbalance in link prediction by defining a weighting factor $\zeta\in[0,1]$ to
be used as a multiplier for the positive class in the cross-entropy loss
formulation. For a single example $\mathbf{a}_{i}$ and its reconstructed
output $\mathbf{\hat{a}}_{i}$, we compute the MBCE loss as follows:

$\displaystyle\mathcal{L}_{\textsc{bce}}=-\mathbf{a}_{i}\log\left(\sigma\left(\mathbf{\hat{a}}_{i}\right)\right)\cdot\zeta-(1-\mathbf{a}_{i})$
$\displaystyle\log\left(1-\sigma\left(\mathbf{\hat{a}}_{i}\right)\right),$
$\displaystyle\mathcal{L}_{\textsc{mbce}}=\frac{\sum_{i}\mathbf{m}_{i}\odot\mathcal{L}_{\textsc{bce}}}{\sum_{i}\mathbf{m}_{i}}.$

Here, $\mathcal{L}_{\textsc{bce}}$ is the balanced binary cross-entropy loss
with weighting factor $\zeta=1-\frac{\text{\\# positive links}}{\text{\\#
negative links}}$, $\sigma(\cdot)$ is the sigmoid function, $\odot$ is the
Hadamard (element-wise) product, and $\mathbf{m}_{i}$ is the Boolean mask:
$\mathbf{m}_{i}=1$ if $\mathbf{a}_{i}\neq\textsc{unk}$, else
$\mathbf{m}_{i}=0$.

Semi-Supervised Node Classification The MTGAE model can also be used to
perform efficient information propagation on graphs for the task of semi-
supervised node classification. For a given augmented adjacency vector
$\mathbf{\bar{a}}_{i}$, the model learns the corresponding node embedding
$\mathbf{z}_{i}$ to obtain an optimal reconstruction. Intuitively,
$\mathbf{z}_{i}$ encodes a vector of latent features derived from the
concatenation of both graph and node features, and can be used to predict the
label of the $i$th node. For multi-class classification, we decode
$\mathbf{z}_{i}$ using the softmax activation function to produce a
probability distribution over node labels. More precisely, we predict node
labels via the following transformation:
$\mathbf{\hat{y}}_{i}=\text{softmax}(\mathbf{\tilde{z}}_{i})=\frac{1}{\mathcal{Z}}\exp(\mathbf{\tilde{z}}_{i})$,
where $\mathcal{Z}=\sum_{i}\exp(\mathbf{\tilde{z}}_{i})$ and
$\mathbf{\tilde{z}}_{i}=\mathbf{U}\boldsymbol{\cdot}\text{ReLU}\left(\mathbf{W}^{\text{T}}\mathbf{z}_{i}+\mathbf{b}^{(3)}\right)+\mathbf{b}^{(5)}$.

Multi-Task Learning In many applications, such as knowledge base completion
and network analysis, the input graph is partially observed with an incomplete
set of edges and a small fraction of labeled nodes. Thus, it is desirable for
a model to predict the labels of missing links and nodes _simultaneously_ in a
multi-task learning setting. We achieve multi-task learning on graphs by
training the MTGAE model using a joint loss function that combines the masked
categorical cross-entropy loss for semi-supervised node classification with
the MBCE loss for unsupervised link prediction:

$\mathcal{L}_{\textsc{multi-task}}=\overbrace{-\textsc{mask}_{i}\sum_{c\in
C}\mathbf{y}_{ic}\log(\mathbf{\hat{y}}_{ic})}^{\text{semi-supervised
classification}}+\mathcal{L}_{\textsc{mbce}},$

where $C$ is the set of node labels, $\mathbf{y}_{ic}$ is the binary indicator
if node $i$ belongs to class $c$, $\mathbf{\hat{y}}_{ic}$ is the softmax
probability that node $i$ belongs to class $c$, $\mathcal{L}_{\textsc{mbce}}$
is the loss defined for unsupervised link prediction, and $\textsc{mask}_{i}$
is the Boolean variable: $\textsc{mask}_{i}=1$ if node $i$ has a label, else
$\textsc{mask}_{i}=0$.

The training complexity of the MTGAE model is $\mathcal{O}((N+F)DI)$, where
$N$ is the number of nodes, $F$ is the dimensionality of node features, $D$ is
the size of the hidden layer, and $I$ is the number of iterations. In
practice, $F$, $D\ll N$, and $I$ are independent of $N$. Thus, the overall
complexity of MTGAE is $\mathcal{O}(N)$, linear in the number of nodes.

##  2 Empirical Evaluation

Table 1: Summary of datasets (_left_) and baselines (_right_) used in
empirical evaluation. See [12, 14] for dataset details. The notation
$|\mathit{O}^{+}|$:$|\mathit{O}^{-}|$ denotes the ratio of positive to
negative edges and is a measure of class imbalance. Label rate is defined as
the number of nodes labeled for training divided by the total number of nodes.
Acronyms: AUC – Area Under ROC Curve; AP – Average Precision.

Dataset Nodes Average $|\mathit{O}^{+}|$:$|\mathit{O}^{-}|$ Node Node Label
Degree Ratio Features Classes Rate Pubmed 19,717 4.5 $1:4384$ 500 3 0.003
Citeseer 3,327 2.8 $1:1198$ 3,703 6 0.036 Cora 2,708 3.9 $1:\hskip 4.5pt694$
1,433 7 0.052 Arxiv-GRQC 5,242 5.5 $1:\hskip 4.5pt947$ – – – BlogCatalog
10,312 64.8 $1:\hskip 4.5pt158$ – – –

Baseline Evaluation Task Metric SDNE [14] Reconstruction Precision@$k$ VGAE
[8] Link Prediction AUC, AP GCN [7] Node Classification Accuracy

Implementation Details We closely follow the experimental protocols described
in [7, 8] to train and evaluate our MTGAE model for LPNC. For link prediction,
we form disjoint test, and validation, sets containing 10, and 5, percent of
randomly sampled positive links and the same number of negative links,
respectively, while utilizing all node features. For node classification, we
split the data into disjoint test, and validation, sets of 1,000, and 500,
examples, respectively and use only 20 examples per class for semi-supervised
learning. In comparison to the baselines, we evaluate our MTGAE model on the
same data splits over 10 runs with random weight initialization and report
mean AUC/AP scores for link prediction and accuracy scores for node
classification. We also compare the representation capacity of our MTGAE model
against the related autoencoder-based SDNE model on the network reconstruction
task. We use the ranking metric precision@$k$ to evaluate the model’s ability
to retrieve positive edges as part of the reconstruction.

Hyper-parameter tuning is performed on the validation set. Key hyper-
parameters include mini-batch size, dimensionality of the hidden layers, and
the percentage of dropout regularization [13]. In all experiments, the
dimensionality of the hidden layers in the MTGAE architecture is fixed at
$N$-256-128-256-$N$. We train for 100 epochs using Adam [6] gradient descent
with a fixed learning rate of $0.001$ on mini-batches of 64 examples.

We implement the MTGAE architecture using Keras [3] on top of the GPU-enabled
TensorFlow [1] backend. The diagonal elements of the adjacency matrix are set
to $1$ with the interpretation that every node is connected to itself. We
apply mean-variance normalization after each ReLU activation layer to help
improve link prediction performance, where it compensates for noise between
train and test instances by normalizing the activations to have zero mean and
unit variance. During training, we implement several regularization techniques
to mitigate overfitting, including dropout for highly sparse graphs and early
stopping as a form of regularization in time when the model shows signs of
overfitting on the validation set. We initialize weights according to the
Glorot scheme described in [4]. We do not apply weight decay regularization.

Results and Analysis Results of the reconstruction task for the Arxiv-GRQC and
BlogCatalog network datasets are illustrated in Figure 2. In comparison to
SDNE, we show that our MTGAE model achieves better precision@$k$ performance
for all $k$ values, up to $k=10,000$ for Arxiv-GRQC and $k=100,000$ for
BlogCatalog, when trained on the complete datasets. We also systematically
test the capacity of the MTGAE model to reconstruct the original networks when
up to 80 percent of the edges are randomly removed, akin to the link
prediction task. We show that the MTGAE model only gets worse precision@$k$
performance than SDNE on the Arxiv-GRQC dataset when more than 40 percent of
the edges are missing. On the BlogCatalog dataset, the MTGAE model achieves
better precision@$k$ performance than SDNE for large $k$ values even when 80
percent of the edges are missing at random. This experiment demonstrates the
superior representation capacity of our MTGAE model when compared to SDNE,
which is attributed to parameter sharing in the architecture.

Figure 2: Comparison of precision@$k$ performance between our MTGAE model and
the related autoencoder-based SDNE model for the reconstruction task on the
Arxiv-GRQC and BlogCatalog network datasets. The parameter $k$ indicates the
total number of retrieved edges.

Lastly, we report LPNC results obtained by our MTGAE model in the MTL
scenario. The model takes as input an incomplete graph with 10 percent of the
positive edges, and the same number of negative edges, missing at random and
all available node features to simultaneously predict labels for the nodes and
missing edges. Table 2 shows the efficacy of the MTGAE model for MTL when
compared against recent state-of-the-art _task-specific_ link prediction and
node classification models, which require the complete adjacency matrix as
input. For link prediction, MTGAE significantly outperforms the best VGAE
model on Cora and Citeseer. For node classification, MTGAE is the best
performing model on the Citeseer and Pubmed datasets, which have very low node
label rates.

Table 2: Comparison of LPNC performances between our MTGAE model and recent
state-of-the-art graph embedding methods. Link prediction performance is
reported as the combined average of AUC and AP scores. Accuracy is used for
node classification performance.

Method Cora Citeseer Pubmed Link Prediction MTGAE 0.946 0.949 0.944 VGAE [8]
0.920 0.914 0.965 Node Classification MTGAE 0.790 0.718 0.804 GCN [7] 0.815
0.703 0.790 Planetoid [15] 0.757 0.647 0.772

Future Work Further research will explore inductive reasoning on out-of-
network nodes and mitigate $\mathcal{O}(N)$ complexity for improved
scalability on large, dynamic graphs.

## References

  * [1] Abadi, M. et al.: TensorFlow: Large-scale machine learning on heterogeneous systems. https://github.com/tensorflow/tensorflow (2015)
  * [2] Caruana, R.: Multitask learning: a knowledge-based source of inductive bias. In _Proceedings of the Tenth International Conference on Machine Learning_ (1993)
  * [3] Chollet, F. et al.: Keras. https://github.com/keras-team/keras (2015)
  * [4] Glorot, X., Bengio, Y.: Understanding the difficulty of training deep feedforward neural networks. In _AISTATS_ 9 (2010)
  * [5] Hassan, K., Bahareh, B.: A survey on statistical relational learning. In _Proceedings of the 23rd Canadian Conference on Advances in Artificial Intelligence_ , 256–268 (2010)
  * [6] Kingma, D.P., Ba, J.L.: Adam: a method for stochastic optimization. In _International Conference on Learning Representations (ICLR)_ (2015)
  * [7] Kipf, T.N., Welling, M.: Semi-supervised classification with graph convolutional networks. _arXiv:1609.02907_ (2016)
  * [8] Kipf, T.N., Welling, M.: Variational graph auto-encoders. _arXiv:1611.07308_ (2016)
  * [9] Kuchaiev, O., Ginsburg, B.: Training deep autoencoders for collaborative filtering. _arXiv:1708.01715_ (2017)
  * [10] Rossi, R., McDowell, L., Aha, D., Neville, J.: Transforming graph data for statistical relational learning. _Journal of Artificial Intelligence Research_ 45, 363–441 (2012)
  * [11] Sedhain, S., Menon, A.K., Sanner, S., Xie, L.: AutoRec: Autoencoders meet collaborative filtering. In _Proceedings of the 24th International Conference on World Wide Web_ (2015)
  * [12] Sen, P., Namata, G.M., Bilgic, M., Getoor, L., Gallagher, B., Eliassi-Rad, T.: Collective classification in network data. _AI Magazine_ 29, 93–106 (2008)
  * [13] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R.: Dropout: A simple way to prevent neural networks from overfitting. _Journal of Machine Learning Research_ 15, 1929–1958 (2014)
  * [14] Wang, D., Cui, P., Zhu, W.: Structural deep network embedding. In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD)_ (2016)
  * [15] Yang, Z., Cohen, W., Salakhutdinov, R.: Revisiting semi-supervised learning with graph embeddings. In _International Conference on Machine Learning (ICML)_ (2016)
  * [16] Yang, Y., Hospedales, T.M.: Trace norm regularised deep multi-task learning. _arXiv:1606.04038_ (2017)