# Production of meson molecules in ultra-peripheral heavy ion collisons

F.C. Sobrinho1, L.M. Abreu1,2, C.A. Bertulani3,4, F.S. Navarra1
1Instituto de Física, Universidade de São Paulo, Rua do Matão 1371 - CEP
05508-090, Cidade Universitária, São Paulo, SP, Brazil
2Instituto de Física, Universidade Federal da Bahia, Campus Ondina, Salvador,
Bahia 40170-115, Brazil
3Department of Physics and Astronomy, Texas A&M University-Commerce, Commerce,
Texas 75429, USA
4Institut für Kernphysik, Technische Universität Darmstadt, 64289 Darmstadt,
Germany

###### Abstract

In this work we present the first calculation of exotic charmonium production
in ultra-peripheral collisions, in which the exotic state is explicitly
treated as a meson molecule. Our formalism is general but we focus on the
lightest possible exotic charmonium state: a $D^{+}D^{-}$ molecular bound
state. It was proposed some time ago and it has been object of experimental
searches. Here we study the production of the open charm pair in the process
$\gamma\gamma\to D^{+}D^{-}$. Then we use a prescription to project the free
pair $|D^{+}D^{-}\rangle$ onto a bound state at the amplitude level and
compute the cross section of the process $\gamma\gamma\to B$ (where $B$ is the
bound state). Finally, we convolute this last cross section with the
equivalent photon distributions coming from the projectile and target in an
ultra-peripheral collision and find the $AA\to AAB$ cross section, which, for
$Pb-Pb$ collisions at $\sqrt{s_{NN}}=5.02$ TeV, is of the order of
$3\,\mu\mbox{b}$.

Quantum Chromodynamics, Exotic Charmonium, Ultra-Peripheral Collisions

###### pacs:

12.38.-t, 24.85.+p, 25.30.-c

##  I Introduction

One of the most important research topics in modern hadron physics is the
study of the exotic heavy quarkonium states exo ; exo2 . These new mesonic
states are not conventional $c\bar{c}$ configurations and their minimum quark
content is $c\bar{c}q\bar{q}$. This leads us to the main question in the
field: are these multiquark states compact tetraquarks or are they meson
molecules? So far there is no conclusive answer. One can try to address this
question with the help of experiment and study the observables: masses, decay
widths and production rates. How can multiquark states be produced? In $B$
decays and in $e^{+}e^{-}$, proton-proton, proton-nucleus and nucleus-nucleus
collisions. We will focus on the latter, which can be divided into central
(and semi-central) and ultra-peripheral (UPCs) BKN05 . In UPCs the nuclei do
not overlap and there are only few particles produced. In these collisions the
elementary processes which contribute to particle production are photon-
photon, photon-Pomeron and Pomeron-Pomeron fusion. The advantage of UPCs is
the low particle-production multiplicity, thus with a reduced background if
proper detection techniques are used. Such features have been explored at the
large hadron collider (LHC) at CERN and at the relativistic heavy ion collider
(RHIC) at Brookhaven.

In this work we will study exotic charmonium production in photon-photon
processes in nucleus-nucleus collisions. Coming back to the question
formulated above, the strategy to get the answer is to compute the cross
section for production (in UPCs) of a given exotic charmonium state assuming
that it is i) a tetraquark and also assuming that it is ii) a meson molecule.
We believe that the resulting cross sections are very different from each
other and hence, just looking at the production rate, one could experimentally
discriminate between the two configurations. Here we will address only the
production of molecules. The study of tetraquark production is in progress.

The production of hadron molecules has been discussed in the context of $B$
decays pi , in $e^{+}e^{-}$ collisions, in proton-proton pp1 ; pp2 , in
proton-nucleus and in central nucleus-nucleus collisions aa . In this work we
present the first study of meson molecule production in UPCs. The method
employed here is applicable to all molecular states. We start with the
lightest charm meson molecule: the $D^{+}D^{-}$ state (also called
$D\bar{D}$). It was predicted in the study of meson-meson interactions in the
charm sector in 15 , where it was found to be bound by about 20 MeV. The state
was confirmed in subsequent theoretical studies 16 ; 17 . More recently it was
also found in lattice calculations 18 . In 19 , it was shown that the peak in
the $D\bar{D}$ invariant mass, observed by the BELLE collaboration be , could
be well explained by the existence of a hidden charm scalar resonance below
the threshold 15 . An updated experimental work was performed in 20  and,
again, support for the $D\bar{D}$ state in the reaction $e^{+}e^{-}\to
D\bar{D}$ (and also in $\gamma\gamma\to D\bar{D}$) was found. Recent analyses
of these data were published in 21 ; 22 . A more refined theoretical work of
these reactions was performed in mb ; 23 ; ao , claiming again evidence for
this bound state.

In the next Section we present the formalism employed to describe $D^{+}D^{-}$
pair production; in Section III we present the prescription to create the
bound state; in Section IV we discuss the equivalent photon spectrum; in
Section V, performing a low energy approximation, we derive an analytical
formula for the cross section of bound state production. In the final section
we present numerical results and discussion.

##  II Production of free $D^{+}D^{-}$ pairs

There are two ways to produce a $D^{+}D^{-}$ from two photons. In the first,
one of the photons splits directly into the pair $\gamma\to D^{+}D^{-}$, where
one of the mesons is already on the mass shell, and the second photon brings
the other $D$ to the mass shell. This process can be described by a well known
hadronic effective Lagrangian, from which we obtain the pair production
amplitude. This amplitude is subsequently projected onto the amplitude for
bound state formation. If the properties of the bound state are known, the
only unknown in this formalism is the form factor, which must be attached to
the vertices to account for the finite size of the hadrons.

In the second way to produce the pair, one photon splits into a $c\bar{c}$
pair which, after interacting with the second photon, hadronizes into the
$D^{+}D^{-}$ pair. Then, using a coalescence prescription, we obtain a model
for the production of the bound state. The hadronization process involves
uncertainties related to its non-perturbative nature. Here we can not
automatically use fragmentation functions, which require a hard scale.
Moreover, the coalescence prescription contains some inherent arbitrariness.

While the relation between these two mechanisms (and whether they are
complementary or equivalent) remains to be explored, we choose to work with
the hadronic formalism. Along this line, we will study the process
$\gamma\gamma\to D^{+}D^{-}$ with the Lagrangian densities lag

$\mathcal{L}=(D_{\mu}\phi)^{*}(D^{\mu}\phi)-m_{D}^{2}\phi^{*}\phi-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}\,,$
(1)

and

$\mathcal{L}=-ig_{\gamma
D^{+}D^{*-}}F_{\mu\nu}\epsilon^{\mu\nu\alpha\beta}(D^{*-}_{\alpha}\overset{\leftrightarrow}{\partial_{\beta}}D^{+}+D^{-}\overset{\leftrightarrow}{\partial_{\beta}}D^{*+}_{\alpha})\,,$
(2)

where

$D_{\mu}\phi=\partial_{\mu}\phi+ieA_{\mu}\phi\,\,,\hskip
28.45274ptF_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}\,,$

and $\phi$, $D^{*}$ and $A_{\mu}$ represent the $D^{+}$ (or $D^{-}$), the
$D^{*+}$ (or $D^{*-}$) and the photon fields, respectively. The Feynman rules
can be derived from the interaction terms and they yield the Feynman diagrams
for the process $\gamma\gamma\to D^{+}D^{-}$ shown in Fig. 1. In the figure we
also show the quadrimomenta of the incoming photons $k^{\mu}=(E_{p},0,0,k)$,
${k^{\prime}}^{\mu}=(E_{k^{\prime}},0,0,k^{\prime})$ and of the outgoing
mesons $p^{\mu}=(E_{p},0,0,p)$,
${p^{\prime}}^{\mu}=(E_{p^{\prime}},0,0,p^{\prime})$. The total amplitude is
given by:

$iM=iM_{(a)}+iM_{(b)}+iM_{(c)}iM_{(d)}+iM_{(e)}\,,$ (3)

where

Figure 1: Feynman diagrams for the process $\gamma\gamma\rightarrow
D^{+}D^{-}$.

$\displaystyle iM_{(a)}$
$\displaystyle=2ie^{2}g_{\mu\nu}F(\bar{q}^{2})F(\bar{q}^{2})\varepsilon^{*\mu}(k)\varepsilon^{*\nu}(k^{\prime})\,,$
(4) $\displaystyle iM_{(b)}$
$\displaystyle=\varepsilon^{*\mu}(k)ieF(\hat{t})(-2p_{\mu}+k_{\mu})\frac{i}{(k-p)^{2}-m_{D}^{2}}ieF(\hat{t})(2p^{\prime}_{\nu}-k^{\prime}_{\nu})\varepsilon^{*\nu}(k^{\prime})\,,$
(5) $\displaystyle iM_{(c)}$
$\displaystyle=\varepsilon^{*\mu}(k^{\prime})ieF(\hat{u})(-2p_{\mu}+k^{\prime}_{\mu})\frac{i}{(k^{\prime}-p)^{2}-m_{D}^{2}}ieF(\hat{u})(2p^{\prime}_{\nu}-k_{\nu})\varepsilon^{*\nu}(k)\,,$
(6) $\displaystyle iM_{(d)}$
$\displaystyle=\varepsilon^{*}_{\mu}(k)[-2g\epsilon^{\sigma\mu\alpha\rho}k_{\sigma}(k_{\rho}-2p_{\rho})F(\hat{t})]\left[\frac{-i(g_{\alpha\beta}-\frac{(k-p)_{\alpha}(k-p)_{\beta}}{m_{D^{*}}^{2}})}{(k-p)^{2}-m_{D^{*}}^{2}}\right][2g\epsilon^{\delta\nu\beta\lambda}k^{\prime}_{\delta}(-k^{\prime}_{\lambda}+2p^{\prime}_{\lambda})F(\hat{t})]\varepsilon^{*}_{\nu}(k^{\prime})\,,$
(7) $\displaystyle iM_{(e)}$
$\displaystyle=\varepsilon^{*}_{\mu}(k^{\prime})[-2g\epsilon^{\sigma\mu\alpha\rho}k^{\prime}_{\sigma}(k^{\prime}_{\rho}-2p_{\rho})F(\hat{u})]\left[\frac{-i(g_{\alpha\beta}-\frac{(k^{\prime}-p)_{\alpha}(k^{\prime}-p)_{\beta}}{m_{D^{*}}^{2}})}{(k^{\prime}-p)^{2}-m_{D^{*}}^{2}}\right][2g\epsilon^{\delta\nu\beta\lambda}k_{\delta}(-k_{\lambda}+2p^{\prime}_{\lambda})F(\hat{u})]\varepsilon^{*}_{\nu}(k)\,,$
(8)

where $\bar{q}^{2}=[(k-p)^{2}+(k^{\prime}-p)^{2}]/2$ and $g=g_{\gamma
D^{+}D^{*-}}=-0.035$ lag . We have introduced the Mandelstam variables of the
elementary process, which are $\hat{s}=(k+k^{\prime})^{2}$,
$\hat{t}=(k-p)^{2}$ and $\hat{u}=(k^{\prime}-p)^{2}$. As usual, we have
included form factors, $F(q)$, in the vertices of the above amplitudes. We
shall follow kk  and use the monopole form factor given by

$F(q^{2})=\frac{\Lambda^{2}-m^{2}_{D^{(*)}}}{\Lambda^{2}-q^{2}}\,,$ (9)

where $q$ is the 4-momentum of the exchanged meson and $\Lambda$ is a cut-off
parameter. This choice has the advantage of yielding automatically
$F(m^{2}_{D})=1$ and $F(m^{2}_{D^{*}})=1$ when the exchanged meson is on-
shell. The above form is arbitrary but there is hope to improve this
ingredient of the calculation using QCD sum rules to calculate the form
factor, as done in ff , thereby reducing the uncertainties. Taking the square
of the amplitude Eq. (3) and the average over the photon polarizations it is
straigthforward to calculate the differential cross section:

$\frac{d\sigma}{d\Omega}=\frac{1}{64\pi^{2}}\frac{1}{E_{CM}^{2}}\frac{|{\bf
p}|}{|{\bf k}|}|M(\gamma\gamma\rightarrow D^{+}D^{-})|^{2}\,,$ (10)

In the center-of-mass reference frame we have ${\bf k}=-{\bf k^{\prime}}$ and
hence ${\bf p}=-{\bf p^{\prime}}$, $E_{CM}=E_{k}+E_{k^{\prime}}=2|{\bf k}|$
and $E_{CM}=E_{p}+E_{p^{\prime}}=2\sqrt{|{\bf p}|^{2}+m_{D}^{2}}$. It is then
easy to see that:

$\frac{|{\bf p}|}{|{\bf
k}|}=\sqrt{\frac{(E_{CM}^{2}-4m_{D}^{2})/4}{E_{CM}^{2}/4}}=\sqrt{1-\frac{4m_{D}^{2}}{E_{CM}^{2}}}\,,$
(11)

Inserting Eq. (11) into Eq. (10) and using $E^{2}_{CM}=\hat{s}$ we find:

$\sigma=\frac{1}{64\pi^{2}}\frac{1}{\hat{s}}\sqrt{1-\frac{4m_{D}^{2}}{\hat{s}}}\int|M(\gamma\gamma\rightarrow
D^{+}D^{-})|^{2}d\Omega\,.$ (12)

The angular integral can be done using the relations:

$\hat{t}=m_{D}^{2}-\frac{\hat{s}}{2}+\left(\sqrt{\hat{s}(\frac{\hat{s}}{4}-m_{D}^{2})}\right)\cos(\theta),\quad\hat{u}=m_{D}^{2}-\frac{\hat{s}}{2}-\left(\sqrt{\hat{s}(\frac{\hat{s}}{4}-m_{D}^{2})}\right)\cos(\theta),$

where $\theta$ is the angle between ${\bf k}$ and ${\bf p}$. We emphasize that
the only unknown in our calculation is the cut-off parameter $\Lambda$. In
what follows, we will determine it fitting our cross section to the LEP data
on the process $e^{+}e^{-}\to e^{+}e^{-}c\bar{c}$.

##  III Production of bound states

Now we describe the method to construct a bound state (denoted $B$) from the
$D^{+}D^{-}$ pair. As in pp1 , we impose phase space constraints on the
mesons, forcing them to be “close together”. Here we do this through the
prescription discussed in pes . The bound state $|B\rangle$ is defined as

$\frac{|B\rangle}{\sqrt{2E_{B}}}\equiv\int\frac{d^{3}q}{(2\pi)^{3}}\tilde{\psi^{*}}({\bf
q})\frac{1}{\sqrt{2E_{q}}}\frac{1}{\sqrt{2E_{-q}}}|{{\bf q},-{\bf q}\rangle},$
(13)

where $E_{B}$ is the bound state energy, ${\bf q}$ is the relative three
momentum between $D^{+}$ and $D^{-}$ in the state $B$, $E_{\pm q}$ are the
energies of $D^{+}$ and $D^{-}$ and $\tilde{\psi}({\bf q})$ is the bound state
wave function in momentum space, which has the following properties:

$\tilde{\psi}({\bf q})=\int d^{3}xe^{i{\bf q}\cdot{\bf x}}\psi({\bf x});\hskip
28.45274pt\int\frac{d^{3}q}{(2\pi)^{3}}|\tilde{\psi}({\bf q})|^{2}=1.$ (14)

From Eq. (13), we can write the following relation between the amplitudes:

$\frac{M(\gamma\gamma\rightarrow
B)}{\sqrt{2E_{B}}}=\int\frac{d^{3}q}{(2\pi)^{3}}\tilde{\psi}^{*}({\bf
q})\frac{1}{\sqrt{2E_{D^{+}}}}\frac{1}{\sqrt{2E_{D^{-}}}}M(\gamma\gamma\rightarrow
D^{+}D^{-}),$ (15)

We assume that the ${\bf p}\simeq{\bf p^{\prime}}$ and hence $E_{D^{+}}\simeq
E_{D^{-}}=E_{D}$ and also ${\bf q}={\bf p}-{\bf p^{\prime}}\simeq 0$.
Therefore the energy $E_{D}$ and the amplitude $M(\gamma\gamma\rightarrow
D^{+}D^{-})$ can be taken out of the integral. Moreover, since the binding
energy is small we have $E_{B}\simeq 2E_{D}$ and hence

$\displaystyle\frac{M(\gamma\gamma\rightarrow B)}{\sqrt{2E_{B}}}$
$\displaystyle=\frac{M(\gamma\gamma\rightarrow
D^{+}D^{-})}{E_{B}}\int\frac{d^{3}q}{(2\pi)^{3}}\tilde{\psi}^{*}({\bf q}),$
$\displaystyle M(\gamma\gamma\rightarrow B)$
$\displaystyle=\sqrt{\frac{2}{E_{B}}}M(\gamma\gamma\rightarrow
D^{+}D^{-})\int\frac{d^{3}q}{(2\pi)^{3}}\int d^{3}x\,\psi^{*}({\bf x})e^{i{\bf
q}\cdot{\bf x}},$
$\displaystyle=\sqrt{\frac{2}{E_{B}}}M(\gamma\gamma\rightarrow D^{+}D^{-})\int
d^{3}x\,\psi^{*}({\bf x})\,\delta^{(3)}({\bf x}),$
$\displaystyle=\psi^{*}(0)\sqrt{\frac{2}{E_{B}}}M(\gamma\gamma\rightarrow
D^{+}D^{-})\,.$ (16)

With the amplitude above we calculate the cross section for bound state
production:

$d\sigma=\frac{1}{H}\frac{d^{3}p_{B}}{(2\pi)^{3}}\frac{1}{2E_{B}}(2\pi)^{4}\delta^{(4)}(k+k^{\prime}-p_{B})|M(\gamma\gamma\rightarrow
B)|^{2},$ (17)

where $p_{B}$ is the momentum of the produced bound state and $H$ is the flux
factor. Now we will work in the center of mass frame of the $AA\to AAB$
collision, in which the momenta of the incoming photons may be different. In
this frame we have

$k=(\omega_{1},0,0,\omega_{1})\,,\hskip
28.45274ptk^{\prime}=(\omega_{2},0,0,-\omega_{2})\,,\hskip
28.45274ptp_{B}\equiv p+p^{\prime}=(E_{B},0,0,\omega_{1}-\omega_{2})\,,$ (18)

where $E_{B}=\sqrt{(\omega_{1}-\omega_{2})^{2}+m_{B}^{2}}$ and $\omega_{1}$
and $\omega_{2}$ are the energies of the colliding photons. The flux factor is
then given by

$H=4\sqrt{(k\cdot k^{\prime})^{2}-m_{k}^{2}m_{k^{\prime}}^{2}}=4k\cdot
k^{\prime}=4(k_{0}k^{\prime}_{0}-{\bf k}\cdot{\bf
k^{\prime}})=4(\omega_{1}\omega_{2}-\omega_{1}(-\omega_{2}))=2(4\omega_{1}\omega_{2})\,,$
(19)

Inserting this expression into Eq.(17) and integrating, the cross section
reads

$\displaystyle\sigma(\omega_{1},\omega_{2})$
$\displaystyle=\frac{2\pi}{2(4\omega_{1}\omega_{2})}\int\frac{d^{3}p_{B}}{2E_{B}}\delta(E_{CM}-E_{B})\delta^{(3)}({\bf
k}+{\bf k^{\prime}}-{\bf
p}_{B})\left[\frac{2}{E_{B}}|\psi(0)|^{2}|M(\gamma\gamma\rightarrow
D^{+}D^{-})|^{2}\right],$
$\displaystyle=\frac{\pi|\psi(0)|^{2}}{4\omega_{1}\omega_{2}E_{B}^{2}}|M(\gamma\gamma\rightarrow
D^{+}D^{-})|^{2}\frac{4\omega_{1}^{2}+m_{B}^{2}}{8\omega_{1}^{2}}\delta(\omega_{2}-\frac{m_{B}^{2}}{4\omega_{1}})\,,$
(20)

where we have used that $E^{2}_{CM}=4\omega_{1}\omega_{2}$.

To proceed with the calculation we need to know the bound state wave function
at the origin $|\psi(0)|^{2}$. Fortunately, in go  a similar bound state of
open charm mesons was studied with the Bethe-Salpeter equation and an
expression for the wave function was derived. In the first part of their paper
the authors present a formalism which is general and can be adapted to our
system. Formally, the Bethe-Salpeter equation reads $T=V+VGT$, where $T$ is
the two-body amplitude, $V$ is a matrix with elements $V_{ij}$ which are the
amplitudes of the $i\to j$ transitions and which are calculated from a given
effective Lagrangian. Finally $G$ is a loop function, which can be regularized
with a cut-off. Here we will just quote the main formulas needed to calculate
$\psi(0)$, which is given by

$\psi(0)=\frac{g}{(2\pi)^{3/2}}G\,,$ (21)

where

$G=-8\mu\pi\left(\Lambda_{0}-\gamma\arctan\left(\frac{\Lambda_{0}}{\gamma}\right)\right)\,,\hskip
28.45274pt\gamma=\sqrt{2\mu E_{b}}\,,\hskip
28.45274ptg^{2}=\frac{\gamma}{8\pi\mu^{2}(\arctan(\frac{\Lambda_{0}}{\gamma})-\frac{\gamma\Lambda_{0}}{\gamma^{2}+\Lambda_{0}^{2}})}.$
(22)

In the above expressions $\mu$ is the reduced mass ($\mu=m_{D}/2$),
$\Lambda_{0}$ is a cut-off parameter and $E_{b}$ is the binding energy. We
shall follow mb  and assume that $\Lambda_{0}=1$ GeV. From the above equations
we see that one can compute the (dynamically generated) mass of a bound state
and then determine its binding energy. Knowing $\mu$, $E_{b}$ and fixing
$\Lambda_{0}$, we can use the above formulas to calculate $\psi(0)$. In what
follows our reference value will be obtained using $m_{D}=1870$ MeV and the
mass of the bound state equal to $M_{B}=3723$ MeV, as found in mb . With these
numbers we get $E_{b}=17$ MeV and $|\psi(0)|^{2}=0.008$ GeV3. These will be
the values used to obtain all results, unless stated otherwise.

##  IV Equivalent photon approximation and the number of photons

The equivalent photon approximation is well known and it is described in
several papers epa ; epa2 . In general, when the photon source is a nucleus
one has to use form factors and the calculation becomes somewhat complicated.
Here we will follow epa  and define an UPC in momentum space. The distribution
of equivalent photons generated by a moving particle with the charge $Ze$ is
epa :

$n({\bf q})d^{3}q=\frac{Z^{2}\alpha}{\pi^{2}}\frac{({\bf
q}_{\perp})^{2}}{\omega\,q^{4}}d^{3}q=\frac{Z^{2}\alpha}{\pi^{2}\omega}\frac{({\bf
q}_{\perp})^{2}}{\left(({\bf
q}_{\perp})^{2}+(\omega/\gamma)^{2}\right)^{2}}d^{3}q\,,$ (23)

where $\alpha=e^{2}/(4\,\pi)$, $q$ is the photon 4-momentum, ${\bf q}_{\perp}$
is its transverse component, $\omega$ is the photon energy and $\gamma$ is the
Lorentz factor of the photon source ($\gamma=\sqrt{s}/2m_{p}$ and $m_{p}$ is
the proton mass). To obtain the equivalent photon spectrum, one has to
integrate this expression over the transverse momentum up to some value
$\hat{q}$. The value of $\hat{q}$ is given by $\hat{q}=\hbar c/2R$, where $R$
is the radius of the projetile. For Pb, $R\approx 7$ fm and hence
$\hat{q}\approx 0.014$ GeV. After the integration over the photon transverse
momentum the equivalent photon energy spectrum is given by:

$n(\omega)d\omega=\frac{2Z^{2}\alpha}{\pi}\ln\left(\frac{\hat{q}\gamma}{\omega}\right)\frac{d\omega}{\omega},$
(24)

Because of the approximations the above distribution is valid when the
condition $\omega\ll\hat{q}\gamma$ is fullfiled. Using Eq. (24) we can compute
the cross sections of free pair production, $\sigma_{P}$, and of bound state
production, $\sigma_{B}$. They are given by:

$\displaystyle\sigma_{P}(A\,A\rightarrow A\,A\,D^{+}D^{-})$
$\displaystyle=\int\limits_{m_{D}^{2}/\hat{q}\gamma}^{\hat{q}\gamma}d\omega_{1}\int\limits_{m_{D}^{2}/\omega_{1}}^{\hat{q}\gamma}d\omega_{2}\,\sigma_{P}(\omega_{1},\omega_{2})\,n(\omega_{1})\,n(\omega_{2}),$
(25) $\displaystyle\sigma_{B}(A\,A\rightarrow A\,A\,B)$
$\displaystyle=\int\limits_{m_{D}^{2}/\hat{q}\gamma}^{\hat{q}\gamma}d\omega_{1}\int\limits_{m_{D}^{2}/\omega_{1}}^{\hat{q}\gamma}d\omega_{2}\,\sigma_{B}(\omega_{1},\omega_{2})\,n(\omega_{1})\,n(\omega_{2}),$
(26)

where $\sigma_{P}(\omega_{1},\omega_{2})$ and
$\sigma_{B}(\omega_{1},\omega_{2})$ are given by Eqs. (12) (with
$\hat{s}=4\omega_{1}\omega_{2}$) and (20) respectively.

##  V The low energy approximation

###  V.1 Free pairs

At low photon energies and close to the $D^{+}D^{-}$ threshold, the produced
mesons are non-relativistic and we can use the approximation
$k-p\approx(0,0,0,m_{D})$ in the heavy meson propagator, i.e.:

$\frac{1}{(k-p)^{2}-m_{D}^{2}}\approx\frac{1}{0-m_{D}^{2}-m_{D}^{2}}=\frac{-1}{2m_{D}^{2}}\,.$

An analogous expression can be written for the $D^{*}$ propagator. From the
above relation we can see that in this low energy regime the amplitudes with
propagators are proportional to $\propto 1/m_{D}^{2}$ (Figs. 1b and 1c) and
$\propto 1/m_{D}^{*2}$ (Figs. 1d and 1e) and can be neglected when compared to
the amplitudes without propagators, such as the one of the contact interaction
in Fig. 1a. With this approximation the amplitude for $D^{+}D^{-}$ production
in the process $\gamma\gamma\rightarrow D^{+}D^{-}$ is given by:

$iM(\gamma\gamma\rightarrow D^{+}D^{-})\approx
2ie^{2}F^{2}(-m_{D}^{2})g_{\mu\nu}\varepsilon^{*\mu}(k)\varepsilon^{*\nu}(k^{\prime})\,.\\\
$ (27)

Taking the square and performing the average over the photon polarizations we
have:

$\displaystyle|M(\gamma\gamma\rightarrow D^{+}D^{-})|^{2}$
$\displaystyle=\frac{1}{4}\sum_{pol}2ie^{2}F^{2}(-m_{D}^{2})g_{\mu\nu}\varepsilon^{*\mu}(k)\varepsilon^{*\nu}(k^{\prime})\,\,(-2ie^{2})F^{2}(-m_{D}^{2})g_{\sigma\rho}\varepsilon^{\sigma}(k)\varepsilon^{\rho}(k^{\prime})\,,$
$\displaystyle=e^{4}F^{4}(-m_{D}^{2})g_{\mu\nu}g_{\sigma\rho}g^{\mu\sigma}g^{\nu\rho}=4e^{4}F^{4}(-m_{D}^{2})\,.$
(28)

Inserting this amplitude into Eq. (10) we find:

$\frac{d\sigma}{d\Omega}=\frac{e^{4}F^{4}(-m_{D}^{2})}{16\pi}\frac{1}{E_{CM}^{2}}\sqrt{1-\frac{4m_{D}^{2}}{E_{CM}^{2}}}.$
(29)

Performing the integral over the solid angle and using the definitions
$\alpha\equiv e^{2}/4\pi$ and $E_{CM}^{2}=4\omega_{1}\omega_{2}$, we find

$\sigma(\omega_{1},\omega_{2})=\frac{\pi\alpha^{2}F^{4}(-m_{D}^{2})}{\omega_{1}\omega_{2}}\sqrt{1-\frac{m_{D}^{2}}{\omega_{1}\omega_{2}}}\,,$
(30)

which is then substituted in Eq. (25) to give the final cross section for
$A\,A\rightarrow A\,A\,D^{+}D^{-}$.

###  V.2 Bound states

In the low energy approximation the produced bound state is non-relativistic
and then Eq. (16) reduces to:

$M(\gamma\gamma\rightarrow
B)=\psi^{*}(0)\sqrt{\frac{2}{m_{B}}}M(\gamma\gamma\rightarrow D^{+}D^{-})\,.$
(31)

Inserting Eq. (28) into the above equation and then using it in Eq. (20) we
have:

$\sigma(\omega_{1},\omega_{2})=\frac{32\pi^{3}\alpha^{2}F^{4}(-m_{D}^{2})|\psi(0)|^{2}}{m_{B}}\frac{1}{\omega_{1}\omega_{2}}\delta(4\omega_{1}\omega_{2}-E_{B}^{2}).$
(32)

Substituting the above expression into Eq. (26) and integrating we obtain the
final analytical expression:

$\sigma_{B}(A\,A\rightarrow
A\,A\,B)=\frac{256\pi|\psi(0)|^{2}Z^{4}\alpha^{4}F^{4}(-m_{D}^{2})}{3m_{B}^{5}}\left[\ln\left(\frac{s\hat{q}^{2}}{m_{p}^{2}m_{B}^{2}}\right)\right]^{3}\,.$
(33)

We emphasize that “low energy” here refers to the energy released by the
projectiles, i.e. the invariant mass of the photon pair. The nuclear
projectiles themselves may have very high energies.

Figure 2: Cross section for the process $e^{+}e^{-}\to c\bar{c}$ as a function of the energy $\sqrt{s}$ measured by the LEP Collaborations. Data are from Refs. lep . Purple stars from TASSO, pale red single star from JADE, bright orange stars from TOPAZ, AMY and VENUS, triangles from L3, squares and the green star from ALEPH, single diamond from DELPHI and circle from OPAL. The curves are calculated with Eqs. (12) and (25). |
---|---
(a) | (b)
Figure 3: Cross sections for free $D^{+}D^{-}$ pair production as a function of the energy $\sqrt{s}$. a) Dependence on $\hat{q}$ for fixed $\Lambda$. b) Dependence on $\Lambda$ for fixed $\hat{q}$. |  |
---|---|---
(a) | (b) | (c)
Figure 4: Cross sections for $D^{+}D^{-}$ bound state production as a function of the energy $\sqrt{s}$. a) Dependence on $\hat{q}$ for fixed $\Lambda$. b) Dependence on $\Lambda$ for fixed $\hat{q}$. c) Dependence on the binding energy for fixed $\Lambda$ and $\hat{q}$.  |  |
---|---|---
(a) | (b) | (c)
Figure 5: Cross sections as a function of $\sqrt{s}$. a) Comparison between
the free $D^{+}D^{-}$ pair and bound state $B$ cross sections. b) Comparison
between the complete numerical solution and the approximate analytical cross
section. c) Ratio between the cross section.

##  VI Numerical results and discussion

Having derived all the main formulas and discussed the numerical inputs, now
we present our numerical results. In Fig. 2 we show the cross sections for
free pair production and compare it to the existing experimental data from LEP
lep . In fact, the LEP data are for $e^{+}\,e^{-}\to
e^{+}\,e^{-}\,c\,\bar{c}$, i.e., the measured final states are $D^{+}D^{-}$
and $D^{0}\bar{D}^{0}$. We assume that these two final states have the same
cross section and, in order to compare with the data, we multiply our cross
section $\sigma(e^{+}\,e^{-}\to e^{+}\,e^{-}\,D^{+}\,D^{-})$ by a factor two.
In order to fit these data we will adapt expression (25) to electron-positron
collisions. The $\gamma\gamma\to D^{+}D^{-}$ cross section is same but the
photon flux from the electron (and also from the positron) and the integration
limits are different. The adaptation of Eq. (25) is performed in Appendix A.
Comparing our formula with these data, we determine the only parameter in the
calculation, which is the cut-off $\Lambda$. In the figure, the curves are
obtained substituting Eqs. (12) and (42) into (25). In the latter
$\hat{q}=m_{e}$. We did not attempt to perform a least chi square fit. Instead
we will carry on some uncertainty and work with the band $0.35<\Lambda<0.49$
GeV.

In Fig. 3 we show the cross section for $D^{+}D^{-}$ production. The black
solid lines show the result with our central parameter choice. Fig. 3a shows
the sensitivity of the result to the value of $\hat{q}$. In Fig. 3b we vary
the values of $\Lambda$ in the range defined in Fig. 2. In this sense we
propagate the uncertainty from the fit of the data to our results. Taking this
as the error in our result, the obtained cross section for the reaction
$Pb\,Pb\rightarrow Pb\,Pb\,D^{+}D^{-}$ at $\sqrt{s_{NN}}=5.02$ TeV is:

$\sigma(Pb\,Pb\rightarrow
Pb\,Pb\,D^{+}D^{-})=0.75^{+0.4}_{-0.4}\,\,\mbox{mb}.$ (34)

Assuming that the reaction $Pb\,Pb\rightarrow Pb\,Pb\,D^{0}\bar{D}^{0}$ has
the same cross section as the one given above for charged states, the total
cross section for charm production in photon-photon exclusive processes, we
have

$\sigma_{exclusive}(Pb\,Pb\rightarrow
Pb\,Pb\,c\bar{c})=1.5^{+0.4}_{-0.4}\,\,\mbox{mb}.$ (35)

In Appendix B, using crude approximations, we have arrived at the following
identity for inclusive cross sections:

$\sigma_{QED}(Pb\,Pb\rightarrow
Pb\,Pb\,c\,\bar{c}\,X)\approx\sigma_{QCD}(p\,p\rightarrow
p\,p\,c\,\bar{c}\,X)\,.$ (36)

This suggests that, for a given high energy, electromagnetic interactions in
Pb-Pb are as efficient as strong interactions in p-p for charm production.
Very recently, two independent analyses of the LHC data obtained estimates for
the inclusive charm production cross section in proton proton collisions at
$\sqrt{s}=5.02$ TeV, which are dominated by the strong interaction. In yg  the
authors find:

$\sigma_{inclusive}(p\,p\rightarrow
c\,\bar{c}\,X)=8.43^{+1.05}_{-1.16}\,\,\mbox{mb}.$ (37)

and a quite similar result was obtained in lu . According to (36), (35) should
be smaller (because it is exclusive) but of the order of magnitude of (37).
Considering the uncertainties, we believe that this is approximately true.

In Fig. 4 we present the cross section for bound state production and study
its dependence on $\hat{q}$ (Fig. 4a), on $\Lambda$ (Fig. 4b) and on the
binding energy $E_{b}$ (Fig. 4c). As expected, it is much smaller than the
cross section for open free pair production. However, it is encouraging to see
that at $\sqrt{s_{NN}}\approx 5.02$ TeV we have:

$\sigma(Pb\,Pb\rightarrow Pb\,Pb\,B)=3.0^{+0.8}_{-1.2}\,\,\mu\mbox{b}\,.$ (38)

This number should be compared with results found in br  and in fa . In those
papers, the production cross section of scalar states $X(3940)$ and $X(3915)$
in Pb Pb ultra-peripheral collisions at $\sqrt{s_{NN}}=5.02$ TeV were
calculated and the results were in the range

$5\leq\sigma(Pb\,Pb\rightarrow Pb\,Pb\,R)\leq 11\,\,\mu b\,.$ (39)

where $R$ stands for $X(3940)$ or $X(3915)$. The works br  and fa  are
relevant for us because there the states $R$ were also treated as molecules.
However there is an important difference. In br  and fa  the authors used the
Low formula, which connects the $\gamma\gamma\to R$ cross section with the
$R\to\gamma\gamma$ decay width (used as input). In very few cases this width
was measured and in some other very few cases the width was estimated with the
help of a formalism valid for dynamically generated (and hence molecular)
states. Here we propose a method to form the molecular state which is more
general and independent of the knowledge of the decay width. Another
difference is that the states $X(3940)$ and $X(3915)$ are significantly
heavier than the $D^{+}D^{-}$ molecule, whose mass is $3723$ MeV. Moreover, in
br  and fa  the equivalent photon calculation was done in the impact parameter
space. In spite of these differences the obtained cross sections are of the
same order of magnitude.

For completeness, in Fig. 5a we compare the cross sections for free pair and
bound state production and in Fig. 5b we compare the exact numerical
evaluation of $\sigma_{B}$ with the approximate analytical expression,
Eq.(33). We observe that the cross section obtained with the analytical
formula is accurate only at low energies. At higher energies it becomes larger
than the complete numerical formula. We can understand this behavior noticing
that in Eq.(33) we assumed that both $E_{B}$ and the form factor $F(q^{2})$
did not depend on $\omega_{1}$ nor on $\omega_{2}$ at low energies and
therefore resulted in a smaller denominator ($m_{B}$ where it should have been
$E_{B}$) and a constant argument of the form factor, which, as we can see from
Fig. 4b, is crucial to our numerical results. Nevertheless, the exact and the
analytical formula differ essentially only by a multiplicative factor close to
10. Dividing Eq. (33) by 10, it reproduces the exact formula within 20 %
accuracy in the relevant LHC range and can thus be useful for practical
applications. This is shown Fig. 5c.

To summarize, we have calculated, for the first time in the literature, the
cross section for the production of a heavy meson molecule in ultra-peripheral
collisions. We have combined a effective Lagrangian to compute the amplitude
of the process $\gamma\gamma\to D^{+}D^{-}$ with a prescription to project
this amplitude onto the amplitude for bound state formation. The resulting
$\gamma\gamma\to B$ cross section was then convoluted with the equivalent
photon fluxes from the projetile and target and the final cross section
$\sigma_{B}(AA\to AAB)$ was obtained. For $\sqrt{s_{NN}}=5.02$ TeV it is
$3.0^{+0.8}_{-1.2}\,\,\mu b$. This number is consistent with the results
obtained for other scalar exotic charmonium molecules in br  and fa . The
parameters of the calculation are $\Lambda$, $\hat{q}$ and $E_{b}$, which are
the hadronic form factor cut-off, the maximum momentum of an emitted photon
and the binding energy, respectively. All these parameters can be constrained
by experimental information and by calculations. Thus, we believe that in the
future it will be possible to increase the precision of our calculation.

###### Acknowledgements.

We are deeply indebted to K. Khemchandani and A. Martinez Torres for
instructive discussions. This work was partially financed by the Brazilian
funding agencies CNPq, CAPES, FAPESP, FAPERGS and INCT-FNA (process number
464898/2014-5). F.S.N. gratefully acknowledges the support from the Fundação
de Amparo à Pesquisa do Estado de São Paulo (FAPESP). C.A.B. acknowledges
support by the U.S. DOE Grant DE-FG02-08ER41533 and the Helmholtz Research
Academy Hesse for FAIR.

##  Appendix A Cross section of the process $e^{+}\,e^{-}\to
e^{+}\,e^{-}\,D^{+}\,D^{-}$

In this appendix we will adapt Eq. (25) to the $e^{+}e^{-}\rightarrow
e^{+}e^{-}c\bar{c}$ process. We start from Eq.(23) with $Z=1$:

$n({\bf q})d^{3}q=\frac{\alpha}{\pi^{2}\omega}\frac{({\bf
q}_{\perp})^{2}}{\left(({\bf
q}_{\perp})^{2}+(\omega/\gamma)^{2}\right)^{2}}d^{3}q\,.$ (40)

First we recall that $d^{3}q=dq_{x}dq_{y}dq_{z}=q_{\perp}dq_{\perp}d\theta
dq_{z}=1/2dq_{\perp}^{2}d\theta dq_{z}\to\pi dq_{\perp}^{2}dq_{z}$. Then we
make the following change of variables:

$dq_{\perp}^{2}dq_{z}=\frac{\omega}{\sqrt{\omega^{2}-q_{\perp}^{2}}}dq_{\perp}^{2}d\omega\,.$

After changing the variables we integrate Eq. (40) over $q_{\perp}^{2}$:

$n(\omega)=\frac{\alpha}{\pi}\int\limits_{0}^{\omega^{2}}\frac{q_{\perp}^{2}}{[q_{\perp}^{2}+(\omega/\gamma)^{2}]^{2}}\frac{1}{\sqrt{\omega^{2}-q_{\perp}^{2}}}dq_{\perp}^{2}\,.$
(41)

The solution of this integral is:

$n(\omega)=\frac{\alpha}{\pi}\frac{1}{\omega}\frac{\gamma}{(1+\gamma^{2})^{3/2}}\left[2\gamma^{2}\text{arcsinh}(\gamma)+\text{arcsinh}(\gamma)-\gamma\sqrt{1+\gamma^{2}}\right]\,.$
(42)

After these changes in Eq. (25), we can write the cross section for the
process $e^{+}\,e^{-}\rightarrow e^{+}\,e^{-}\,D^{+}\,D^{-}$ inserting Eq.
(42) into Eq. (25) and recalling that for electrons we use
$\gamma=\sqrt{s}/2m_{e}$ and also $\hat{q}=m_{e}$.

##  Appendix B $PbPb$ in QED versus $pp$ in QCD

In this paper we have been presenting predictions for quantities which are
poorly known experimentally. In order to know, at least, what to expect and to
have an idea of the order of magnitude of the cross sections we present below
an estimate of the cross sections of the QED process $Pb\,Pb\to
Pb\,Pb\,c\,\bar{c}$ and of the QCD process $p\,p\to p\,p\,c\,\bar{c}$. The
latter can be calculated with the simple convolution formula of the parton
model:

$\sigma(p\,p\rightarrow
p\,p\,c\,\bar{c})=\int\limits_{\frac{4m_{c}^{2}}{s}}^{1}dx_{1}f(x_{1})\int\limits_{\frac{4m_{c}^{2}}{x_{1}s}}^{1}dx_{2}f(x_{2})\sigma(x_{1},x_{2})\,,$
(43)

where $x_{1}$ and $x_{2}$ are the proton momentum fractions of the colliding
partons. In the above expression the integration limits come from the
kinematical constraint $x_{1}x_{2}s\geq 4m_{c}^{2}$. We know that this
reaction is domintated by the elementary process $g\,g\to c\,\bar{c}$. In a
rough approximation the gluon momentum distributions and the elementary
$\sigma(g\,g\to c\,\bar{c})$ cross section are given by:

$f(x_{1})=1/x_{1};\quad
f(x_{2})=1/x_{2};\quad\sigma(x_{1},x_{2})=\frac{\alpha_{s}^{2}}{x_{1}x_{2}s}\,.$
(44)

With these choices the integral above can be easily performed and yields:

$\sigma(pp\rightarrow
ppc\bar{c})=\frac{\alpha_{s}^{2}}{4m_{c}^{2}}\left[\ln\left(\frac{s}{4m_{c}^{2}}\right)-\left(1-\frac{4m_{c}^{2}}{s}\right)\right]\,.$
(45)

In the case of charm production in an UPC of Pb-Pb we have an analogous
convolution formula written in terms of the energies $\omega_{1}$ and
$\omega_{2}$ of the colliding photons. Assuming that the maximum energy
carried by one emitted photon is $\sqrt{s}/2$, the cross section is written
as:

$\sigma(Pb\,Pb\rightarrow
Pb\,Pb\,c\,\bar{c})=\int\limits_{\frac{2m_{c}^{2}}{\sqrt{s}}}^{\frac{\sqrt{s}}{2}}d\omega_{1}n(\omega_{1})\int\limits_{\frac{m_{c}^{2}}{\omega_{1}}}^{\frac{\sqrt{s}}{2}}d\omega_{2}n(\omega_{2})\sigma(\omega_{1},\omega_{2})\,.$
(46)

In the above expression the integration limits come from the kinematical
constraint $\hat{s}=(k_{1}+k_{2})^{2}=2(k_{1}\cdot
k_{2})=2(\omega_{1}\omega_{2}-(\omega_{1})(-\omega_{2}))=4\omega_{1}\omega_{2}\geq
4m_{c}^{2}$. The number of equivalent photons with energy $\omega$,
$n(\omega)$, and the photon-photon fusion cross section into an object with
invariant mass $4\omega_{1}\omega_{2}$ can be roughly approximated by

$n(\omega_{1})=\frac{Z^{2}\alpha}{\omega_{1}};\quad
n(\omega_{2})=\frac{Z^{2}\alpha}{\omega_{2}};\quad\sigma(\omega_{1},\omega_{2})=\frac{\alpha^{2}}{4\omega_{1}\omega_{2}}\,.$
(47)

We note the similarity between the above expressions and (44). As before this
integral can be easily solved and we find:

$\sigma(Pb\,Pb\rightarrow
Pb\,Pb\,c\,\bar{c})=\frac{Z^{4}\alpha^{4}}{4m_{c}^{2}}\left[\ln\left(\frac{s}{4m_{c}^{2}}\right)-\left(1-\frac{4m_{c}^{2}}{s}\right)\right]\,.$
(48)

Not surprisingly, (45) and (48) are identical except for the pre-factors.
Using $Z=82$, $\alpha=1/137$ and $\alpha_{s}^{2}=0.1$, we find

$\frac{\sigma(Pb\,Pb\rightarrow Pb\,Pb\,c\,\bar{c})}{\sigma(p\,p\rightarrow
p\,p\,c\,\bar{c})}=\frac{Z^{4}\alpha^{4}}{\alpha_{s}^{2}}\approx\frac{0.128}{0.1}=1.28\,.$
(49)

We could assume that $\alpha_{s}^{2}=0.2$, which is also a reasonable value.
Then the above ratio would have been $0.64$.

In Eq. (44) we could improve the approximation for $f(x)$. At increasingly
higher energies the more singular behavior of the gluon distribution can be
represented by $f(x)\simeq 1/x^{1+\delta}$, with $\delta\simeq 0.5$.
Analogously, in Eq. (47) we could improve the approximation for $n(\omega)$
including the $ln(\omega)$ correction. This would change both cross sections
in (49) in the same direction.

From this exercise we conclude that, for charm inclusive production at the
same nucleon-nucleon center of mass energy, we have:

$\sigma_{QED}(Pb\,Pb\rightarrow
Pb\,Pb\,c\,\bar{c})\approx\sigma_{QCD}(p\,p\rightarrow p\,p\,c\,\bar{c})\,.$
(50)

The above approximate identity is, of course, very crude but it tells us that
the two reactions have comparable cross sections.

## References

  * (1) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C.-P. Shen, C. E. Thomas, A. Vairo, and C.-Z. Yuan, Phys. Rept. 873, 1 (2020); R. M. Albuquerque, J. M. Dias, K. P. Khemchandani, A. Martinez Torres, F. S. Navarra, M. Nielsen and C. M. Zanetti, J. Phys. G 46, 093002 (2019).
  * (2) L. Meng, B. Wang, G. J. Wang and S. L. Zhu, Phys. Rept. 1019, 1 (2023).
  * (3) C. A. Bertulani, S. Klein and J. Nystrand, Ann. Rev. Nuc. Part. Sci. 55, 271 (2005).
  * (4) D. Marietti, A. Pilloni and U. Tamponi, Phys. Rev. D 106, 094040 (2022). Pioneering studies were presented in S. J. Brodsky and F. S. Navarra, Phys. Lett. B 411, 152 (1997).
  * (5) P. Artoisenet and E. Braaten, Phys. Rev. D 83, 014019 (2011); Phys. Rev. D 81, 114018 (2010); C. Bignamini, B. Grinstein, F. Piccinini, A. D. Polosa and C. Sabelli, Phys. Rev. Lett. 103, 162001 (2009).
  * (6) A. Esposito, E. G. Ferreiro, A. Pilloni, A. D. Polosa and C. A. Salgado, Eur. Phys. J. C 81, 669 (2021)
  * (7) H. Zhang, J. Liao, E. Wang, Q. Wang and H. Xing, Phys. Rev. Lett. 126, 012301 (2021); B. Wu, X. Du, M. Sibila and R. Rapp, Eur. Phys. J. A 57, 122 (2021).
  * (8) D. Gamermann, E. Oset, D. Strottman, and M. J. Vicente Vacas, Phys. Rev. D 76, 074016 (2007).
  * (9) J. Nieves and M. P. Valderrama, Phys. Rev. D 86, 056004 (2012).
  * (10) C. Hidalgo-Duque, J. Nieves, and M. P. Valderrama, Phys. Rev. D 87, 076006 (2013).
  * (11) S. Prelovsek, S. Collins, D. Mohler, M. Padmanath, and S. Piemonte, JHEP 06, 035 (2021).
  * (12) D. Gamermann and E. Oset, Eur. Phys. J. A 36, 189 (2008).
  * (13) P. Pakhlov et al. [Belle], Phys. Rev. Lett. 100, 202001 (2008).
  * (14) K. Chilikin et al. (Belle), Phys. Rev. D 95, 112003 (2017).
  * (15) E. Wang, W.-H. Liang, and E. Oset, Eur. Phys. J. A 57, 38 (2021).
  * (16) E. Wang, H.-S. Li, W.-H. Liang, and E. Oset, Phys. Rev. D 103, 054008 (2021).
  * (17) C. W. Xiao and E. Oset, Eur. Phys. J. A 49, 52 (2013).
  * (18) O. Deineka, I. Danilkin, and M. Vanderhaeghen, Phys. Lett. B 827, 136982 (2022).
  * (19) P. C. S. Brandão, J. Song, L. M. Abreu and E. Oset, Phys. Rev. D 108, 054004 (2023).
  * (20) X. H. Cao, M. L. Du and F. K. Guo, [arXiv:2401.16112 [hep-ph]].
  * (21) P. Lebiedowicz, O. Nachtmann and A. Szczurek, Phys. Rev. D 98, 014001 (2018).
  * (22) M. E. Bracco, M. Chiapparini, F. S. Navarra and M. Nielsen, Prog. Part. Nucl. Phys. 67, 1019 (2012).
  * (23) An introduction to quantum field theory, M. Peskin and M. Schroeder, Addison-Wesley (1996), p. 150.
  * (24) D. Gamermann, J. Nieves, E. Oset and E. Ruiz Arriola, Phys. Rev. D 81, 014029 (2010).
  * (25) C. A. Bertulani and G. Baur,Phys. Rept. 163, 299 (1988).
  * (26) G. Baur, K. Hencken, D. Trautmann, S. Sadovsky and Y. Kharlov, Phys. Rept. 364, 359 (2002).
  * (27) W. Da Silva [DELPHI], Nucl. Phys. B Proc. Suppl. 126, 185 (2004); A. Csilling [OPAL], AIP Conf. Proc. 571, 276 (2001). [arXiv:hep-ex/0010060 [hep-ex]].
  * (28) Y. Yang and A. Geiser, PoS EPS-HEP2023, 367 (2024), arXiv:2311.07523 [hep-ph]].
  * (29) C. Bierlich, J. Wilkinson, J. Sun, G. Manca, R. G. de Cassagnac and J. Otwinowski, arXiv:2311.11426 [hep-ph]].
  * (30) B. D. Moreira, C. A. Bertulani, V. P. Goncalves and F. S. Navarra, Phys. Rev. D 94, 094024 (2016).
  * (31) R. Fariello, D. Bhandari, C. A. Bertulani and F. S. Navarra, Phys. Rev. C 108, 044901 (2023).