# Robust constraints on feebly interacting particles using XMM-Newton

Pedro De la Torre Luque<EMAIL_ADDRESS>The Oskar Klein
Centre, Department of Physics, Stockholm University, Stockholm 106 91, Sweden
Shyam Balaji<EMAIL_ADDRESS>Laboratoire de Physique Théorique et
Hautes Energies (LPTHE), UMR 7589 CNRS & Sorbonne Université, 4 Place Jussieu,
F-75252, Paris, France  Institut d’Astrophysique de Paris, UMR 7095 CNRS &
Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France  Pierluca
Carenza<EMAIL_ADDRESS>The Oskar Klein Centre, Department of
Physics, Stockholm University, Stockholm 106 91, Sweden

###### Abstract

During galactic Supernova (SN) explosions, a large amount of feebly
interacting particles (FIPs) may be produced. In this work we analyze
electrophilic FIPs with masses in the MeV-range that escape from SN and decay
into electron-positron pairs, causing an exotic leptonic injection. This
contribution adds up to known components, leading to an unexpected excess of
X-ray fluxes generated by inverse-Compton scattering of the injected particles
on low-energy photon backgrounds. For the first time in the context of FIPs,
we use XMM-Newton X-ray measurements to obtain the strongest and most robust
bounds on electrophilic FIPs produced by SN in our Galaxy.

_Electrophilic FIPs._ Massive stars explode in Supernovae (SN), reaching
nuclear densities and very high temperatures in the core. These extreme
environmental conditions are ideal to abundantly produce light Feebly
Interacting Particles (FIPs), if they exist. Currently, a lot of attention is
devoted to FIP phenomenology from the point of view of theoretical and
experimental investigation Antel  _et al._ (2023). The large family of FIPs
includes, among others, axions and axion-like particles (ALPs) Raffelt and
Seckel (1988); Keil  _et al._ (1997); Carenza  _et al._ (2019); Caputo  _et
al._ (2022), sterile neutrinos Kolb  _et al._ (1996); Raffelt and Zhou (2011);
Mastrototaro  _et al._ (2020), light CP-even scalars Balaji  _et al._ (2022,
2023) and dark photons (DPs) Chang  _et al._ (2017). Once produced in SN,
because of their weak interaction with matter, FIPs may escape the stellar
core and compete with neutrinos in subtracting energy to the SN. This in turn
leads to a shortening of the neutrino burst signal Raffelt (2012); Chang  _et
al._ (2018). Based on the SN 1987A neutrino observation, roughly in agreement
with a SN model cooling via only neutrino emission, it is possible to set
stringent constraints on FIP properties Raffelt and Seckel (1988). The
importance of SN for FIP phenomenology is not limited to indirect consequences
on the neutrino burst, but also for direct signatures of FIP interactions and
decays outside the star. In the case of ALPs, the phenomenology connected with
their conversions into photons in astrophysical magnetic fields, or radiative
decays, is well explored for the case of SN 1987A, future galactic or
extragalactic SN. For some examples, see Refs. Payez  _et al._ (2015); Hoof
and Schulz (2023); Meyer  _et al._ (2017); Meyer and Petrushevska (2020);
Crnogorčević  _et al._ (2021); Calore  _et al._ (2023, 2020, 2022a) for light
ALPs and Refs. Giannotti  _et al._ (2011); Jaeckel  _et al._ (2018); Hoof and
Schulz (2023); Müller  _et al._ (2023a, b) for heavy ones.

Figure 1: Limits on the production of electrons-positrons from FIP decays at
$95\%$ Confidence Level (C.L.). We compare the limits derived from the fit to
Voyager-1 data (blue), SPI latitude (green) and longitude (cyan) profiles,
from Ref. De la Torre Luque  _et al._ , with the ones obtained in this work
from observations of the MOS detector (red). The main systematic uncertainties
in the evaluation of these limits are represented as hatched bars, which, in
the case of MOS, correspond to uncertainties on electron/positron
reacceleration by plasma turbulence. We also evaluate the limit as obtained
for different values of the Alfvèn speed, ranging from $\sim 10$ km/s to $40$
km/s. These limits are valid for FIP masses below $\sim 20$ MeV.

In this work we analyze FIPs with masses in the MeV-range, produced in SN, and
decaying into electron-positron pairs outside the progenitor star. We refer to
these particles as ‘electrophilic’ Appelquist and Shrock (2002); Asaka  _et
al._ (2005); Abel  _et al._ (2008); Jaeckel and Ringwald (2010); Alekhin  _et
al._ (2016). Several motivated FIPs belong to this class of exotic particles,
such as ALPs, sterile neutrinos and DPs. The exotic injection of electrons and
positrons by this mechanism summed up with the standard channels, mainly
pulsar wind nebulae, SN remnants Ginzburg and Syrovatskii (1969); Gabici  _et
al._ (2019) and cosmic ray (CR) interactions with the interstellar gas De la
Torre Luque  _et al._ (2023), leads to important observational consequences
that we will explore. For instance, positrons would annihilate with background
electrons and contribute to the monochromatic signal at 511 keV Dar  _et al._
(1987). Requiring that positrons injected by FIP decays do not produce a
signal exceeding the one measured by the Spectrometer on INTEGRAL (SPI)
instrument Bouchet  _et al._ (2010); Siegert  _et al._ (2016), bounds on
various FIP models were determined DeRocco  _et al._ (2019); Calore  _et al._
(2021, 2022b).

We have explored, for the first time, a different phenomenology and probe of
FIPs. We characterize the inverse-Compton (IC) emission generated by the
injected leptons when scattering on the photon background (see Ref. Cirelli
_et al._ (2021) for an application to Weakly Interacting Massive Particles,
WIMPs). We use the X-ray observations from XMM-Newton in the range $2.5$-$8$
keV to place robust constraints on electrophilic FIPs, setting the stage for
future studies of these data to probe FIPs. As we see from Fig. 1, where we
show the upper limits in the number on electrons/positrons produced in the
Galaxy by SN-injected FIPs for different observables, XMM-Newton data (in
particular, from the MOS detector, in red) may lead to the most stringent
constraints up to date, even stronger than the emission line at $511$ keV
(longitudinal profile, cyan).

_FIP-induced electron-positron flux injection._ Electrophilic FIPs are defined
as FIPs coupling mainly to electrons and positrons. Importantly for our
discussion, when their mass exceeds $\sim 1$ MeV the, dominant decay channel
is into electron-positron pairs. Common examples are ALPs ($a$) decaying as
$a\to e^{+}e^{-}$ or sterile neutrinos, where a possible decay channel is
$\nu_{s}\to\nu_{\mu}e^{+}e^{-}$. In this second case, even though the FIP
produces other particles in the decay, we can model it in a similar way as any
electrophilic FIP. There are many other representative cases, like scalars,
Kaluza-Klein gravitons and supersymmetric particles. Therefore, electrophilic
FIPs are commonly found in Standard Model extensions and it motivates our
discussion on their phenomenology.

Figure 2: Comparison of the predicted X-ray signals generated from
electrophilic FIPs from $1$ to $8$ keV, compared to the data from the MOS
detector (XMM-Newton) for the rings 3 (the most constraining), in the left
panel, and ring 30 (the less constraining), in the right panel. We show the
predicted FIP X-ray emission for FIP masses ranging from $2$ to $20$ MeV,
where FIPs produce roughly the same emission. The region below $2.5$ keV, not
considered in our analysis, is shaded.

After being generated during galactic SN explosions, FIPs escape the parent
star and undergo a decay process that populates the diffuse electron-positron
Galactic background. The spatial distribution of this injection roughly
follows the SN distribution, being affected by propagation effects. The
electron/positron flux caused by FIPs has an energy distribution inherited by
the parent FIP. Assuming that electrons/positrons have an energy $E_{e}$ which
is half of the decaying FIP, the injected flux can be parametrized as van
Bibber  _et al._ (1989)

$\begin{split}\frac{dN_{e}}{dE_{e}}&=N_{e}C_{0}\left(\frac{4E_{e}^{2}-m_{X}^{2}}{E_{0}^{2}}\right)^{\beta/2}e^{-(1+\beta)\frac{2E_{e}}{E_{0}}}\,,\\\
C_{0}&=\frac{2\sqrt{\pi}\left(\frac{1+\beta}{2m_{X}}\right)^{\frac{1+\beta}{2}}E_{0}^{\frac{\beta-1}{2}}}{K_{\frac{1+\beta}{2}}\left((1+\beta)\frac{m_{X}}{E_{0}}\right)\Gamma\left(1+\frac{\beta}{2}\right)}\,,\end{split}$
(1)

where FIPs are emitted with a modified blackbody spectrum with $E_{0}$
parameter related to the FIP average energy, a FIP mass $m_{X}>2m_{e}$ and a
spectral index $\beta$. In the normalization factor, $K_{\frac{1+\beta}{2}}$
is the modified Bessel function of the second kind of order $(1+\beta)/2$ and
$\Gamma$ is the Euler-Gamma function. This flux is normalized such that

$\int_{m_{X}/2}^{\infty}dE_{e}\frac{dN_{e}}{dE_{e}}=N_{e}\,,$ (2)

where with $N_{e}$ we denote the number of electrons, equal to the number of
positrons, produced in SN explosions by FIP decays, i.e.
$N_{e}=N_{e^{+}}=N_{e^{-}}$. Note that the simple prescription in Eq. (1) does
not depend on the type of FIP model, and the examples that we mentioned can be
parametrized in this model independent way.

The lepton fluxes injected by FIPs are affected by propagation effects in the
Galaxy. This is taken into account by solving the diffusion equation of these
particles with the DRAGON2
111https://github.com/cosmicrays/DRAGON2-Beta_version code Evoli  _et al._
(2017, 2018), with the assumption that FIPs decay close to their source, at a
distance smaller than a few kpc Calore  _et al._ (2021). Thus, we simulate the
leptonic injection as happening directly at the SN with an energy spectrum
defined as in Eq. (1) and following the spatial SN distribution evaluated in
Ref. Ferriere  _et al._ (2007) convolved with the Steiman-Cameron distribution
Steiman-Cameron  _et al._ (2010) of the spiral arms (four-arm model).

The injected flux, i.e. the number of electrons/positrons per unit area, time
and energy, is written as

$\frac{d\Phi_{e}}{dE_{e}}=\frac{\Gamma_{cc}}{4\pi
d_{SN}^{2}}\frac{dN_{e}}{dE_{e}}\,,$ (3)

where $\Gamma_{cc}=2$ SN per century is the galactic rate of core-collapse SN
explosions Rozwadowska  _et al._ (2021) and $d_{SN}=10.2$ kpc is an effective
length-scale that allows us to calculate the total area of emission from the
SN and is roughly the average distance between galactic SN calculated from the
source distribution employed in this work. This parameter is obtained by
imposing that the injected number of particles per unit of energy in Eq. (1)
is equal to the integral of the total flux density of particles (obtained
convolving Eq. (1) with the spatial distribution of sources) over the volume
of the Galaxy. Regarding the propagation, we employ a spatially-constant
diffusion coefficient derived from the analyses of secondary-to-primary flux
ratios in Ref. de la Torre Luque  _et al._ (2022) and adapted in Ref. De la
Torre Luque  _et al._ (2023) for the spiral arm structure of the source, gas
and magnetic field galactic distributions (B/C best-fit model). In particular,
for the diffusion coefficient defined in Eq. (3.3) of Ref. de la Torre Luque
_et al._ (2022), the parameters are shown in Tab. 1. We adopt the magnetic
field model derived by Ref. Pshirkov  _et al._ (2011), with a normalization of
the disk, halo and turbulent magnetic field intensities set to the values
found in Ref. Di Bernardo  _et al._ (2013) from the study of synchrotron
radiation. The energy density distribution of the radiation fields has been
taken from Ref. Porter  _et al._ (2008). A key point to assume a stationary
and smooth lepton injection is the observation that electrons and positrons
propagate and interact in the galactic environment on a timescale of Gyr,
extremely long compared to the SN explosions rate. Therefore, we can model the
lepton injection as time-independent and smoothly following the SN
distribution.

Norm. Energy | $E_{0}$ |  $4$ GeV
---|---|---
Diffusion coeff. | $D_{0}$ |  $1.02\times 10^{29}$ cm2 s-1
Diffusion index | $\delta$ | $0.49$
Break energy | $E_{b}$ |  $312$ GeV
Index break | $\Delta\delta$ | $0.20$
Smooth. param. | $s$ | $0.04$
$\beta$ exponent | $\eta$ | $-0.75$
Halo height | $H$ |  $8$ kpc
Alfvèn velocity | $V_{A}$ |  $13.4$ km s-1
Table 1: Main propagation parameters used in our analysis De la Torre Luque
_et al._ (2023).

Our results are robust at the level of tens of percent with respect to
uncertainties in the source distribution, parameters in the diffusion
coefficient and Galactic halo height, the main uncertainties come from the
determination of the Alfvèn speed, which constitutes an important variable for
our predictions. The Alfvèn speed is the main parameter characterizing the
reacceleration of charged particles when interacting with plasma turbulences
in the Galaxy Heinbach and Simon (1995); Strong and Moskalenko (1998);
Cardillo (2019); Bustard and Oh (2022). Different CR analyses obtain very
different values of this parameter, depending on the spallation cross sections
employed or the datasets analyzed, ranging anywhere from V${}_{A}\sim 0$ km/s
to V${}_{A}\sim 40$ km/s de la Torre Luque  _et al._ (2022); Luque  _et al._
(2021); Luque (2021). This makes our predicted signals uncertain by orders of
magnitude around the GeV scale. We have shown this in Fig. 8 of our companion
work De la Torre Luque  _et al._ . Something similar was observed in Ref.
Boudaud  _et al._ (2017) for the case of low mass WIMPs annihilating or
decaying via electrophilic channels. The effect of reacceleration in the
electron/positron signals produced by these particles was shown to cause
differences of several orders of magnitude in the extracted limits on the
annihilation rate for different reacceleration setups.

_XMM-Newton constraints on the exotic lepton injection._ During their
propagation in the Milky Way the electrons/positrons injected by FIPs interact
with the interstellar gas, the interstellar radiation fields (ISRFs) – which
mainly consist of the optical and ultraviolet light from stars, infrared light
from dust and the CMB – and the galactic magnetic field. The IC emission
produced from the boost of the low energy ISRF photons produce hard X-rays, at
keV scale, and low-energy $\gamma$-rays, at the MeV scale. For example, a
typical IR photon with an energy around $0.1$ eV would gain an energy factor
of $(E_{e}/m_{e})^{2}\sim 10^{4}$ when interacting with an electron of $E\sim
50$ MeV, acquiring an energy around the keV-scale. At slightly higher
energies, the IR field dominates the emission, and then the optical and UV
fields (with smaller energy density though) take over. As we show in Fig. 4 of
our companion work De la Torre Luque  _et al._ , one can benefit from MeV
$\gamma$-ray observations as well as X-ray observations in order to probe
these signals and the injected electron/positron population, setting
constraints on the injected number of electrons $N_{e}$.

In this study, we model the secondary diffuse emission due to the interaction
of the injected leptons using the HERMES code Dundovic  _et al._ (2021), which
performs a numerical integration of the emission along the line of sight at
each location in the Galaxy and energy, using detailed gas maps and updated
ISRF models Vernetto and Lipari (2021, 2016). The effect of reacceleration in
the electron spectrum causes a significant increase in the X-ray emission from
the IC radiation of ISRF photons. In this work, we show how reacceleration
boosts these signals and we derive constraints on the injection of
electrons/positrons using data from the MOS detector Turner  _et al._ (2001)
on board the XMM-Newton telescope Jansen  _et al._ (2001).

X-ray observations allow us to constrain these signals, benefiting from lower
systematic uncertainties. In particular, we use the data from the MOS detector
in the energy range from $1$ keV to above $10$ keV, which was provided by Ref.
Foster  _et al._ (2021) divided in galactocentric rings around the Galactic
center that are $6^{\circ}$ wide. We observe that our constraints are stronger
towards the inner rings and weaker in the external rings. However, the limits
derived from the two inner most rings (rings 1 and 2) are weaker than the
limit derived from ring 3 (the strongest limit), as it was observed in Ref.
Cirelli  _et al._ (2023). In Fig. 2 we compare the calculated FIP signals, for
FIP masses between $2$ and $20$ MeV, with MOS data at the most constraining
ring (ring 3) and the outermost ring (ring 30). The X-ray flux produced from
FIPs of such masses differ by no more than $\sim 10\%$, whereas the energy
trend of the emission is roughly the same. As we see, while an injection of
electrons of $N_{e}=5\times 10^{53}$ already lies above MOS data in the
$2$-$4$ keV range for ring 3 (left panel), the outermost ring only allows us
to constrain $N_{e}$ to be smaller than $\sim 10^{56}$. In addition, the
spread of the data for the outer-most ring (ring 30) makes it much less
constraining than the inner rings. In particular, we derive bounds on $N_{e}$
at $95\%$ C.L. from a $\chi^{2}$ analysis of the data from $2.5$ to $8$ keV,
as performed by Refs. Foster  _et al._ (2021); Cirelli  _et al._ (2023), due
to the instrumental noise in the detector at lower and higher energies
respectively. In Tab. 2 we show the limits obtained for a $m_{X}=10$ MeV FIP
for different rings and Alfvèn velocities. On top of this, combining data and
predictions for all the rings allow us to set a constraint that a factor $\sim
2$ stronger than the one obtained from the ring 3 alone.

In Fig. 3, we show the predicted X-ray signal generated from the
electrons/positrons produced by FIPs for different reacceleration levels,
corresponding to Alfvèn speeds of $13$ (our reference value), $20$, $30$ and
$40$ km/s. These are the fluxes predicted for an injection of $N_{e}=2.5\times
10^{52}$, which is roughly the most optimistic limit obtained from the
analysis of the longitude profile of the $511$ keV line emission. This already
shows that for an Alfvèn speed above $\sim 30$ km/s, the predicted X-ray
emission will overproduce the MOS data for ring 3.

Figure 3: Similar to what is shown in Fig. 2 for different levels of
reacceleration, corresponding to Alfvèn speeds of $13$, $20$, $30$ and $40$
km/s, for a FIP mass $m_{X}=10$ MeV and $N_{e}=2.5\times 10^{52}$. These
limits are compared to the data from ring 3 of the MOS detector.

|  $N_{e}(\times 10^{53})$ $95\%$ C.L.
---|---
Best-fit VA (ring 1) | $37$
Best-fit VA (ring 3) | $4.7$
Best-fit VA (ring 4) | $6.1$
Best-fit VA (ring 30) | $1015$
Best-fit VA (Combined) | $2.8$
V${}_{A}=20$ km/s (Comb.) | $1.1$
V${}_{A}=30$ km/s (Comb.) | $0.18$
V${}_{A}=40$ km/s (Comb.) | $0.056$

Table 2: Limits on $N_{e}$ at the $95\%$ C.L. derived from MOS data between
$2.5$ keV and $8$ keV. Results for various rings and Alfvèn velocities are
shown.

This, in addition to the fact that MOS measurements and its uncertainties are
much more robust than those evaluated for the $511$ keV line makes the MOS
instrument exceptionally valuable for probing FIP physics.

Finally, we may comment again on Fig. 1, where we compare the limits obtained
for the number of electrons/positrons injected by FIPs with masses ranging
from $\sim 1$ MeV to $\sim 20$ MeV from different observables. We added
hatched regions indicating the uncertainties in our predictions (see Ref. De
la Torre Luque  _et al._ for a complete description of the Voyager-1 and SPI
limits). We also indicate the limit obtained for different values of the
Alfvèn speed, ranging from $\sim 10$ km/s to $40$ km/s. Notably, the data from
the MOS detector (combining all the rings) leads to constraints on $N_{e}$
ranging from $\sim 3\times 10^{53}$ to $\sim 5\times 10^{53}$, which are
compatible and even better, in the case of high reacceleration, to those
derived from the longitude profile of the $511$ keV line. With the latter
being less robust than MOS constraints due to sizeable systematic
uncertainties in the extraction of the data, as discussed in more detail in
Ref. De la Torre Luque  _et al._ .

_XMM-Newton constraints on FIP properties._ In Ref. De la Torre Luque  _et
al._ we derived constraints on specific FIP models, corresponding to our
bounds on $N_{e}$. Here, we briefly discuss how we expect the constraints on
FIPs to improve with the proposed bounds and we refer the interested reader to
Ref. De la Torre Luque  _et al._ for a comprehensive discussion.

An useful simplification is to consider constraints in the weak coupling
regime, i.e. where the positron production per SN grows as $g^{2}$ with $g$
FIP-electron coupling. Thus, we neglect the strong coupling regime, with an
exponential sensitivity to $g$ since FIPs starts to decay inside the SN and
part of the produced electron-positron pairs remains trapped inside the star.
Under this assumption, the number of electrons/positrons per SN is
proportional to some power of the coupling, $N_{e}\sim g^{\alpha}$, where
$\alpha$ might be different from $2$ in order to schematically take into
account the effect of type Ib/c SNe. Precisely, in this discussion we will
consider sterile neutrinos mixed with muon and tau neutrinos, and DPs, all of
them with a mass $m_{X}=10$ MeV. In these cases, $g=|U_{s\mu}|$ (sterile-muon
neutrino mixing parameter) and $\alpha=1.64$, $g=|U_{s\tau}|$ and
$\alpha=1.52$, $g=\epsilon$ (kinetic mixing parameter with the photon) and
$\alpha=2$, respectively De la Torre Luque  _et al._ . In the case of ALPs we
can only determine a bound in the strongly coupled regime because of our
assumption of FIPs decaying close to the production point, that is violated in
the weak coupling regime. Thus, we do not discuss ALP constraints here.

Now it is easy to understand how the limiting values of $N_{e}$ can be
reflected as constraints on the FIP-electron couplings: $g\sim
N_{e}^{1/\alpha}$. In Fig. 4 we show the improvements of the bound on $g$ as a
function of $N_{e}$, with respect to a nominal bound placed with
$N_{e}=10^{55}$, for sterile neutrinos (black and red lines) and DPs (blue
line). In Fig. 4 we highlight as colored regions, the bounds discussed in this
work with the same color code of Fig. 1, showing only the most stringent
constraints related to the 511 keV line. The red hatched region is the
uncertainty band of the MOS constraint that, under optimistic assumptions,
improves the leading one by a factor $\sim 3$ on $N_{e}$. In the considered
models, bounds on the FIP couplings improve up to a factor $\sim 2$ thanks to
the MOS analysis.

_Conclusions._ In conclusion, we have explored for the first time, the
capability of the hard X-ray observations from the XMM-Newton experiment to
set constraints on the production of electrons/positrons from the decay of
FIPs. In addition, we evaluated the impact of reacceleration on the
electromagnetic signals generated from FIPs. We find that observations from
the XMM-Newton mission allow us to set one of the strongest limits on the
production of electrons and positrons from FIPs in the Galaxy, possibly even
more stringent than those derived from the $511$ keV emission. This study
opens the path for future searches of particles beyond Standard Model coupled
with electrons and positrons in astrophysical X-ray spectra, showing a
remarkable complementary with cosmological bounds.

Figure 4: Improvement of the bounds on $g$, for various FIP models discussed
in the text, as a function of the maximal $N_{e}$ allowed by observations. The
511 keV line bounds, for the latitudinal and longitudinal distributions are
shown in green and cyan, respectively. The Voyager-1 bound is shown in blue
and the red hatched region refers to the uncertainty band on the maximal
$N_{e}$ associated with MOS observations.

###### Acknowledgements.

_Acknowledgements._ SB would like to thank Jordan Koechler and Marco Cirelli
for useful discussions regarding XMM-Newton data. This article is based upon
work from COST Action COSMIC WISPers CA21106, supported by COST (European
Cooperation in Science and Technology). SB is supported by funding from the
European Union’s Horizon 2020 research and innovation programme under grant
agreement No. 101002846 (ERC CoG “CosmoChart”) as well as support from the
Initiative Physique des Infinis (IPI), a research training program of the Idex
SUPER at Sorbonne Université. The work of PC is supported by the European
Research Council under Grant No. 742104 and by the Swedish Research Council
(VR) under grants 2018-03641 and 2019-02337. PDL is supported by the European
Research Council under grant 742104 and the Swedish National Space Agency
under contract 117/19. This project used computing resources from the Swedish
National Infrastructure for Computing (SNIC) projects 2021/3-42, 2021/6-326,
2021-1-24 and 2022/3-27 partially funded by the Swedish Research Council
through grant no. 2018-05973.

## References

  * Antel  _et al._ (2023) C. Antel _et al._ , in  _Workshop on Feebly-Interacting Particles_ (2023) arXiv:2305.01715 [hep-ph] .
  * Raffelt and Seckel (1988) G. Raffelt and D. Seckel, Phys. Rev. Lett. 60, 1793 (1988).
  * Keil  _et al._ (1997) W. Keil, H.-T. Janka, D. N. Schramm, G. Sigl, M. S. Turner, and J. R. Ellis, Phys. Rev. D 56, 2419 (1997), arXiv:astro-ph/9612222 .
  * Carenza  _et al._ (2019) P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo, and A. Mirizzi, JCAP 10, 016 (2019), [Erratum: JCAP 05, E01 (2020)], arXiv:1906.11844 [hep-ph] .
  * Caputo  _et al._ (2022) A. Caputo, G. Raffelt, and E. Vitagliano, Phys. Rev. D 105, 035022 (2022), arXiv:2109.03244 [hep-ph] .
  * Kolb  _et al._ (1996) E. W. Kolb, R. N. Mohapatra, and V. L. Teplitz, Phys. Rev. Lett. 77, 3066 (1996), arXiv:hep-ph/9605350 .
  * Raffelt and Zhou (2011) G. G. Raffelt and S. Zhou, Phys. Rev. D 83, 093014 (2011), arXiv:1102.5124 [hep-ph] .
  * Mastrototaro  _et al._ (2020) L. Mastrototaro, A. Mirizzi, P. D. Serpico, and A. Esmaili, JCAP 01, 010 (2020), arXiv:1910.10249 [hep-ph] .
  * Balaji  _et al._ (2022) S. Balaji, P. S. B. Dev, J. Silk, and Y. Zhang, JCAP 12, 024 (2022), arXiv:2205.01669 [hep-ph] .
  * Balaji  _et al._ (2023) S. Balaji, M. E. Ramirez-Quezada, J. Silk, and Y. Zhang, Phys. Rev. D 107, 083038 (2023), arXiv:2301.02258 [hep-ph] .
  * Chang  _et al._ (2017) J. H. Chang, R. Essig, and S. D. McDermott, JHEP 01, 107 (2017), arXiv:1611.03864 [hep-ph] .
  * Raffelt (2012) G. G. Raffelt, Proc. Int. Sch. Phys. Fermi 182, 61 (2012), arXiv:1201.1637 [astro-ph.SR] .
  * Chang  _et al._ (2018) J. H. Chang, R. Essig, and S. D. McDermott, JHEP 09, 051 (2018), arXiv:1803.00993 [hep-ph] .
  * Payez  _et al._ (2015) A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi, and A. Ringwald, JCAP 02, 006 (2015), arXiv:1410.3747 [astro-ph.HE] .
  * Hoof and Schulz (2023) S. Hoof and L. Schulz, JCAP 03, 054 (2023), arXiv:2212.09764 [hep-ph] .
  * Meyer  _et al._ (2017) M. Meyer, M. Giannotti, A. Mirizzi, J. Conrad, and M. A. Sánchez-Conde, Phys. Rev. Lett. 118, 011103 (2017), arXiv:1609.02350 [astro-ph.HE] .
  * Meyer and Petrushevska (2020) M. Meyer and T. Petrushevska, Phys. Rev. Lett. 124, 231101 (2020), [Erratum: Phys.Rev.Lett. 125, 119901 (2020)], arXiv:2006.06722 [astro-ph.HE] .
  * Crnogorčević  _et al._ (2021) M. Crnogorčević, R. Caputo, M. Meyer, N. Omodei, and M. Gustafsson, Phys. Rev. D 104, 103001 (2021), arXiv:2109.05790 [astro-ph.HE] .
  * Calore  _et al._ (2023) F. Calore, P. Carenza, C. Eckner, M. Giannotti, G. Lucente, A. Mirizzi, and F. Sivo, (2023), arXiv:2306.03925 [astro-ph.HE] .
  * Calore  _et al._ (2020) F. Calore, P. Carenza, M. Giannotti, J. Jaeckel, and A. Mirizzi, Phys. Rev. D 102, 123005 (2020), arXiv:2008.11741 [hep-ph] .
  * Calore  _et al._ (2022a) F. Calore, P. Carenza, C. Eckner, T. Fischer, M. Giannotti, J. Jaeckel, K. Kotake, T. Kuroda, A. Mirizzi, and F. Sivo, Phys. Rev. D 105, 063028 (2022a), arXiv:2110.03679 [astro-ph.HE] .
  * Giannotti  _et al._ (2011) M. Giannotti, L. D. Duffy, and R. Nita, JCAP 01, 015 (2011), arXiv:1009.5714 [astro-ph.HE] .
  * Jaeckel  _et al._ (2018) J. Jaeckel, P. C. Malta, and J. Redondo, Phys. Rev. D 98, 055032 (2018), arXiv:1702.02964 [hep-ph] .
  * Müller  _et al._ (2023a) E. Müller, F. Calore, P. Carenza, C. Eckner, and M. C. D. Marsh, (2023a), arXiv:2304.01060 [astro-ph.HE] .
  * Müller  _et al._ (2023b) E. Müller, P. Carenza, C. Eckner, and A. Goobar, (2023b), arXiv:2306.16397 [astro-ph.HE] .
  * (26) P. De la Torre Luque, S. Balaji, and P. Carenza, arXiv:23XX.XXXX [hep-ph] .
  * Appelquist and Shrock (2002) T. Appelquist and R. Shrock, Phys. Lett. B 548, 204 (2002), arXiv:hep-ph/0204141 .
  * Asaka  _et al._ (2005) T. Asaka, S. Blanchet, and M. Shaposhnikov, Phys. Lett. B 631, 151 (2005), arXiv:hep-ph/0503065 .
  * Abel  _et al._ (2008) S. A. Abel, M. D. Goodsell, J. Jaeckel, V. V. Khoze, and A. Ringwald, JHEP 07, 124 (2008), arXiv:0803.1449 [hep-ph] .
  * Jaeckel and Ringwald (2010) J. Jaeckel and A. Ringwald, Ann. Rev. Nucl. Part. Sci. 60, 405 (2010), arXiv:1002.0329 [hep-ph] .
  * Alekhin  _et al._ (2016) S. Alekhin _et al._ , Rept. Prog. Phys. 79, 124201 (2016), arXiv:1504.04855 [hep-ph] .
  * Ginzburg and Syrovatskii (1969) V. L. Ginzburg and S. I. Syrovatskii,  _The origin of cosmic rays_ (1969).
  * Gabici  _et al._ (2019) S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, E. Orlando, A. Strong, and A. Vittino, Int. J. Mod. Phys. D 28, 1930022 (2019), arXiv:1903.11584 [astro-ph.HE] .
  * De la Torre Luque  _et al._ (2023) P. De la Torre Luque, M. N. Mazziotta, and F. Loparco, (2023), arXiv:2305.02958 [astro-ph.HE] .
  * Dar  _et al._ (1987) A. Dar, J. Goodman, and S. Nussinov, Phys. Rev. Lett. 58, 2146 (1987), [Erratum: Phys.Rev.Lett. 59, 1871 (1987)].
  * Bouchet  _et al._ (2010) L. Bouchet, J.-P. Roques, and E. Jourdain, Astrophys. J. 720, 1772 (2010), arXiv:1007.4753 [astro-ph.HE] .
  * Siegert  _et al._ (2016) T. Siegert, R. Diehl, G. Khachatryan, M. G. H. Krause, F. Guglielmetti, J. Greiner, A. W. Strong, and X. Zhang, Astron. Astrophys. 586, A84 (2016), arXiv:1512.00325 [astro-ph.HE] .
  * DeRocco  _et al._ (2019) W. DeRocco, P. W. Graham, D. Kasen, G. Marques-Tavares, and S. Rajendran, JHEP 02, 171 (2019), arXiv:1901.08596 [hep-ph] .
  * Calore  _et al._ (2021) F. Calore, P. Carenza, M. Giannotti, J. Jaeckel, G. Lucente, and A. Mirizzi, Phys. Rev. D 104, 043016 (2021), arXiv:2107.02186 [hep-ph] .
  * Calore  _et al._ (2022b) F. Calore, P. Carenza, M. Giannotti, J. Jaeckel, G. Lucente, L. Mastrototaro, and A. Mirizzi, Phys. Rev. D 105, 063026 (2022b), arXiv:2112.08382 [hep-ph] .
  * Cirelli  _et al._ (2021) M. Cirelli, N. Fornengo, B. J. Kavanagh, and E. Pinetti, Phys. Rev. D 103, 063022 (2021), arXiv:2007.11493 [hep-ph] .
  * van Bibber  _et al._ (1989) K. van Bibber, P. M. McIntyre, D. E. Morris, and G. G. Raffelt, Phys. Rev. D 39, 2089 (1989).
  * Evoli  _et al._ (2017) C. Evoli, D. Gaggero, A. Vittino, G. Di Bernardo, M. Di Mauro, A. Ligorini, P. Ullio, and D. Grasso, JCAP 02, 015 (2017), arXiv:1607.07886 [astro-ph.HE] .
  * Evoli  _et al._ (2018) C. Evoli, D. Gaggero, A. Vittino, M. Di Mauro, D. Grasso, and M. N. Mazziotta, JCAP 07, 006 (2018), arXiv:1711.09616 [astro-ph.HE] .
  * Ferriere  _et al._ (2007) K. Ferriere, W. Gillard, and P. Jean, Astron. Astrophys. 467, 611 (2007), arXiv:astro-ph/0702532 .
  * Steiman-Cameron  _et al._ (2010) T. Y. Steiman-Cameron, M. Wolfire, and D. Hollenbach, Astrophys. J. 722, 1460 (2010).
  * Rozwadowska  _et al._ (2021) K. Rozwadowska, F. Vissani, and E. Cappellaro, New Astron. 83, 101498 (2021), arXiv:2009.03438 [astro-ph.HE] .
  * de la Torre Luque  _et al._ (2022) P. de la Torre Luque, M. N. Mazziotta, A. Ferrari, F. Loparco, P. Sala, and D. Serini, JCAP 07, 008 (2022), arXiv:2202.03559 [astro-ph.HE] .
  * Pshirkov  _et al._ (2011) M. S. Pshirkov, P. G. Tinyakov, P. P. Kronberg, and K. J. Newton-McGee, Astrophys. J. 738, 192 (2011), arXiv:1103.0814 [astro-ph.GA] .
  * Di Bernardo  _et al._ (2013) G. Di Bernardo, C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, JCAP 03, 036 (2013), arXiv:1210.4546 [astro-ph.HE] .
  * Porter  _et al._ (2008) T. A. Porter, I. V. Moskalenko, A. W. Strong, E. Orlando, and L. Bouchet, Astrophys. J. 682, 400 (2008), arXiv:0804.1774 [astro-ph] .
  * Heinbach and Simon (1995) U. Heinbach and M. Simon, Astrophys. J. 441, 209 (1995).
  * Strong and Moskalenko (1998) A. W. Strong and I. V. Moskalenko, Astrophys. J. 509, 212 (1998), arXiv:astro-ph/9807150 [astro-ph] .
  * Cardillo (2019) M. Cardillo, Galaxies 7, 49 (2019), arXiv:1904.12476 [astro-ph.HE] .
  * Bustard and Oh (2022) C. Bustard and S. P. Oh, Astrophys. J. 941, 65 (2022), arXiv:2208.02261 [astro-ph.HE] .
  * Luque  _et al._ (2021) P. D. L. T. Luque, M. N. Mazziotta, F. Loparco, F. Gargano, and D. Serini, JCAP 07, 010 (2021), arXiv:2102.13238 [astro-ph.HE] .
  * Luque (2021) P. D. L. T. Luque, JCAP 11, 018 (2021), arXiv:2107.06863 [astro-ph.HE] .
  * Boudaud  _et al._ (2017) M. Boudaud, J. Lavalle, and P. Salati, Phys. Rev. Lett. 119, 021103 (2017), arXiv:1612.07698 [astro-ph.HE] .
  * Dundovic  _et al._ (2021) A. Dundovic, C. Evoli, D. Gaggero, and D. Grasso, Astron. Astrophys. 653, A18 (2021), arXiv:2105.13165 [astro-ph.HE] .
  * Vernetto and Lipari (2021) S. Vernetto and P. Lipari, PoS ICRC2021, 923 (2021).
  * Vernetto and Lipari (2016) S. Vernetto and P. Lipari, Phys. Rev. D 94, 063009 (2016), arXiv:1608.01587 [astro-ph.HE] .
  * Turner  _et al._ (2001) M. J. L. Turner _et al._ , Astron. Astrophys. 365, L27 (2001), arXiv:astro-ph/0011498 .
  * Jansen  _et al._ (2001) F. Jansen _et al._ (XMM), Astron. Astrophys. 365, L1 (2001).
  * Foster  _et al._ (2021) J. W. Foster, M. Kongsore, C. Dessert, Y. Park, N. L. Rodd, K. Cranmer, and B. R. Safdi, Phys. Rev. Lett. 127, 051101 (2021), arXiv:2102.02207 [astro-ph.CO] .
  * Cirelli  _et al._ (2023) M. Cirelli, N. Fornengo, J. Koechler, E. Pinetti, and B. M. Roach, JCAP 07, 026 (2023), arXiv:2303.08854 [hep-ph] .


# Coordinating Cross-modal Distillation for Molecular Property Prediction
††thanks: * These authors contributed equally.

Hao Zhang2,* Nan Zhang1,*, Ruixin zhang2, Lei Shen2, Yingyi Zhang2, Meng Liu3
1 Academy for Engineering and Technology, Fudan University, China
2 Tencent Youtu Lab
3 Shanghai Jiao Tong University

###### Abstract

In recent years, molecular graph representation learning (GRL) has drawn much
more attention in molecular property prediction (MPP) problems. The existing
graph methods have demonstrated that 3D geometric information is significant
for better performance in MPP. However, accurate 3D structures are often
costly and time-consuming to obtain, limiting the large-scale application of
GRL. It is an intuitive solution to train with 3D to 2D knowledge distillation
and predict with only 2D inputs. But some challenging problems remain open for
3D to 2D distillation. One is that the 3D view is quite distinct from the 2D
view, and the other is that the gradient magnitudes of atoms in distillation
are discrepant and unstable due to the variable molecular size. To address
these challenging problems, we exclusively propose a distillation framework
that contains global molecular distillation and local atom distillation. We
also provide a theoretical insight to justify how to coordinate atom and
molecular information, which tackles the drawback of variable molecular size
for atom information distillation. Experimental results on two popular
molecular datasets demonstrate that our proposed model achieves superior
performance over other methods. Specifically, on the largest MPP dataset
PCQM4Mv2 served as an ”ImageNet Large Scale Visual Recognition Challenge” in
the field of graph ML, the proposed method achieved a 6.9% improvement
compared with the best works. And we obtained fourth place with the MAE of
0.0734 on the test-challenge set for OGB-LSC 2022 Graph Regression Task. We
will release the code soon.

##  I Introduction

In recent years, molecular property prediction (MPP) has played a crucial role
in new materials and drug discovery. Since atoms and bonds in a molecule are
naturally in a graph structure, graph neural networks (GNN) [1] have shown
great vitality in MPP. Many works [2, 3, 4] take 2D molecular structures as
input of GNN for training and testing. Their performances are limited because
they do not consider the 3D atom coordinates that determine certain chemical
and physical functionalities of molecules. Methods [5, 6, 7] that exploit 3D
structures exhibit advantages in accuracy, especially properties related to
quantum mechanics, such as single-point energy, atomic forces, and dipole
moments. However, obtaining accurate 3D structures is expensive and time-
consuming. Calculating 3D information with DFT takes ${10^{5}}$ times longer
than model feedforward [2], which seriously impedes the application of
3D-based methods.

Some works attempt to augment 2D information with 3D information through
training to enhance the prediction accuracy while keeping inference speed. One
direction is developing 2D-3D contrastive tasks. Another way is to generate
original conformations through reconstruction tasks. Usually, these tasks
require several different 3D conformers to create positive and negative pairs
for self-supervised learning, but obtaining different conformers is time-
consuming for large-scale data. Also, different conformers from one SMILES[8]
code may represent different property values, which lower the performance of
contrastive-based methods. Lastly, some applications do not support the
requirement of multi conformers. With the force field given, one SMILES code
usually can only generate one 3D conformer structure.

Distillation is another intuitive solution that utilizes 3D knowledge in
training while keeping inference speed. Also, distillation methods do not have
the limitation of multiple conformers requirements. But some challenging
problems remain for distillation methods. On the one hand, the 3D view is
quite distinct from the 2D view, making the top-level-feature-only
distillation less effective. On the other hand, which tokens participate in
distillation is still an open question for the MPP. As self-attention
calculates information of the virtual token from all-atom features, it is
natural to include all-atom tokens together in distillation. But simply
distilling atom tokens and virtual token together leads to poor performance.
We found the gradient magnitudes of atom tokens distillation are inconsistent
and unstable due to the variable molecular size. So distilling local atom
information without correction will cause degradation problems.

Motivated by the above observation, we propose a novel coordinating cross-
modal distillation (CCMD) framework. We formulate a coordinating weight, which
dynamically scales and balances the gradients of global molecular token and
local atom tokens according to the number of atoms. Also, the distillation of
virtual token and atom tokens is carried in all layers to enhance the 2D-3D
distillation performance. Lastly, we introduce absolute position encoding to
integrate adjacent edge information, achieving better performance in both
views of 2D and 3D. Experimental results on two popular molecular datasets
demonstrate that our proposed model achieves state-of-the-art performance.

In conclusion, our main contributions to this work can be summarized as
follows:

  * •

We exclusively develop an enhanced cross-modal distillation on all layers for
molecular property prediction, which distills global molecular information and
local atom information to transfer 3D geometric information to 2D view.

  * •

For coordinating global molecular and local atom information, we formulate a
novel weight according to the number of atoms in a molecular to scale gradient
magnitudes, which boosts the distillation performance. Specific, we
demonstrate the coordinating weight is $f(N)=\frac{1}{{N}^{2}}$ for
transformer and $f(N)=\frac{1}{{N}}$ for GIN. N is the number of atoms in a
molecular.

  * •

We design an absolute position encoding integrated with adjacent edge
information, which could distinguish atoms same atom id in the input layer and
increase the performance of both 2D and 3D.

  * •

CCMD can enable molecular attention to more atoms in the 2D view. Extensive
experiments on the PCQM4Mv2 dataset [9] demonstrate that our proposed method
consistently outperforms other state-of-the-art methods in most situations.

##  II Related Work

###  II-A Molecular Property Prediction

According to different dimensions, molecular expressions contain 1D SMILES
code, 2D graph, and 3D graph. Due to the natural capability for representing
the molecular structure, graph neural networks (GNN) have attracted attention
in molecular property prediction, especially with the input of 2D or 3D
expressions.

####  II-A1 2D Methods

Many works focus on applying GNN to 2D graphs. For example, Gilmer et al.[2]
proposed Message Passing Neural Networks, which combine message passing and
aggregating algorithms with GNN, forming a successful framework for MPP. Liu
et al.[3] introduced a simple unsupervised representation for molecules that
embeds the vertices in the molecule graph and constructs a compact feature by
assembling the vertex embeddings. Ying et al.[10] proposed Graphormer, which
utilizes a Transformer [11] in the GNN by effectively encoding the graph
structure information. Although these methods have shown promising
capabilities, their potential is limited as they do not use 3D coordinates
that better represent molecular structure and energy.

####  II-A2 3D Methods

3D molecular structures contain spatial information important to molecular
property prediction, such as bond angles and lengths. Many recent methods try
to answer how to utilize this 3D information in the GNN model. For instance,
Klicpera et al.[5] proposed to let message embeddings interact based on the
distance between atoms and the angle between directions to improve quantum
mechanical property prediction. Lu et al.[12] proposed a Multilevel Graph
Convolutional neural Network that extracts features from the conformation and
spatial information with multilevel interactions. Liu et al.[7] proposed
SphereNet, a 3D graph network framework with spherical message passing.
3D-based methods achieve significant performance improvements. But obtaining
3D structures is time-costing, which reduces the value of 3D-based methods in
large-scale applications.

###  II-B Cross Modality Training Methods

Since training is relatively insensitive to time-consuming, an intuitive
approach is to use 2D and 3D information in training and only use 2D in
prediction.

####  II-B1 Contrastive Methods

These methods use 3D conformers and 2D structures to construct positive and
negative pairs while applying a contrastive loss during training. For example,
Liu et al.[13] proposed the Graph Multi-View Pretraining framework, in which
they perform self-supervised learning by leveraging the correspondence and
consistency between 2D topological structures and 3D geometric views. Li et
al.[14] developed GeomGCL, a dual-view geometric message passing network that
utilizes the molecular geometry across 2D and 3D views. Hu et al.[15] proposed
a strategy of pre-training GNN at the level of individual nodes as well as
entire graph, so the GNN can learn useful local and global information. Stark
et al.[16] proposed 3D Infomax, which pre-trains a GNN by maximizing the
mutual information between its 2D and 3D embeddings. The drawbacks of
contrastive methods lie in the requirements of multiple conformers. Large-
scale conformer generation is costly, while one SMILES code usually generates
only one 3D conformer with a specific force field.

####  II-B2 Knowledge Distilling Methods

Knowledge distilling (KD) is an intuitive solution for 3D to 2D cross-modality
training. As our work relies on Transformer-based methods such as Graphormer,
we only discuss Transformer related KD methods. Many works [17, 18, 19, 20,
21] have attempted to apply KD to the Transformer-based model. Upon molecular
property prediction, Zhu et al.[22] proposed ST-KD, an end-to-end Transformer
KD framework, bridging the knowledge transfer between graph-based and SMILES-
based models. But ST-KD only KD the atom to atom attentions, neglecting the
informative embeddings in atom tokens. Based on our experiment, gradient
magnitudes of atom-token distillation are discrepant and unstable due to the
variable molecular size. Distilling atom-token embeddings without correction
will cause degradation problems.

To address this issue, we propose a theoretical insight to justify how to
coordinate atom token embeddings and molecular features and formulate an auto-
correcting distilling weight to dynamically tune and balance the gradients
according to the number of atoms.

##  III Method

Figure 1: The illustration of of our CCMD framework with global molecular and
local atom information for molecular property prediction. $\mathcal{L}_{m}$
represents the global molecular distillation, $\mathcal{L}_{a}$ represents the
local atom information distillation, and coordinating denotes the coordinating
weight scaling gradient magnitudes.

###  III-A Overview

In this section, we present the key components of CCMD. As illustrated in Fig.
1, the whole graph transformer framework consists of a 3D model and a 2D
model, and the backbone is based on the graphormer [10]. Firstly, the 3D model
is trained with 3D information as a teacher model. Then, frozing the
parameters of the teacher, we train the 2D student model with a distilling
loss and a supervised target loss. The distilling module consists of the
global molecular distillation and local atom distillation, which can further
improve the performance of 2D molecular graphs. At last, for coordinating
global and local information, we formulate a novel coordinating weight to
scale and balance the gradients optimally.

###  III-B Preliminary

For convenience, let $\mathcal{D}=\\{(X,Y,E)\\}$ denotes a 2D molecular graph,
the $X$, $Y$ and $E$ denote the input of atoms, property labels and edges
respectively. In addition, we use $\widetilde{}$ mark the 3D view, so 3D
molecular graph can be denoted as
$\widetilde{\mathcal{D}}=\\{(\widetilde{X},\widetilde{Y},\widetilde{E})\\}$,
where $\widetilde{E}=f(\widetilde{P})$, and $\widetilde{P}$ represents the 3D
coordinates of atoms. For instance, the $i$-th atom embedding in $l$-th layer
can be denoted as $X_{i}^{l}$, where $i\in[0,N+1]$. $i=0$ denotes the virtual
token representing the molecular representation, $N$ denotes the total number
of atoms. $l\in[0,L]$ denotes the order number of the layer.

###  III-C Input Representations

####  III-C1 Absolute Position Encoding

Following the previous works (Park et al. 2022; Ying et al. 2021), we
integrate an absolute position encoding (APE) into graphormer [10] illustrated
in Fig.2. The APE is the adjacent edge information of the node, which could
distinguish some atoms with the same atom id easily in the input layer. Let
$X^{0}$ denote the input tokens which are integrates with the APE, $X^{0}$ is
calculated as:

$X^{0}=X\oplus E,$ (1)

where $\oplus$ is used to gather the information from neighbors. In specific,
the $i$-th input token is calculated as:

$X_{i}^{0}=MLP(\sum_{i=1}^{N+1}e_{i,j}).$ (2)

For the 2D view, ${e}_{i,j}$ is the feature of original Chemical bonds. For
the 3D view, we utilize radial basis functions (RBF) [23] as the transformed
function $\widetilde{e}_{i,j}=RBF(d_{ij})$ to obtain the edge information,
where ${d}_{ij}$ is the distance between atom i and atom j. This 3D embedding
can encode diverse geometric factors

###  III-D Backbone

The backbone used is graphormer [10]. All tokens in the $l+1$ layer is updated
by layer normalization (LN) and multi-head attention (MHA) on the tokens in
the $l$ layer

$X^{{}^{\prime}(l+1)}=MHA(LN(X^{l}))+X^{l},$ (3)

then the feed-forward blocks (FFN) is applied to the tokens in the $l+1$
layer:

$X^{l+1}=FFN(LN(X^{l+1}))+X^{{}^{\prime}(l+1)}.$ (4)

Moreover, the final molecular representation $R$ is obtained from the virtual
token of the last layer by an MLP:

$R=MLP(X_{0}^{L}).$ (5)

At last, for the supervised task, we use $L1$ loss. So the supervised term
$\mathcal{L}_{3D}$ is defined as:

$\mathcal{L}_{3D}(\widetilde{R},\widetilde{Y})=L_{1}Loss(\widetilde{R},\widetilde{Y}),$
(6)

The supervised term $\mathcal{L}_{2D}$ is defined as:

$\mathcal{L}_{2D}(R,Y)=L_{1}Loss({R},{Y}),$ (7)

###  III-E Distillation

####  III-E1 Challenges

Usually, for most of the existing distillation methods, the total loss is
defined as:

$\mathcal{L}(R,Y,\widetilde{R})={L}_{2D}(R,Y)+\mathcal{W}*L_{1}Loss({X},\widetilde{X}),$
(8)

${L}_{2D}(R,Y)$ is the task-related loss, and $L_{1}Loss({X},\widetilde{X})$
is the distilling loss on all tokens and ${W}$ is the hyper-parameter.
However, the virtual token ${X}_{0}$ is the target, which represents the
global representation of the molecular, and the other tokens ${X}_{i}(i>0)$
are the atom representations that are auxiliary to the molecular
representation. So the strategy of distillation should be divided into two
modules: global molecular distillation and local atom distillation. In
addition, We should select a proper weight to balance the global molecular
distillation and local atom distillation.

####  III-E2 Global molecular distilling

The virtual token ${X}_{0}$ represents the molecular embedding and usually be
used to predict molecular property, so distilling the virtual token embedding
provides a global view of learning on 3D geometry structural knowledge.
Considering that 2D graphs are seriously distinct from 3D conformers, we
conduct the global molecular distilling from 3D to 2D on all layers which
could learn more geometric information at different levels step by step. At
last, the target is conducted by minimizing the following loss function given
any pairs of $\widetilde{X}_{l,0}$ and $X_{l,0}$ in the training batch:

$\mathcal{L}_{m}=\sum_{l=1}^{L}L_{1}Loss(\widetilde{X}_{0}^{l},X_{0}^{l})$ (9)

####  III-E3 Local Atom Distilling

The virtual token ${X}_{0}$ is the global view of molecular, and the token
${X}_{i}(i>0)$ is the local view of molecular. The virtual token is calculated
from atom tokens by self-attention with the updating process. Considering the
relationship between the virtual token and atom tokens, we use the local atom
distilling to further boost the performance of the virtual token. One possible
way is to calculate the similarity (distance) of each pairwise atom embedding,
then the local atom calibration distillation loss is defined as:

$\mathcal{L}_{a}=\sum_{l=1}^{L}\sum_{j=1}^{N+1}L_{1}Loss(\widetilde{X}_{j}^{l},X_{j}^{l})$
(10)

Figure 2: The illustration of the backbone with absolute position encoding and
relative position encoding.

####  III-E4 Coordinating weight

For boosting the performance of distilling with global and local information,
we formulate a novel coordinating weight according to the number of molecular
size to scale gradient magnitudes. The analyzing and reducing process from two
perspectives is as follows.

Coordinating the total loss Comparing the global molecular loss
$\mathcal{L}_{m}$ and local atom loss $\mathcal{L}_{a}$, we can draw the
conclusion that $\mathcal{L}_{a}$ is unstable due to the sum operation on all
tokens. This means that when the molecular is very large(with big N), the
final loss ${L}$ is exploding. So, we should add the function $f(N)$ according
to the variable number of atoms to scale and calibrate the local atom loss
$\mathcal{L}_{a}$. The formulation is as follows:

$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{L}_{2D}+(\sum_{l=1}^{L}L_{1}Loss(\widetilde{X}_{0}^{l},X_{0}^{l})$
(11)
$\displaystyle+f(N)*\sum_{l=1}^{L}\sum_{j=1}^{N+1}L_{1}Loss(\widetilde{X}_{j}^{l},X_{j}^{l}),$

$L$ is fixed and $N$ is variable which represents the number of atoms. In
order to eliminate the effect of sum operation on all tokens, it’s intuitive
to set $f(N)=\frac{1}{N}$. In addition, we couldn’t add $\frac{1}{N}$ to the
molecular loss, which will dissipate the loss of the virtual token as N
increases, so it’s necessary to separate global and local loss, and only scale
the atom loss ${{L}_{a}}$. Generally, to make the total loss not explode, we
should add a function $f(N)\propto{O}(\frac{1}{N})$ or
$f(N)\propto{o}(\frac{1}{N})$. ${O}$ means “is of the same order as” and ${o}$
means “is ultimately smaller than”.

Coordinating the molecular gradient Since the virtual token is the target
embedding, its gradient is much important. So we will calculate the
relationship between its gradient and N. To simplify the forward of
transformer, $X_{i}^{l+1}$ in $(l+1)$ layer is formed as:

$X_{i}^{l+1}=MHA(X^{l})\\\ \&=\sum_{j=0}^{N+1}({W}^{l}_{ji}*{X}^{l}_{j})$ (12)

The ${W}^{l}_{ji}$ is the attention weight. In the back-propagation, the
gradient of the virtual token in $l$th layer is as formulated:

$\displaystyle{grad(X_{0}^{l})}$ $\displaystyle=\nabla_{X_{0}^{l}}{L}$ (13)
$\displaystyle=\sum_{k=0}^{N+1}((\nabla_{X_{k}^{l+1}}{L})*(\nabla_{X_{0}^{l}}{X_{k}^{l+1}}))$

Finally, for transformer, the virtual token gradient $grad(X_{0}^{l+1})$ in
${l}$ layer can be described in follow:

$\displaystyle{grad(X_{0}^{l})}$ $\displaystyle\propto{{O}({N}^{2})},$ (14)

Different from multi-head attention in transformer, a simple GNN always use
${MLP}$ to calculate ${X}^{l+1}_{i}$ with ${X}^{l}_{0}$ . So
$(\nabla_{X_{0}^{l}}{X_{k}^{l+1}})$ is constant, then the result is:

$\displaystyle{grad(X_{0}^{l})}$ $\displaystyle\propto{{O}({N})},$ (15)

For detailed proof, please refer to the Section1 in Supplementary Material.
For transformer, the formula indicates that the gradient of virtual token
$grad(X_{0}^{l})$ will be scaled by ${N}^{2}$, and N is the number of atoms.
So scaled with $\frac{1}{N}$ couldn’t solve the unstable gradient of the
virtual token. However, for a simple GNN, $\frac{1}{N}$ is enough.
Conclusion To this end, for transformer, we proved the weight
$f(N)=\frac{1}{{N}^{2}}$ is suitable to scale the losses and balance the
gradient magnitudes of the global molecular information and local atom
information. So the total objective function is formulated as:

$\mathcal{L}=\mathcal{L}_{2D}+(\mathcal{L}_{m}+\frac{1}{{N}^{2}}*\mathcal{L}_{a}).$
(16)

##  IV Experiments

To evaluate our method and investigate the effectiveness of the proposed
components, we carry out extensive experiments on two popular molecular
property prediction tasks, especially in PCQM4Mv2 [9]. In this section, we
introduce the detailed configuration in our experiments and present the
detailed experimental results.

###  IV-A Datasets

####  IV-A1 PCQM4Mv2 dataset

PCQM4Mv2 dataset [9] served as an ”ImageNet Large Scale Visual Recognition
Challenge” in the field of graph ML is a quantum chemistry dataset originally
curated under the PubChemQC project. It contains 3.8 million molecular graphs.
The 3D structure for training molecules is calculated by DFT, which is
equipped with only one conformer for each SMILES code.

####  IV-A2 MolHIV dataset

MolHIV dataset [24] is adopted from the MOLECULENET [25] with 41K molecular
graphs for molecular property prediction. All the molecules are pre-processed
using RDKIT [26]. Each graph represents a molecule, where nodes are atoms, and
edges are chemical bonds. Input node features are 9-dimensional, containing
atomic number and chirality, as well as other additional atom features such as
formal charge and whether the atom is in the ring. Input edge features are
3-dimensional, containing bond type, bond stereochemistry as well as an
additional bond feature indicating whether the bond is conjugated.

###  IV-B Experimental Setup

Our experiments are conducted on eight NVIDIA Tesla V100 GPUs. The
implementation is based on PyTorch [27]. We use the Graphormer architecture
[10] as the based model (L = 12, d = 768) for both 2D and 3D branches. Our
CCMD method contains improved baseline with APE, global molecular
distillation, local atom distillation, searched weight and coordinating
weight, denoted as baseline, $\mathcal{L}_{m}$, $\mathcal{L}_{a}$,
${W}_{search}$ and ${W}_{coordinating}$ respectively. In order to prevent the
explosion of a total loss, the local atom distillation loss is calculated by
means denoted as $\overline{\mathcal{L}_{a}}$, where
$\overline{\mathcal{L}_{a}}=\frac{1}{N}*\mathcal{L}_{a}$. So we set the
coordinating weight as $\frac{1}{N}$ for experiments on graphormer and $1$ for
experiments on GIN [28].

For the settings and metrics on the PCQM4Mv2 dataset, the previous methods use
the 2D view for training and validation. We add the 3D structure calculated by
DFT for training and distillation, and also evaluate on the validation with 2D
view (the test labels are no longer available and results are given over the
validation set). And split of dataset is followed by previous work [10]. We
set the batch size as 512 with 70 epochs and use the Adam algorithm [29] with
an initialized learning rate of 0.0002 and a momentum of 0.9. The results are
evaluated in terms of Mean Absolute Error (MAE), and a lower MAE indicates
better performance.

For the settings and metrics on the MolHIV dataset, we add the 3D structure
calculated by RDKIT for 3D view training and distillation. And we evaluate on
the test dataset with 2D view. And split of dataset is followed by previous
work [10]. We pretrain 3D and 2D branches on the PCQM4Mv1 dataset following
the previous works [10]. When finetune on the MolHIV dataset, we set the batch
size as 64 with 10 epochs and use Adam algorithm [29] with an initialized
learning rate of 0.0002 and a momentum of 0.9. The results are evaluated in
terms of Area Under the ROC Curve (AUC), and a higher AUC indicates better
performance.

###  IV-C Experimental Results

####  IV-C1 The comparison on PCQM4Mv2 dataset

Table I presents the results for molecular property prediction task on the
OGB-LSC PCQM4Mv2 datasets [9]. What’s worth mentioning is that contrastive
methods are hard applicable because the PCQM4Mv2 is equipped with only one
conformer for each SMILES code. Compared with established methods, our model
is currently the best on the PCQM4Mv2 leaderboard. As shown in Table I, our
approach achieves the significant MAE improvement by 6% compared with the
state-of-the-art method, indicating the effectiveness of our proposed method.

Figure 3: The learning curves of 2D supervised performance and distilling performance during training process. TABLE I: The comparison with state-of-the-art methods on the validation set of PCQM4Mv2 dataset. Our approach outperforms the best model on average by 6%. Method |  MAE
---|---
GCN [30] |  0.1379
GCN-VN[2, 30] |  0.1153
Graphormer [10] |  0.0864
EGT [31] |  0.0869
GRPE-Standard [32] |  0.0890
GPS [33] |  0.0858
Ours |  0.0809

####  IV-C2 The comparison on MOLHIV dataset

We further investigate our proposed method on the MOLHIV dataset. The
comparison results are shown in Table II. From Table II, we can see that our
method outperforms other methods, achieving prominent improvement by 1.39% on
AUC compared with baseline. In addition, the origin MOLHIV dataset is not
attached with 3D information, so we use the rdkit to generate the 3D position
as the 3D coordinates of atoms. It should be noted that the generated 3D
coordinates with randomly selected 3D conformers are not matched the real
coordinates, so the generated force fields are different from the real force
fields. The weakness of the 3D teacher model limits the performance of the
student. Our approach can still obtain improvement with biased conformers,
outperforming state-of-the-art methods and software (e.g.rdkit) generated 3D
coordinates-based results. These phenomenon demonstrate the effectiveness of
our distilling framework.

TABLE II: The results of our method compared with other methods on the test set of MOLHIV dataset. Method | AUC
---|---
DeeperGCN-FLAG[30] | 79.42 ± 1.20
PNA [30] | 79.05 ± 1.32
DGN [28] | 79.70 ± 0.97
PHC-GNN[30] | 79.34 ± 1.16
(baseline)Graphormer [10] | 80.51 ± 0.53
GraphMVP [13] | 77.0
GeomGCL[14] | 80.6 ± 0.009
EGT [31] | 80.60 ± 0.65
GRPE-Standard [32] | 81.39 ± 0.49
GPS [33] | 78.80 ± 0.0101
Ours(baseline+$\mathcal{L}_{m}$) | 81.4±0.11
Ours(baseline+$\overline{\mathcal{L}_{a}}$) | 79.0±0.01
Ours(baseline+$\mathcal{L}_{m}$ +($\mathcal{W}_{search}*\overline{\mathcal{L}_{a}})$) | 81.47±0.03
Ours(baseline+$\mathcal{L}_{m}$ +($\mathcal{W}_{coordinating}*\overline{\mathcal{L}_{a}})$) | 81.9±0.01

###  IV-D Ablation Study

We investigated the impact of these components of our proposed method on the
validation set of PCQM4Mv2 dataset, and results are shown in Table III.

TABLE III: The ablation study of our method on the validation set of the PCQM4Mv2 dataset. L denotes the conduction on the last layer and ”all” denotes distilling with the features, which are not divided into global and local respectively. Our contributions have a significant impact on performance. Method |  MAE
---|---
Graphormer [10] |  0.0864
baseline(reproduce+APE) |  0.0845
baseline+3D |  0.040
(L)baseline+all |  0.0853
(L)baseline+$\mathcal{L}_{m}$ |  0.0842
(L)baseline+$\overline{\mathcal{L}_{a}}$ |  0.0874
(L)baseline $+$ $\mathcal{L}_{m}$ $+$ $\mathcal{W}_{coordinating}$ $*$($\overline{\mathcal{L}_{a}}$) |  0.0837
baseline+$\mathcal{L}_{m}$ |  0.0822
baseline+$\overline{\mathcal{L}_{a}}$ |  0.087
baseline$+$ $\mathcal{L}_{m}$ $+$ $\mathcal{W}_{search}*$($\overline{\mathcal{L}_{a}}$) |  0.0818
baseline$+$ $\mathcal{L}_{m}$ $+$ $\mathcal{W}_{coordinating}*$($\overline{\mathcal{L}_{a}}$) |  0.0809

##### The validation of absolute position encoding (APE)

We train our graphormer with APE for 80 epochs as the baseline in Table III.
Compared to the original graphormer training with 300 epochs, our baseline
obtains a 1.52% improvement in the second row.

##### The validation of dividing local atom and global molecular

We observed that distilling the feature encoded with virtual token and atom
token in the fourth row in Table III generates a negative gain decreased by
0.94% compared with baseline. The phenomenon indicates binding virtual token
and atom tokens damages the performance of the distillation.

##### The validation of global molecular

Compared to the baseline model, the global molecular distillation
$\mathcal{L}_{m}$ slightly increases the performance by distilling the virtual
token embedding which is encoded with the structure information, reaching
0.36% improvement in the fifth row in Table III.

##### The validation of layers

In order to expand the effect of distilling molecular representation, we
conduct the distillation on all layers, which generates 2.72% improvement
compared with conducting on the last layer. This result demonstrates that
distilling on all layers superiors to distilling on the last layer.

##### The validation of local atoms

In addition, we employ local atom loss $\mathcal{L}_{a}$ on all layers. We
have observed one interesting phenomenon: when only leveraging the local atom
information distillation, the network decreased by 3.43% in sixth row. This
phenomenon demonstrates that distilling local atom information without
correction will cause degradation problems.

##### The validation of search weight

In the seventh row, the network further increases 3.2% compared with baseline
by distilling local atom and global molecular representation with coordinating
weight searched manually denoted as $\mathcal{W}_{search}$. The result
indicates that local atoms with coordinating weight can boost the performance
of global molecular. But the optimal coordinating weight is hard to search
manually due to the variable molecular size.

##### The validation of coordinating weight

Finally, when we incorporate the global molecular and local atoms with
coordinating weight $\mathcal{W}_{coordinating}$ to balance the grads
according to the number of atoms, our method in the eleventh row achieves the
best result with 0.809 MAE which achieves large improvement by 4.26% compared
with baseline. These results validate that our proposed coordination
distillation with global molecular and local atom information is effective for
molecular property prediction.

###  IV-E Visionlization of effectiveness

To understand and analyze the availability of our approach, we show the
variation of attention weight in Fig. 4. We observed that the correlation
between points is more prominent and notable compared with the baseline model.
Compared with the baseline that all atoms only focus on the 2-hop nearest
neighbors, both the global molecular token and local atom tokens could focus
on almost all the atoms with CCMD. Thus, CCMD enriched the student model with
more structured information. This phenomenon promulgates the significance of
distillation. To reveal the significance of coordinating weight, we
certificate the validation of coordinating distilling strategy in Fig. 5,
which shows that coordinating weight is superior to any manually searched
weight. To analyze the behavior of the model after adding the distillation
loss term during training, we present the learning curves of 2D supervised
performance and distilling performance in Fig. 3. We can observe that the
convergence of MAE is more stable compared with only 2D supervised loss.

Figure 4: The visionlization of student model according to the attention
weight. We compare the attention of atom tokens in yellow dashed fields and
virtual token in the blue dashed boxes.

Figure 5: The comparison between coordinating weight and manually searched
weights. The purple points are the results with manually searched weights, and
the green point is the result with coordinating weight according to the number
of atoms.

###  IV-F Effectiveness with Other Backbone

To evaluate the effectiveness of our approach, we change the gragh transformer
backbone of 2D branch into GIN [28] architecture (L = 3, d = 768) on the
PCQM4Mv2 dataset in Table IV. The 3D model based on gragh transformer is a
heavy model which has 12 blocks and 57.65 MB parameters, and the 2D model
based on GIN is a light GNN model which has 12 blocks and 1.923 MB parameters.
GIN and gragh transformer are very different in structure, making 2D-3D
distilling much more challenging, while our method in the fourth row still
increases MAE by 3.36% compared with baseline in the first row. What’s more,
the distillation with coordinating weight as $\frac{1}{N}$ in the fourth row
increases 0.67% compared with coordinating weight as $\frac{1}{{N}^{2}}$ in
the third row. This phenomenon indicates that coordinating weight with
$\frac{1}{N}$ is enough to scale the gradients for the GNN network. These
results in Table IV confirm the scalability and effectiveness of the
distilling strategy.

TABLE IV: The results with GIN backbone on the validation set of PCQM4Mv2 dataset.  Method |  MAE
---|---
baseline |  0.122
baseline+$\mathcal{L}_{m}$ |  0.1198
baseline+$\mathcal{L}_{m}$\+ $\frac{1}{N^{2}}*(\mathcal{L}_{a})$ |  0.1187
baseline+$\mathcal{L}_{m}$ \+ $\frac{1}{N}*\mathcal{L}_{a}$ |  0.1179

##  V Conclusion

In this work, we exclusively propose a novel CCMD framework for molecular
property prediction. In particular, we integrate global molecular token
distillation and local atom information distillation with a novel coordinating
weight to balance and scale the gradient magnitudes. In addition, we establish
an absolute position encoding to integrate adjacent edge information.
Experimental results on two popular molecular property prediction tasks
demonstrate the superiority of our model. In the future, we will further
experiment on the datasets from other areas to demonstrate the generalization
ability of our method and explore more superior distilling strategy.

## References

  * [1] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.  The graph neural network model.  IEEE transactions on neural networks, 20(1):61–80, 2008.
  * [2] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In International conference on machine learning, pages 1263–1272. PMLR, 2017.
  * [3] Shengchao Liu, Mehmet F Demirel, and Yingyu Liang.  N-gram graph: Simple unsupervised representation for graphs, with applications to molecules.  Advances in neural information processing systems, 32, 2019.
  * [4] Minghao Xu, Hang Wang, Bingbing Ni, Hongyu Guo, and Jian Tang.  Self-supervised graph-level representation learning with local and global structure.  In International Conference on Machine Learning, pages 11548–11558. PMLR, 2021.
  * [5] Johannes Klicpera, Janek Groß, and Stephan Günnemann.  Directional message passing for molecular graphs.  arXiv preprint arXiv:2003.03123, 2020.
  * [6] Zhuoran Qiao, Matthew Welborn, Animashree Anandkumar, Frederick R Manby, and Thomas F Miller III.  c.  The Journal of chemical physics, 153(12):124111, 2020.
  * [7] Yi Liu, Limei Wang, Meng Liu, Xuan Zhang, Bora Oztekin, and Shuiwang Ji.  Spherical message passing for 3d graph networks.  arXiv preprint arXiv:2102.05013, 2021.
  * [8] Andrey A Toropov, Alla P Toropova, Dilya V Mukhamedzhanoval, and Ivan Gutman.  Simplified molecular input line entry system (smiles) as an alternative for constructing quantitative structure-property relationships (qspr).  2005\.
  * [9] Weihua Hu, Matthias Fey, Hongyu Ren, Maho Nakata, Yuxiao Dong, and Jure Leskovec.  Ogb-lsc: A large-scale challenge for machine learning on graphs.  arXiv preprint arXiv:2103.09430, 2021.
  * [10] Chengxuan Ying, Tianle Cai, Shengjie Luo, Shuxin Zheng, Guolin Ke, Di He, Yanming Shen, and Tie-Yan Liu.  Do transformers really perform badly for graph representation?  Advances in Neural Information Processing Systems, 34:28877–28888, 2021.
  * [11] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [12] Chengqiang Lu, Qi Liu, Chao Wang, Zhenya Huang, Peize Lin, and Lixin He.  Molecular property prediction: A multilevel quantum interactions modeling perspective.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 1052–1060, 2019.
  * [13] Shengchao Liu, Hanchen Wang, Weiyang Liu, Joan Lasenby, Hongyu Guo, and Jian Tang.  Pre-training molecular graph representation with 3d geometry.  arXiv preprint arXiv:2110.07728, 2021.
  * [14] Shuangli Li, Jingbo Zhou, Tong Xu, Dejing Dou, and Hui Xiong.  Geomgcl: Geometric graph contrastive learning for molecular property prediction.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 4541–4549, 2022.
  * [15] Weihua Hu, Bowen Liu, Joseph Gomes, Marinka Zitnik, Percy Liang, Vijay Pande, and Jure Leskovec.  Strategies for pre-training graph neural networks.  arXiv preprint arXiv:1905.12265, 2019.
  * [16] Hannes Stärk, Dominique Beaini, Gabriele Corso, Prudencio Tossou, Christian Dallago, Stephan Günnemann, and Pietro Liò.  3d infomax improves gnns for molecular property prediction.  In International Conference on Machine Learning, pages 20479–20502. PMLR, 2022.
  * [17] Xiaoqi Jiao, Yichun Yin, Lifeng Shang, Xin Jiang, Xiao Chen, Linlin Li, Fang Wang, and Qun Liu.  Tinybert: Distilling bert for natural language understanding.  arXiv preprint arXiv:1909.10351, 2019.
  * [18] Siqi Sun, Yu Cheng, Zhe Gan, and Jingjing Liu.  Patient knowledge distillation for bert model compression.  arXiv preprint arXiv:1908.09355, 2019.
  * [19] Wenhui Wang, Furu Wei, Li Dong, Hangbo Bao, Nan Yang, and Ming Zhou.  Minilm: Deep self-attention distillation for task-agnostic compression of pre-trained transformers.  Advances in Neural Information Processing Systems, 33:5776–5788, 2020.
  * [20] Yonglong Tian, Dilip Krishnan, and Phillip Isola.  Contrastive representation distillation.  arXiv preprint arXiv:1910.10699, 2019.
  * [21] Jin Dong, Marc-Antoine Rondeau, and William L Hamilton.  Distilling structured knowledge for text-based relational reasoning.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 6782–6791, 2020.
  * [22] Wenhao Zhu, Ziyao Li, Lingsheng Cai, and Guojie Song.  Stepping back to smiles transformers for fast molecular representation inference.  arXiv preprint arXiv:2112.13305, 2021.
  * [23] Martin D Buhmann.  Radial basis functions.  Acta numerica, 9:1–38, 2000.
  * [24] Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec.  Open graph benchmark: Datasets for machine learning on graphs.  Advances in neural information processing systems, 33:22118–22133, 2020.
  * [25] Zhenqin Wu, Bharath Ramsundar, Evan N Feinberg, Joseph Gomes, Caleb Geniesse, Aneesh S Pappu, Karl Leswing, and Vijay Pande.  Moleculenet: a benchmark for molecular machine learning.  Chemical science, 9(2):513–530, 2018.
  * [26] Greg Landrum et al.  Rdkit: A software suite for cheminformatics, computational chemistry, and predictive modeling.  Greg Landrum, 2013.
  * [27] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. Devito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.  Automatic differentiation in pytorch.  2017\.
  * [28] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  arXiv preprint arXiv:1810.00826, 2018.
  * [29] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [30] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  arXiv preprint arXiv:1609.02907, 2016.
  * [31] Md Shamim Hussain, Mohammed J Zaki, and Dharmashankar Subramanian.  Global self-attention as a replacement for graph convolution.  2022\.
  * [32] Wonpyo Park, Woong-Gi Chang, Donggeon Lee, Juntae Kim, et al.  Grpe: Relative positional encoding for graph transformer.  In ICLR2022 Machine Learning for Drug Discovery, 2022.
  * [33] Ladislav Rampášek, Mikhail Galkin, Vijay Prakash Dwivedi, Anh Tuan Luu, Guy Wolf, and Dominique Beaini.  Recipe for a general, powerful, scalable graph transformer.  arXiv preprint arXiv:2205.12454, 2022.

##  VI section 1

In this part, we will prove the relation of gradient on the virtual token and
atom number ${N}$.

For transformer, to simplify the proof, $X_{i}^{l+1}$ in $l$ layer is formed
as:

$X_{i}^{l+1}=MHA(X^{l})=\sum_{j=0}^{N+1}({W}^{l}_{ji}*{X}^{l}_{j})$ (17)

The ${W}^{l}_{ji}$ is the attention weight, and ${W}^{l}_{ji}$ is formulated
as:

${W}^{l}_{ji}=(\frac{e^{{X}^{l}_{i}*{X}^{l}_{j}}}{\sum_{j=0}^{N+1}(e^{{X}^{l}_{i}*{X}^{l}_{j})}})$
(18)

So,

$\displaystyle{grad(X_{0}^{l})}$ $\displaystyle=\nabla_{X_{0}^{l}}{L}$ (19)
$\displaystyle=\sum_{k=0}^{N+1}((\nabla_{X_{k}^{l+1}}{L})*(\nabla_{X_{0}^{l}}{X_{k}^{l+1}}))$
$\displaystyle=\sum_{k=0}^{N+1}((\nabla_{X_{k}^{l+1}}{L})*(\nabla_{X_{0}^{l}}(\sum_{j=0}^{N+1}{W}_{jk}*{X}^{l}_{j}))$
$\displaystyle=\sum_{k=0}^{N+1}((\nabla_{X_{k}^{l+1}}{L})*(\nabla_{X_{0}^{l}}((\sum_{j=1}^{N+1}{W}_{jk}*{X}^{l}_{j})+{W}_{0k}*{X}^{l}_{0}))$
$\displaystyle=\sum_{k=0}^{N+1}((\nabla_{X_{k}^{l+1}}{L})*((\sum_{j=1}^{N+1}\nabla_{X_{0}^{l}}({{W}_{jk})*{X}^{l}_{j}})+{W}_{0k})$
$\displaystyle\propto(\sum_{k=0}^{N+1}(\nabla_{X_{k}^{l+1}}{L})*(\sum_{j=1}^{N+1}\nabla_{X_{0}^{l}}({{W}_{jk})*{X}^{l}_{j}})+{O}({N})$
$\displaystyle\propto{O}(N)*((\sum_{j=1}^{N+1}\nabla_{{X}_{0}^{l}}{{W}_{jk}})+1)$

Now, we calculate ${\nabla_{{X}_{0}}{{W}_{jk}}}$ of the l-th layer, for
convenience, ${a}$ and ${b}$ denote constants independent of the variable N
and ${X}_{0}$:

$\displaystyle\nabla_{{X}_{0}}{{W}_{jk}}$
$\displaystyle=\nabla_{{X}_{0}}(\frac{e^{{X}_{j}*{X}_{k}}}{\sum_{j=0}^{N+1}e^{{X}_{j}*{X}_{k}}})$
(20)
$\displaystyle\propto(\nabla_{{X}_{0}}(\frac{a}{\sum_{j=0}^{N+1}e^{{X}_{j}*{X}_{k}}}))$
$\displaystyle\propto(\nabla_{{X}_{0}}(\frac{a}{\sum_{j=1}^{N+1}e^{{X}_{j}*{X}_{k}}+e^{{X}_{0}*{X}_{k}}}))$
$\displaystyle\propto(\nabla_{{X}_{0}}(\frac{a}{b+e^{{X}_{0}*{X}_{k}}}))$
$\displaystyle\propto({O}(1))$

Finally, the virtual token gradient $grad(X_{0}^{l})$ in $l$ layer can be
described in follow:

$\displaystyle{grad(X_{0}^{l})}$
$\displaystyle\propto{{O}({N})}*\sum_{j=1}^{N+1}({O}(1))$ (21)
$\displaystyle\propto{{O}({N}^{2})},$

For a simple GNN like GIN:

$\displaystyle{X}^{l+1}_{i}$ $\displaystyle={W}^{l+1}_{0i}*{X}^{l}_{0}$ (22)

So,

$\displaystyle{grad(X_{0}^{l})}$ $\displaystyle=\nabla_{X_{0}^{l}}{L}$ (23)
$\displaystyle=\sum_{k=0}^{N+1}((\nabla_{X_{k}^{l+1}}{L})*(\nabla_{X_{0}^{l}}{X_{k}^{l+1}}))$
$\displaystyle=\sum_{k=0}^{N+1}((\nabla_{X_{k}^{l+1}}{L})*({W}^{l+1}_{0k}))$
$\displaystyle\propto{O}(N)$