# Cost-Efficient Deployment of a Reliable Multi-UAV Unmanned Aerial System

Nithin Babu, , Petar Popovski, , Constantinos B. Papadias Copyright (c) 2015
IEEE. Personal use of this material is permitted. However, permission to use
this material for any other purposes must be obtained from the IEEE by sending
a request to<EMAIL_ADDRESS>This version of the work has been
accepted for publication in VTC2022-Fall, Workshops. This work is supported by
the project PAINLESS which has received funding from the European Union’s
Horizon 2020 research and innovation programme under grant agreement No
812991. N. Babu and C. B. Papadias are with Research, Technology and
Innovation Network (RTIN), The American College of Greece, Greece (e-mail:
<EMAIL_ADDRESS>cpapadias@acg.edu).N. Babu, C. B. Papadias and P. Popovski are
with Department of Electronic Systems, Aalborg University, Denmark (e-mail:
niba@es.aau.dk,cop<EMAIL_ADDRESS>
###### Abstract

In this work, we study the trade-off between the reliability and the
investment cost of an unmanned aerial system (UAS) consisting of a set of
unmanned aerial vehicles (UAVs) carrying radio access nodes, called portable
access points (PAPs)), deployed to serve a set of ground nodes (GNs). Using
the proposed algorithm, a given geographical region is equivalently
represented as a set of circular regions, where each circle represents the
coverage region of a PAP. Then, the steady-state availability of the UAS is
analytically derived by modelling it as a continuous-time birth-death Markov
decision process (MDP). Numerical evaluations show that the investment cost to
guarantee a given steady-state availability to a set of GNs can be reduced by
considering the traffic demand and distribution of GNs.

###### Index Terms:

Portable access points, cost efficiency, Reliability, Markov decision process.

##  I Introduction

Unmanned aerial vehicles (UAVs) carrying radio access nodes, hereafter
referred to as portable access points (PAPs), deployed to provide temporary
cellular service in remote areas or to assist in emergencies has gained large
attention lately [1]. The architecture and Quality-of-Service (QoS)
requirements for such a system have been proposed in the 3GPP item [2].
Regardless of the limited onboard available energy constraint, the mobile
feature of a PAP may enhance the received signal-to-noise ratio (SNR) values
through better communication channels compared to a conventional fixed
infrastructure system. Since the received SNR values are functions of aerial
locations of the PAPs, efficient PAP deployment planning is of paramount
importance. In [1], the authors summarize the works that have considered
UAV(s) placement problems from an energy efficiency perspective, whereas, [3]
outlines the works that position UAV(s) to maximize communication-related
parameters such as coverage area and throughput. [4] proposes a general
probabilistic line-of-sight (LoS)-non-LoS air to ground channel model and
determines the optimal altitude that maximizes the coverage region. The
authors of [5] propose a graph-based algorithm to improve the throughput by
jointly optimizing the user association, UAV altitude, and transmission
direction, whereas in [6] we propose an energy-efficient 3D deployment of a
multi-PAP system. From [1] \- [6], the existing works maximize the respective
performance metrics by considering a constant traffic demand. However, in
practice, the traffic demands from the users are not constant and this varying
nature of the traffic demand can be efficiently utilized by an unmanned aerial
system (UAS) to reduce the investment cost: if the rate at which a new service
request arrives is much lower compared to the time required by a PAP to
complete the current request (low traffic intensity), then instead of
deploying an additional PAP, the mobile feature of the PAP can be used to
assign the same PAP to serve the new request after serving the current request
thereby reducing the investment cost. Hence the main challenge here is to find
the optimal number of PAPs required to ensure the availability of at least one
idle PAP when a new service request generates represented as a reliability
metric called the steady-state availability of the system: the higher the
number of PAPs, the higher the steady-state availability. In this work, we
determine the minimum investment cost, a function of the number of PAPs,
required to guarantee a given availability threshold in different traffic
intensity scenarios by exploiting the trade-off between the number of PAPs and
the system reliability which has, to the best of our knowledge, not been
considered in the literature.

Section II explains the system model and definitions. In Section III, firstly,
the given geographical region is represented as a set of minimum number of
circular cells each representing the coverage region of a PAP using the circle
packing theory; then in Section III-A a reliability analysis of the UAS is
performed to determine the optimal number of PAPs required to guarantee a
given steady state availability threshold by modeling the system as a birth-
death Markov process. All our main findings from the numerical evaluations are
discussed in Section IV.

Figure 1: System setup.

##  II System Model

A set of $N$ stationary ground nodes (GNs) that are uniformly distributed
along a circular geographical region of radius $R$ are considered to be served
by a set of PAPs. As shown in Fig. 1, the known locations of all the GNs and
the PAPs are registered to a centralised controller unit called the UAS
traffic management (UTM) unit [2]. The deployment and positioning of the PAPs
according to the traffic demand from the GNs are controlled by the UTM through
the PAP-UTM links. Moreover, the PAPs are assumed to be equipped with a
directional antenna of half-power beamwidth $2\theta$ with antenna gain in
direction $(\psi,\omega)$ given by,

$\displaystyle G_{\mathrm{p}}$ $\displaystyle=$
$\displaystyle\left\\{\begin{matrix}\dfrac{G_{\mathrm{o}}}{\theta^{2}}&-\theta\leq\psi\leq\theta,-\theta\leq\omega\leq\theta,\\\
\approx 0&\text{otherwise},\end{matrix}\right.$ (1)

where $G_{\mathrm{o}}\approx 2.2846$ [5]. All the GNs are equipped with an
omnidirectional antenna and are allocated with orthogonal channels.
Consequently, the given geographical area is divided into cells of circular
shape representing the possible coverage regions of the PAPs. A cell in which
at least one GN is active is called an active cell, whereas, a cell in which
no GN requires a PAP service is called an idle cell. When a GN from a
previously idle cell requires the service of a PAP, the request is sent to the
UTM through the corresponding GN-UTM link and the UTM sends an activation
command to an idle PAP through the PAP-UTM link. The radius (vertical
coordinate of a PAP) of the coverage region is jointly determined by the
uplink and downlink quality-of-service (QoS) constraints and the horizontal
plane coordinates of the PAPs are assumed to be those of the center of the
cells.

####  II-1 Channel model

The path loss between a GN located at a distance of $r_{i}$ from the center of
a cell and the corresponding PAP hovering at an altitude $h_{\mathrm{p}}$ is
taken as the probabilistic mean of the line-of-sight (LoS) and non-LoS (N-LoS)
path loss components [6]:

$\displaystyle\overline{L}(r_{i},h_{\mathrm{p}})$ $\displaystyle=$
$\displaystyle{P_{\mathrm{l}}\times
L_{\mathrm{l}}+\left(1-P_{\mathrm{l}}\right)\times L_{\text{nl}}},$ (2)
$\displaystyle=$
$\displaystyle\underbrace{\frac{r_{i}^{2}+h_{\mathrm{p}}^{2}}{g_{0}}}_{\text{FSPL}}\times\underbrace{[\eta^{2}_{\text{nl}}+P_{\mathrm{l}}(\eta^{2}_{\mathrm{l}}-\eta^{2}_{\text{nl}})],}_{\eta_{\mathrm{m}}(\phi_{i}),\,\,\text{mean\,\,additional\,\,
path\,\, loss}}$

where $g_{0}$ represents the channel gain at a reference distance of 1m;
$P_{\mathrm{l}}=1/\\{1+a\exp\left[-b(\phi_{i}-a))\right]\\}$ is the LoS
probability with $a$ and $b$ being the environment-dependent parameters given
in [4] and $\phi_{i}=(180/\pi)\text{tan}^{-1}(h_{\mathrm{p}}/{r_{i}})$;
$\eta_{\mathrm{l}}$ and $\eta_{\text{nl}}$ are the mean value of the
additional path loss due to shadowing for LoS and N-LoS links, respectively.

####  II-2 PAP coverage region

A GN is considered to be in the coverage region of a PAP if it has minimum
downlink and uplink SNR values. The SNR received at the $i^{\text{th}}$ GN
from a PAP (downlink) is
$\frac{G_{\mathrm{p}}P}{\overline{L}(r_{i},h_{\mathrm{p}})\sigma^{2}}$, where
$P$ and $\sigma^{2}$ are the power allocated to a GN and noise power,
respectively. From (2),
$\overline{L}(R^{\mathrm{d}}_{\mathrm{p}},h_{\mathrm{p}})\geq\overline{L}(r_{j},h_{\mathrm{p}})\,\forall
r_{j}\leq R^{\mathrm{d}}_{\mathrm{p}}$. Hence, for a given downlink QoS
constraint, $\Gamma^{\mathrm{d}}$, considering the SNR at an edge GN, the
downlink coverage radius is determined as,

$\displaystyle R_{\mathrm{p}}^{\mathrm{d}}(\Gamma^{\mathrm{d}},\theta)$
$\displaystyle=$
$\displaystyle\sqrt{\frac{G_{\mathrm{p}}g_{0}P\text{sin}^{2}\theta}{\Gamma^{\mathrm{d}}\sigma^{2}\eta_{\mathrm{m}}(\theta)}}.$
(3)

For the uplink, each GN chooses its transmit power according to the uplink
power control specified in the 3GPP technical report [7]. Then the transmit
power for the $i^{th}$ GN (in Watts) is given by,

$\displaystyle\overline{P}_{i}$ $\displaystyle=$
$\displaystyle{P_{a}\overline{L}(r_{i},h_{\mathrm{p}})},$ (4)

where $P_{\mathrm{a}}$ is the target arrived power at the PAP. Accordingly,
the received SNR from all the GNs at a PAP will be the same:
$(G_{\mathrm{p}}P_{\mathrm{a}})/\sigma^{2}$. If $\Gamma^{\mathrm{u}}$ is the
minimum uplink QoS threshold, then
$P_{\mathrm{a}}=\left(\Gamma^{\mathrm{u}}\sigma^{2}\right)/G_{\mathrm{p}}$.
Also, the average power transmitted by a GN should be less than the maximum
power $P_{\text{max}}$ that decides the uplink coverage radius
$R^{\mathrm{u}}_{\mathrm{p}}$:

$\displaystyle
P_{\mathrm{a}}\overline{L}(R_{\mathrm{p}}^{\mathrm{u}},h_{\mathrm{p}})\leq
P_{\text{max}},$ $\displaystyle
R_{\mathrm{p}}^{\mathrm{u}}(\Gamma^{\mathrm{u}},\theta)=\sqrt{\frac{G_{\mathrm{p}}g_{0}P_{\text{max}}\text{sin}^{2}\theta}{\Gamma^{\mathrm{u}}\sigma^{2}\eta_{\mathrm{m}}(\theta)}}.$
(5)

Therefore the coverage radius of a PAP is given by,
$R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)=\text{min}\left\\{R_{\mathrm{p}}^{\mathrm{d}}(\Gamma^{\mathrm{d}},\theta),R_{\mathrm{p}}^{\mathrm{u}}(\Gamma^{\mathrm{u}},\theta)\right\\}$
and the corresponding hovering height of the PAP is
$h_{\mathrm{p}}=R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)\text{tan}\theta$.

##  III Cost Vs Reliability

In this section, we analyze the trade-off between the investment cost and
reliability of the UAS. The investment cost is determined by the total number
of PAPs required to guarantee a given reliability threshold to the GNs. The
given circular geographical region needs to be covered by a set of PAP
coverage regions of radius
$R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)$. If
$R>R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)$, we need
multiple cells to cover the given region. With the determined value of
$R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)$ from Section
II, the problem of finding the minimum number of required cells and their
locations is solved using a regular-hexagon-based 7-circle multi-level circle
packing as detailed in Algorithm 1: the maximum radius of a circle which can
be covered by 7 equi-radius circles of radius
$R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)$ is
$R=2R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)$; one of
these 7 smaller circles will be concentric with the given circular region of
radius $R$, and the centers of the remaining 6 circles lie on the vertices of
a regular hexagon of side length
$\sqrt{3}R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)$ [8].
If $(x_{l-1},y_{l-1})$ is the center of the region to be covered, then centers
of the 7 smaller circles of radius $r_{l}$ that covers the given region have
the coordinates $\\{(x_{j},y_{j})\\}$ where,

$\displaystyle x_{j}=x_{l-1}+r_{l}\sqrt{3}\text{cos}\left(\dfrac{2\pi
j}{6}\right)\quad\forall j\in\left\\{0,1,...5\right\\},$ (6) $\displaystyle
y_{j}=y_{l-1}+r_{l}\sqrt{3}\text{sin}\left(\dfrac{2\pi
j}{6}\right)\quad\forall j\in\left\\{0,1,...5\right\\},$ (7)
$\displaystyle(x_{6},y_{6})=(x_{l-1},y_{l-1}).$ (8)

Algorithm 1 arranges the 7-circle packing pattern on multiple levels as shown
in Fig. 2: in the first level, the centers of 7 smaller circles of radius
$R/2$ are determined using step 4 (set of blue circles); if the radius of the
smaller circle is less than or equal to the coverage radius of a PAP, the
packing stops otherwise the packing continues to cover each smaller circle
obtained from the previous level (set of green circles). This continues until
the radius of the circle to be covered is less than or equal to
$R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)$. In step 10,
the circles representing the cells without any uncovered GNs (dotted green
circles) are discarded to obtain the minimum number of cells.

Remark 1: The number of levels required to cover a given geographical region
of radius $R$ using Algorithm 1 is,
$l_{\text{max}}=\left\lceil\text{log}_{2}\left(\dfrac{R}{R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)}\right)\right\rceil.$

It should be noted that, with given $R$ and
$\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}}$, and $\theta$, step 2 to 9 takes a
running time of $O\left(l_{\text{max}}\right)$, whereas step 10 requires
$\sum\limits_{c=0}^{7^{l_{\text{max}}}-1}(7^{l_{\text{max}}}-c)(7^{l_{\text{max}}}-c-1)$
pair checks; hence the complexity of Algorithm 1 is, at worst
$O\left[l_{\text{max}}+7^{2l_{\text{max}}}\right]$.

1 Input: $R$,
$R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta),l=0$,
$\mathcal{C}_{l}=\\{(0,0)\\}$;

2 while  _1_ do

3 $l=l+1;$ $r_{l}=R/2$;

4 For each $\mathbf{c}_{l-1}=(x_{l-1},y_{l-1})\in\mathcal{C}_{l-1}$, find the
set of center of PAP coverage circles using (6) - (8): $\mathcal{C}_{l}$;

5 if  _$r_{l}\leq
R_{\mathrm{p}}(\Gamma^{\mathrm{d}},\Gamma^{\mathrm{u}},\theta)$_ then

6 break;

7 else

8 $R=R/2$;

9 goto Step 3;

10

11

12Remove the circles with zero non-shared GNs;

Output:Total number of cells, $n=|\mathcal{C}_{l}|$.

Algorithm 1 Optimal Number of Cells

Figure 2: Position of cells for R=800m.

###  III-A Reliability Analysis and Optimal Number of PAPs.

In this section, we find the optimal number of PAPs required to guarantee a
given reliability threshold to the GNs distributed along $n$ cells. Let $A(t)$
be an event that a PAP is available for serving the request from a new active
cell at time $t$; hence, the steady-state availability is defined as,

$\displaystyle A$ $\displaystyle=$
$\displaystyle\lim_{t\to\infty}\mathcal{P}[A(t)],$ (9)

where $\mathcal{P}[A(t)]$ is the probability of $A(t)$; here, the
unavailability ($1-A$) is only due to busy PAPs, but not the wireless
environment. We consider $n$ cells determined using Algorithm 1 as $n$ sources
that generate $n$ independent Poisson inputs with common intensity $\lambda$.
Also, the service times of PAPs at cells are considered to be independent,
exponentially distributed random variables with parameter $\kappa$. The
deployment of an additional PAP is only required when a new request originates
from an idle cell. Additionally, the probability that more than one idle cell
becomes active or vice versa is assumed to be negligible. Let $u$ be the
number of PAPs available in the UAS to serve $n$ cells; thus the considered
system can be modelled as a birth-death Markov process as shown in Fig. 3 with
state-space $Z=\\{0,1,2,..u\\}$ representing the number of serving PAPs and
transition rates,

$\displaystyle\lambda_{j}=(n-j)\lambda;\quad\forall j=\\{0,1,2,...(u-1)\\},$
(10) $\displaystyle\kappa_{j}=j\kappa\quad\forall j=\\{1,2,...u\\}.$ (11)

Figure 3: State transition diagram of the UAS.

Let $\\{p_{j},j\in Z\\}$ be the steady-state probability distribution of the
considered Markov process. Then the state equations from Fig. 3 are expressed
as [9],

$\displaystyle\kappa p_{1}=n\lambda p_{0},$ (12) $\displaystyle(n-j+1)\lambda
p_{j-1}+(j+1)\kappa p_{j+1}=\left(j\kappa+(n-j)\lambda\right)p_{j},$
$\displaystyle(n-u+1)\lambda p_{u-1}=u\kappa p_{u}.$ (13)

From (III-A), $p_{1}=n\delta p_{0}$, where $\delta=\lambda/\kappa$. By
successive substitution we get,

$\displaystyle p_{j}$ $\displaystyle=$ $\displaystyle
C(n,j)\delta^{j}p_{0}\quad\forall j\in\mathcal{Z},$ (14)

where $C(n,j)=n!/(j!(n-j)!)$. For a finite $R$ value, the number of cells will
be finite; hence, the stationary state probabilities should satisfy the
normalizing condition $\Sigma_{j=0}^{u}p_{j}=1$:

$\displaystyle p_{0}+\Sigma_{j=1}^{u}C(n,j)\delta^{j}p_{0}=1,$ (15)
$\displaystyle p_{0}=\frac{1}{1+\Sigma_{j=1}^{u}C(n,j)\delta^{j}}.$ (16)

Let $F$ be an event that an active cell does not find an available PAP, then
using (14) and (16),

$\displaystyle P(F)$ $\displaystyle=$ $\displaystyle
p_{u}=\frac{C(n,u)\delta^{u}}{1+\Sigma_{j=1}^{u}C(n,j)\delta^{j}}.$ (17)

Hence, the steady-state availability of the system can be determined as,

$\displaystyle A$ $\displaystyle=$ $\displaystyle
1-P(F)=1-\frac{C(n,u)\delta^{u}}{1+\Sigma_{j=1}^{u}C(n,j)\delta^{j}}.$ (18)

Thus, for a given availability threshold, $\rho$, the optimal number of PAPs
required is obtained using (18) as,

$\displaystyle u_{\text{opt}}$ $\displaystyle=$
$\displaystyle\text{min}\\{(u|A>=\rho),n\\}.$ (19)

The degree of PAP utilization defined as the ratio of the mean number of
deployed PAPs to the total number of PAPs is given by,

$\displaystyle\eta(u_{\text{opt}})$ $\displaystyle=$
$\displaystyle\frac{\Sigma_{j=1}^{u_{\text{opt}}}jp_{j}}{u_{\text{opt}}}=\frac{\Sigma_{j=1}^{u_{\text{opt}}}jC(n,j)\delta^{j}}{u_{\text{opt}}\left(1+\Sigma_{j=1}^{u_{\text{opt}}}C(n,j)\delta^{j}\right)}.$
(20)

##  IV Numerical Evaluation and Conclusion

In this section, we present our main findings from numerical evaluations.The
considered simulation parameters are $g_{0}=1.42\times 10^{-4}$,
$\Gamma^{\mathrm{d}}=\Gamma^{\mathrm{u}}=100$, $\sigma^{2}=1.25\times
10^{-14}W$, $P=1$mW, $P_{\text{max}}=1W$. Fig. 4 shows the variations of the
average number of cells required (rounded to the nearest integer) as a
function of the number of GNs for suburban and urban deployment scenarios; the
corresponding $(a,b,\eta_{\mathrm{l}},\eta_{\mathrm{nl}},\theta)$ parameters
used are $(4.83,0.43,1.01,11.22,70^{o})$, $(9.6,0.16,1.12,10,52^{o})$ [4]. The
averaging is done over the minimum number of cells obtained using Algorithm 1
for 600 independent realizations of GNs. From the figure, for a given number
of GNs, the number of cells required increases as the radius of both the
suburban and urban region increase because of the widely distributed GNs in a
larger geographical region. Additionally, for a given QoS threshold, the
number of cells required to cover any $N$ GNs in the urban region is higher
than that in the suburban region: this is because the radius of the coverage
region of a PAP decreases as we move from suburban to urban region due to an
increase in the number and density of buildings ($\eta_{m}(\theta)$).

Figure 4: Average number of cells Vs number of GNs .

Fig. 5 is plotted considering 10 cells ($R=600$m, $N=15$, urban) and it shows
the variation of the availability and PAP utilization factor, given by (18)
and (20), as a function of the normalized investment cost; the normalized
investment cost is taken as the ratio of $u_{\text{opt}}$ to the total number
of cells. As seen in the figure, the steady-state availability increases with
deploying more PAPs (more investment cost). However, for low traffic intensity
scenarios, $\delta=0.1$, a steady-state availability of 1 is achievable with
PAPs equal to half the number of cells. This benefit of investment cost
reduction becomes less with an increase in the traffic intensity: since the
rate of arrival of the new service request is higher than the service
duration, the probability of finding an idle PAP to serve the request from a
new active cell is very small thereby demanding a deployment scheme that
assigns one dedicated PAP per cell.

Figure 5: Availability Vs cost reduction.

Fig. 6 is plotted by considering a suburban deployment scenario with 25 GNs;
it shows the variation of the average normalized cost incurred to guarantee a
given availability threshold to GNs for different traffic intensity scenarios.
As expected, the average number of cells required increases with an increase
in $R$. Both the availability threshold and the traffic intensity decide the
investment cost; the cost can be reduced either if the intensity of the
traffic demand from the GNs is low or if the availability threshold is light.
Moreover, the cost reduction is more noticeable if the GNs are widely
separated. Hence, from Fig. 5 and Fig. 6, the number of PAPs required should
be economically determined by considering not only the area of the
geographical region but also the traffic intensity and the distribution of the
GNs in the considered region.

Figure 6: Normalized cost Vs radius of the region.

For a given geographical region with a set of GNs, the minimum number of cells
required is determined using Algorithm 1, and depending on the type of service
required by the GNs and the traffic intensity, a suitable availability
threshold could be selected which in turn decides the optimal number of
required PAPs and the corresponding investment cost. The optimal number of PAP
determination by considering a general distribution of the traffic intensity
is left as future work.

## References

  * [1] S. Shakoor, Z. Kaleem, M. I. Baig, O. Chughtai, T. Q. Duong, and L. D. Nguyen, “Role of UAVs in public safety communications: Energy efficiency perspective,” _IEEE Access_ , vol. 7, pp. 140 665–140 679, 2019.
  * [2] “3GPP; Technical Specification Group Services and System Aspects; unmanned aerial system (uas) support,” _Stage 1; Release 17_ , 2017.
  * [3] C. T. Cicek, H. Gultekin, B. Tavli, and H. Yanikomeroglu, “UAV Base Station Location Optimization for Next Generation Wireless Networks: Overview and Future Research Directions,” in _2019 1st International Conference on Unmanned Vehicle Systems-Oman (UVS)_. IEEE, 2019, pp. 1–6.
  * [4] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP Altitude for Maximum Coverage,” _IEEE Wireless Communications Letters_ , vol. 3, no. 6, pp. 569–572, 2014.
  * [5] W. Huang, D. M. Kim, W. Ding, and P. Popovski, “Joint Optimization of Altitude and Transmission Direction in UAV-Based Two-Way Communication,” _IEEE Wireless Communications Letters_ , vol. 8, no. 4, pp. 984–987, 2019\.
  * [6] N. Babu, C. B. Papadias, and P. Popovski, “Energy-Efficient 3D Deployment of Aerial Access Points in a UAV Communication System,” _IEEE Communications Letters_ , vol. 24, no. 12, pp. 2883–2887, 2020.
  * [7] “3GPP, Physical Layer Procedures,” _TR 36.213, Sep. 2015, v 10.12_ , 2017\.
  * [8] G. F. Tóth, “Thinnest Covering of a Circle by Eight, Nine, or Ten Congruent Circles,” _Combinatorial and computational geometry_ , vol. 52, no. 361, p. 59, 2005.
  * [9] F. Beichelt, _Stochastic Processes in Science, Engineering and Finance_. Chapman and Hall/CRC, 2006.


# Femtoscopic correlation function for the $T_{cc}(3875)^{+}$ state

I. Vidaña1, A. Feijoo2, M. Albaladejo2, J. Nieves2 and E. Oset2,3,4 1Istituto
Nazionale di Fisica Nucleare, Sezione di Catania, Dipartimento di Fisica
“Ettore Majorana”, Università di Catania, Via Santa Sofia 64, I-95123 Catania,
Italy  2Instituto de Física Corpuscular, Centro Mixto Universidad de Valencia-
CSIC, Institutos de Investigación de Paterna, Aptdo. 22085, E-46071 Valencia,
Spain  3Departamento de Física Teórica and IFIC, Centro Mixto Universidad de
Valencia-CSIC, Institutos de Investigación de Paterna, Aptdo. 22085, E-46071
Valencia, Spain  4Department of Physics, Guangxi Normal University, Guilin
541004, China

###### Abstract

We have conducted a study of the femtoscopic correlation functions for the
$D^{0}D^{*+}$ and $D^{+}D^{*0}$ channels that build the $T_{cc}$ state. We
develop a formalism that allows us to factorize the scattering amplitudes
outside the integrals in the formulas, and the integrals involve the range of
the strong interaction explicitly. For a source of size of 1 fm, we find
values for the correlation functions of the $D^{0}D^{*+}$ and $D^{+}D^{*0}$
channels at the origin around 30 and 2.5, respectively, and we see these
observables converging to unity already for relative momenta of the order of
200 MeV. We conduct tests to see the relevance of the different contributions
to the correlation function and find that it mostly provides information on
the scattering length, since the presence of the source function in the
correlation function introduces an effective cut in the loop integrals that
makes them quite insensitive to the range of the interaction.

Keywords: Femtoscopy, $T_{cc}$

##  I Introduction

The study of molecular states has long history (see e.g. Refs. Oller00 ; Guo18
; Dong21 ; Guo21 ; Dong21b ) implying that many known mesonic states would
actually contain four quarks (two quarks and two antiquarks), and baryonic
states five quarks (four quarks and one antiquark), hence diverting from the
$q\bar{q}$ standard nature associated to the mesons and $qqq$ for baryons. The
big experimental boost to this idea came from the discovery of the $P_{c}$ and
$P_{cs}$ pentaquark states by LHCb Aaij15 ; Aaij16 ; Aaij19 ; Aaij21 ; Aaij22b
. While these states still contain quantum numbers compatible with a three
quark structure, their large mass precludes them to correspond to excited non
charmed $qqq$ baryons. The cleanest cases, however, were found with the
discovery of mesonic states containing two open quarks of different nature.
One of them was the $X_{0}(2900)$ [$T_{cs}(2900)$] Aaij20 ; Aaij20b , found in
the $D^{-}K^{+}$ mass distribution in the $B^{+}\rightarrow D^{+}D^{-}K^{+}$
decay, with $\bar{c}\bar{s}$-quark content. The other one is the
$T_{cc}(3875)^{+}$ state reported in Aaij22 ; Aaij22c , observed in the
$D^{0}D^{0}\pi^{+}$ mass distribution, which has two open charmed quarks. The
observation of this extremely narrow state, which has a small binding with
respect to the $D^{0}D^{*+}$ threshold $m_{\text{th}}=3875.09$ MeV, has
spurred lots of interest in the hadronic community. Its mass and width are
$m_{\text{th}}+\delta m_{\text{exp}}$, with $\delta m_{\text{exp}}=-360\pm
40^{+4}_{-0}$ keV, and $\Gamma=48\pm 2^{+0}_{-14}$ keV, respectively Aaij22c .

The proximity of the mass of this exotic state to the $D^{0}D^{*+}$ and
$D^{+}D^{*0}$ thresholds has given support to the molecular picture in these
channels Dong21 ; Feijoo21 ; Ling22 ; Fleming21 ; Ren22 ; Chen22 ;
Albaladejo22 ; Du22 ; Baru:2021ldu ; Santowsky22 ; Deng22 ; Ke22 ; Agaev22 ;
Kamiya22 ; Meng22 ; Abreu22 ; Chen22b ; Albaladejo:2022sux , although
tetraquark schemes have also been invoked, even prior to the $T_{cc}$
discovery Ballot83 ; Zouzou86 . Even then, the extreme proximity of the states
to the $D^{0}D^{*+}$ and $D^{+}D^{*0}$ thresholds makes it unavoidable to
include these degrees of freedom in any analysis of the experimental data as
shown in Guo21 .

Additional independent information on this state would be most welcome to
further learn about its properties and nature. Hence, the information that can
be provided by the femtoscopic correlation functions with the observation of
the $D^{0}D^{*+}$ and $D^{+}D^{*0}$ channels in heavy ion collisions, should
be most valuable. The idea of using the correlation functions to learn about
hadron interactions is relatively new, but there are already many theoretical
Morita15 ; Onishi16 ; Morita16 ; Hatsuda17 ; Mihaylov18 ; Haidenbauer19 ;
Morita20 ; Kamiya20 ; Kamiya22b  and experimental Adamczyk15 ; Acharya17 ;
Adam19 ; Acharya19 ; Acharya19b ; Acharya19c ; Acharya20 ; Acharya20b ;
Acharya20c ; Acharya21 ; Acharya21b ; Fabietti21  studies on the subject in
the strangeness sector, and work on the charm sector has already started (see
e.g., the ALICE collaboration analysis on the $Dp$ channel Acharya22 ).

In this work, we take the case of the $T_{cc}(3875)^{+}$ exotic state and
evaluate the correlation function for the channels $D^{0}D^{*+}$ and
$D^{+}D^{*0}$. Work on this issue has already been done in Kamiya22  using
coordinate space wave functions and potentials. Here we divert from this
picture and use a momentum space formalism, as has also been used recently in
Wei23  for the correlations generated from the $DK$ interaction. We rely upon
the framework used in Feijoo21 , where the local hidden gauge approach Bando88
; Harada03 ; Meissner88 ; Nagahiro09 , exchanging vector mesons as a source of
the interaction, has been used. The description of the $T_{cc}$ was done using
the two coupled channels $D^{0}D^{*+}$ and $D^{+}D^{*0}$, and a bound state
was obtained which has essentially111Isospin symmetry is not exact due to the
small difference between the $D^{0}D^{*+}$ and $D^{+}D^{*0}$ thresholds, which
cannot be neglected because the remarkable proximity of the mass of the
$T_{cc}$ to the $D^{+}D^{*0}$ threshold.  isospin $I=0$, as also determined in
the experimental work Aaij22c . We compare our results with those obtained in
Kamiya22  and also make some comments concerning the formalism used in Wei23 .
We take advantage to show that the correlation functions in this case provide
mostly information on the scattering length and very little on the range of
the interaction. Indeed, we discuss how the range of the interaction affects
the correlation functions, concluding that, due to the extension of the
source, the results are rather insensitive to the range of the interaction.

The manuscript is organized as follows. A description of the correlation
function formalism is presented in Sec. II followed by a brief summary of the
study of the $T_{cc}$ state of Ref. Feijoo21  in Sec. III. The results for the
correlation functions of the $D^{0}D^{*+}$ and $D^{+}D^{*0}$ pairs and the
different tests are discussed in Sec. IV. In Sec. V, we obtain the correlation
functions using as an input the $T$-matrix of Ref. Albaladejo22  and compare
with those obtained in Sec. IV. Finally, the main conclusions are presented in
Sec. VI.

##  II Formalism for the correlation function

The correlation function is defined as the ratio of the probability to find
the two interacting particles and the product of the probabilities to find
each individual particle Lisa05 . Within certain approximations, it can be
written as Onishi16

$C(\vec{p}\,)=\int d^{3}\vec{r}S_{12}(\vec{r}\,)|\Psi(\vec{r},\vec{p}\,)|^{2}$
(1)

where $S_{12}(\vec{r}\,)$ is the so-called source function, usually
parametrized as a Gaussian normalized to unity,

$S_{12}(\vec{r}\,)=S_{12}(r)=\frac{1}{(\sqrt{4\pi})^{3}R^{3}}\mbox{exp}\left(-\frac{r^{2}}{4R^{2}}\right)\
.$ (2)

Here $R$ is the size of the source and it takes values in the range $1-5$ fm
depending on the kind of scattering used to produce the particles. Typically,
$R\simeq 1$ fm for proton-proton collisions and $R\simeq 5$ fm in the case of
heavy ion collisions. In addition, $\vec{p}$ in Eq. (1) is the relative
momentum of the two particles, i.e., the momentum of each particle in the
center-of-mass (c.m.) frame of the two particles, while
$\Psi(\vec{r},\vec{p}\,)$ is the c.m. outgoing wave function of the two
particles.222 In the case of dealing with complex optical potentials, it is
useful to recall that in Eq. (1) one should use
$\Psi^{(-)}(\vec{r},\vec{p}\,)$ and that
$\left[\Psi^{(-)}(\vec{r},\vec{p}\,)\right]^{*}$, as a block, corresponds to
an incoming solution of the Schrödinger equation Nieves93  (see also appendix
A of Ref. Kamiya22  reaching the same conclusion). We can thus write the
Schrödinger equation for the interacting wave function $\Psi$ of the pair with
$H=H_{0}+V$ ($H_{0}$ the kinetic energy)

$H\Psi=(H_{0}+V)\Psi=E\Psi\Rightarrow(E-H_{0})\Psi=V\Psi$ (3)

while the free wave function $\Phi$, when $V\equiv 0$, for the same energy $E$
satisfies,

$(E-H_{0})\Phi=0\ .$ (4)

Subtracting Eqs. (3) and (4) we obtain

$\Psi=\Phi+\frac{1}{E-H_{0}}V\Psi\Rightarrow\Psi=\Phi+\frac{1}{E-H_{0}}T\Phi$
(5)

where in the last step of Eq. (5) we have used the definition of the
scattering matrix $T$, $T\Phi\equiv V\Psi$.

Eq. (5) gives us the full wave function $\Psi$ in terms of free one, $\Phi$,
and the scattering matrix. Using the normalization of the state with momentum
$\vec{p}$, used in Gamermann10 , we have

$\langle\vec{p}\,^{\prime}|\vec{p}\,\rangle=\delta^{3}(\vec{p}-\vec{p}\,^{\prime});\,\,\,|\vec{p}\,\rangle\langle\vec{p}\,|=\int
d^{3}\vec{p}$ (6)

which implies that in coordinate space

$\langle\vec{r}\,|\vec{p}\,\rangle=\frac{e^{i\vec{p}\cdot\vec{r}}}{(2\pi)^{3/2}}$
(7)

We are going to use the $T$-matrix obtained as

$T=\frac{V}{1-VG}$ (8)

for a single-channel collision, where $G$ is the loop function of the
intermediate particles Gamermann10

$G(E)=\int_{|\vec{q}\,|<q_{\text{max}}}\frac{d^{3}\vec{q}}{E-\omega_{1}(q)-\omega_{2}(q)+i\eta}$
(9)

with $\omega_{i}(q)=\sqrt{{\vec{q}}\,^{2}+m_{i}^{2}}$ the energy of particle
$i$. Eq. (8) is obtained, as shown in Gamermann10 , starting from a
ultraviolet regularized separable constant (s-wave) potential

$V(\vec{p},\vec{p}\,^{\prime})=V\theta(q_{\text{max}}-|\vec{p}\,|)\theta(q_{\text{max}}-|\vec{p}\,^{\prime}|)\
,$ (10)

which reverts into

$T(E;\vec{p},\vec{p}\,^{\prime})=T(E)\theta(q_{\text{max}}-|\vec{p}\,|)\theta(q_{\text{max}}-|\vec{p}\,^{\prime}|).$
(11)

Eq. (8), called on-shell factorization Nieves:1999bx  of the
Lippmann–Schwinger equation, or Bethe–Salpeter equation if relativistic
energies are used, is common in the study of interactions of particles and can
be equally justified using dispersion relations and neglecting the
contribution of the left hand cut (more properly, including this contribution
as an energy independent term OllerUlf01 ).

Eq. (8) is easily generalized to coupled channels as

$T=\Big{[}1-VG\Big{]}^{-1}V$ (12)

where now $V$ is the transition potential $V_{ij}$ between the channels $i$
and $j$ and $G$ is now the diagonal loop function $G\equiv$ diag$[G_{i}(E)]$,
with $G_{i}(E)$ the loop function of each particular channel $i$. Introducing
the identity operators $|\vec{p}\,^{\prime}\rangle\langle\vec{p}\,^{\prime}|$
and $|\vec{p}\,^{\prime\prime}\rangle\langle\vec{p}\,^{\prime\prime}|$ between
the terms $(E-H_{0})^{-1}T\Phi$ in Eq. (5) we obtain

$\langle\vec{r}\,|\widetilde{\Psi}\rangle=\frac{e^{i\vec{p}\cdot\vec{r}}}{(2\pi)^{3/2}}+\theta(q_{\text{max}}-|\vec{p}\,|)T(E)\int\frac{d^{3}\vec{q}}{(2\pi)^{3/2}}\frac{e^{i\vec{q}\cdot\vec{r}}\,\,\theta(q_{\text{max}}-|\vec{q}\,|)}{E-\omega_{1}(q)-\omega_{2}({q})+i\eta}\equiv\frac{\langle\vec{r}\,|\Psi\rangle}{N}\equiv\frac{\Psi(\vec{r},\vec{p}\,)}{N}$
(13)

and the modulus of the relative momenta of the pair of particles is
$|\vec{p}\,|=\lambda^{1/2}(s,m^{2}_{1},m^{2}_{2})/(2\sqrt{s})$ being
$\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc$ the well known Källen function,
and $\sqrt{s}=E$ the total energy in the c.m. frame. Normalizing the wave
function $\Psi(\vec{r},\vec{p}\,)$ to unit flux, as demanded in the
construction of Eq. (1), we fix $N=(2\pi)^{3/2}$ and find

$\Psi(\vec{r},\vec{p}\,)=e^{i\vec{p}\cdot\vec{r}}+\theta(q_{\text{max}}-|\vec{p}\,|)\,T(E)\int_{|\vec{q}\,|<q_{\text{max}}}\frac{d^{3}\vec{q}\
e^{i\vec{q}\cdot\vec{r}}}{E-\omega_{1}(q)-\omega_{2}({q})+i\eta}$ (14)

In the expansion of $e^{i\vec{q}\cdot\vec{r}}$ in terms of spherical
harmonics,

$e^{i\vec{q}\cdot\vec{r}}=4\pi\sum_{\ell=0}^{+\infty}i^{\ell}j_{\ell}(qr)\sum_{m=-\ell}^{+\ell}Y^{*}_{\ell
m}(\hat{q})Y_{\ell m}(\hat{r})$ (15)

only the $Y^{*}_{00}(\hat{q})$ term contributes in the integral and there,
$e^{i\vec{q}\cdot\vec{r}}$ is replaced by the spherical Bessel function
$j_{0}(qr)$,

$\Psi(\vec{r},\vec{p}\,)=e^{i\vec{p}\cdot\vec{r}}+f_{E}(r)\ $ (16)

where the correction to the plane wave incorporates the information on the
interaction, and it depends on the energy, the momentum regulator and the
relative distance $r$,

$\displaystyle f_{E}(r)$ $\displaystyle=$
$\displaystyle\theta(q_{\text{max}}-|\vec{p}\,|)\,T(E)\times$ (17)
$\displaystyle\times$
$\displaystyle\int_{|\vec{q}\,|<q_{\text{max}}}\frac{d^{3}\vec{q}\,j_{0}(qr)}{E-\omega_{1}(q)-\omega_{2}({q})+i\eta}$

Thus, the modulus square of the wave function can be written simply as

$|\Psi(\vec{r},\vec{p}\,)|^{2}=1+|f_{E}(r)|^{2}+2\mbox{Re}\Big{[}e^{i\vec{p}\cdot\vec{r}}f^{*}_{E}(r)\Big{]}\
.$ (18)

When performing $\int
d^{3}\vec{r}S_{12}(r)2\mbox{Re}\Big{[}e^{i\vec{p}\cdot\vec{r}}f^{*}_{E}(r)\Big{]}$,
we can again expand $e^{i\vec{p}\cdot\vec{r}}$ in terms of spherical harmonics
of $\vec{r}$ as done in Eq. (15) and since $S_{12}(r)$ does not depend on
angular variables, we again pick up the $Y_{00}(\hat{r})$ term and therefore

$\displaystyle\int
d^{3}\vec{r}\,S_{12}(r)|\Psi(\vec{r},\vec{p}\,)|^{2}=1+4\pi\int_{0}^{+\infty}drr^{2}\,S_{12}(r)$
$\displaystyle\hskip
28.45274pt\times\Big{\\{}|j_{0}(pr)+f_{E}(r)|^{2}-j^{2}_{0}(pr)\Big{\\}}$ (19)

obtaining the Koonin–Pratt formula Koonin77 ; Pratt90 ; Bauer92 . Note that in
the derivation of Eq. (19) we added and subtracted $j^{2}_{0}(pr)$.

With respect to the Koonin–Pratt formula we obtain (see Eq. (17)) the
$d^{3}\vec{q}$ integral regulated by the maximum momentum $q_{\text{max}}$
demanded by our approach for the $T$-matrix and the factor
$\theta(q_{\text{max}}-|\vec{p}\,|)$ which will be inoperative in the range of
values of the relative momentum $p$ of interest for the correlation function,
with values of $p$ smaller than $q_{\text{max}}$. Note that $f_{E}(r)$ in Eq.
(17) is already regulated by the Bessel function $j_{0}(qr)$ and we shall
investigate this in connection with the range of the interaction introduced by
our formalism.

In the next step we introduce the effect of coupled channels, similarly as
done in Ref. Haidenbauer19 . Let us assume we have $n$ coupled channels and
that the channel $i$ is the one observed in the final state. We can reach this
state starting from any other channel $j\neq i$, and by analogy to Eq. (14)
the corresponding wave function will be

$\displaystyle\Psi_{j}(\vec{r},\vec{p}\,)=T_{ji}(E)\theta(q_{\text{max}}-|\vec{p}\,|)$
$\displaystyle\times\int_{|\vec{q}\,|<q_{\text{max}}}\,\frac{d^{3}\vec{q}\,e^{i\vec{q}\cdot\vec{r}}}{E-\omega^{(j)}_{1}(q)-\omega^{(j)}_{2}({q})+i\eta}$
(20)

where $\omega^{(j)}_{1}(q)$ and $\omega^{(j)}_{2}(q)$ are the kinetic
energies,including masses, of particles $1,2$ in the channel $j$ and $E$ and
$p$ are related through $E=\omega^{(i)}_{1}(p)+\omega^{(i)}_{2}(p)$. Note that
in the $d^{3}\vec{q}$ integral, one can again replace
$e^{i\vec{q}\cdot\vec{r}}$ by $j_{0}(qr)$. In the Koonin–Pratt formalism the
contribution from these channels is added incoherently to the one of the
observed channel, with a weight $w_{j}$ relative to the weight $w_{i}=1$
implicit in Eq. (19) for the observed channel. We can write in general for the
$i$ observed channel

$\displaystyle\Psi_{j}(\vec{r},\vec{p}\,)=\delta_{ij}j_{0}(pr)+T_{ji}(E)\theta(q_{\text{max}}-|\vec{p}\,|)$
$\displaystyle\times\int_{|\vec{q}\,|<q_{\text{max}}}d^{3}\vec{q}\frac{j_{0}(qr)}{E-\omega_{1}^{(j)}(q)-\omega_{2}^{(j)}(q)+i\eta}$
(21)

and the correlation function will read

$\displaystyle C(\vec{p}\,)=1+4\pi\int_{0}^{+\infty}drr^{2}S_{12}(r)$
$\displaystyle\times\left(\sum_{j}w_{j}|\widetilde{\Psi}_{j}(\vec{r},\vec{p}\,)|^{2}-j^{2}_{0}(pr)\right)$
(22)

with the weight $w_{i}$ of the observed channel $i$ equal to $1$. This formula
is also used in Ref. Wei23  to evaluate the correlation function in the case
of the $DK$ interaction with a different normalization which we address below.

In the $T_{cc}$ study of Ref. Feijoo21  a different, relativistic
normalization is used for the $T$-matrix, common in quantum field theory
approaches. Since $VG$ (or $TG$) is weighted versus $1$ in Eq. (12), the
product $TG$ is independent of normalization conventions. Since the two-hadron
loop function in Feijoo21  is given by Oller97

$\displaystyle
G^{\text{FT}}(\sqrt{s}=E)=\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{\omega_{1}(q)+\omega_{2}(q)}{2\omega_{1}(q)\,\omega_{2}(q)}$
$\displaystyle\times\frac{1}{s-\left[\omega_{1}(q)+\omega_{2}({q})\right]^{2}+i\eta}\,\,,$
(23)

comparing this loop function with the one so far used of Eq. (9) we can
immediately write with the new conventions333In Eq. (24), we should have used
the notation $T^{\text{FT}}$, referring to the scattering amplitude obtained
in Ref. Feijoo21 , but for simplicity, the FT tag is not included. In
addition, we will also omit, from now on, the FT label in the loop function
$G$. The relation of $T^{\text{FT}}$ with the $S-$matrix can be found e.g. in
Ref. Gamermann10  (see also next Sec. III).

$\displaystyle\Psi_{j}(\vec{r},\vec{p}\,)=\delta_{ij}j_{0}(pr)+T_{ji}(E)\theta(q_{\text{max}}-|\vec{p}\,|)\times$
$\displaystyle\int_{|\vec{q}\,|<q_{\text{max}}}\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{\omega_{1}^{(j)}(q)+\omega_{2}^{(j)}(q)}{2\omega_{1}^{(j)}(q)\,\omega_{2}^{(j)}(q)}$
$\displaystyle\times\frac{j_{0}(qr)}{s-\left[\omega_{1}^{(j)}(q)-\omega_{2}^{(j)}(q)\right]^{2}+i\eta}~{},$
(24)

and Eq. (22) equally holds with the new conventions. Under the approximation
$s-(\omega_{1}+\omega_{2})^{2}\sim
2(\omega_{1}+\omega_{2})(\sqrt{s}-\omega_{1}-\omega_{2})$ one recovers the
formulas used in Ref. Wei23 .

##  III The $T_{cc}(3875)^{+}$ states with $D^{0}D^{*+}$ and $D^{+}D^{*0}$
channels

In Ref. Feijoo21  the interaction between the $D^{0}D^{*+}$ and $D^{+}D^{*0}$
channels was studied exchanging vector mesons in the extension of the local
hidden gauge symmetry approach Bando88 ; Harada03 ; Meissner88 ; Nagahiro09
to the charmed sector. It was found that the isospin $I=0$ combination had an
attractive interaction which would bind the system, while the $I=1$ one had a
repulsive interaction. This observation agrees with the conclusions of the
experimental analysis Aaij22c . The study of Ref. Feijoo21  was done with the
explicit two channels where the $T_{ij}(\sqrt{s}=E)$ amplitudes were
evaluated. Only fine tuning of the parameter $q_{\text{max}}$ was done to
obtain the experimental binding of Aaij22c , and the width came out
automatically from the consideration of the $D^{*}\rightarrow D\pi$ decay. We
take advantage to mention, since this information will be used below, that
from this work we can find the scattering lengths and effective ranges

$\displaystyle a_{D^{0}D^{*+}}$ $\displaystyle=$
$\displaystyle(7.78-1.82\,i)\,\,\mbox{fm}$ $\displaystyle a_{D^{+}D^{*0}}$
$\displaystyle=$ $\displaystyle(2.07-1.21\,i)\,\,\mbox{fm}$ $\displaystyle
r_{0,D^{0}D^{*+}}$ $\displaystyle=$ $\displaystyle-3.49\,\,\mbox{fm}$
$\displaystyle r_{0,D^{+}D^{*0}}$ $\displaystyle=$
$\displaystyle(-0.06-2.27\,i)\,\,\mbox{fm}$ (25)

where in the evaluation of $r_{0}$ the very small $D^{*}$ width is
neglected,444This is sufficiently accurate for the purpose of this work.
Nevertheless, a complete study of the three-body dynamics in the $T_{cc}$ due
to the finite life time of the $D^{*}$ can be found in Ref. Du22 . There, the
pion exchange between the $D$ and $D^{*}$ mesons and the finite $D^{*}$ width,
are taken into account simultaneously to ensure that three-body unitarity is
preserved. In that work the low-energy expansion of the $DD^{*}$ scattering
amplitude is also performed and the scattering length and effective range are
extracted maintaining finite the width of the $D^{*}$. and where $a$ and
$r_{0}$ are defined from the Quantum Mechanics amplitude in the effective
range expansion around threshold $\sqrt{s_{\text{th}}}$, where $p$ vanishes

$f_{\text{QM}}\simeq\frac{1}{-\frac{1}{a}+\frac{1}{2}r_{0}p^{2}-ip}\ .$ (26)

It is easy to relate our scattering matrix with this one since from Eqs. (12)
and (23) with a real potential one has

$\mbox{Im}\,T^{-1}=-\mbox{Im}\,G=\frac{p}{8\pi\sqrt{s}}$ (27)

On the other hand

$\mbox{Im}\,f^{-1}_{\text{QM}}=-p\ .$ (28)

Hence

$T=-8\pi\sqrt{s}f_{\text{QM}}\Rightarrow
a=\frac{T(\sqrt{s_{\text{th}}})}{8\pi\sqrt{s}}\ .$ (29)

The values reported in Eq. (25) for the scattering length agree well with
those extracted from experiment in Aaij22c

$\displaystyle a^{\text{exp}}_{D^{0}D^{*+}}$ $\displaystyle=$
$\displaystyle(7.16\pm 0.51-1.85\pm 0.28\,i)\,\,\mbox{fm}$ $\displaystyle
a^{\text{exp}}_{D^{+}D^{*0}}$ $\displaystyle=$
$\displaystyle(1.76-1.82\,i)\,\,\mbox{fm}\ ,$

with similar relative errors expected for $a^{\text{exp}}_{D^{+}D^{*0}}$.
These values also agree with those obtained in Refs. Albaladejo22 ; Du22 .

In Ref. Kamiya22  a scheme is used in which the potentials are written in
coordinate space and the formalism is conducted using explicitly wave
functions as solutions of the Schödinger equation. The range of the potential
is assumed to be given by one pion exchange, which is drastically different
from the one used in Feijoo21 , where the interaction originates from $\rho$
exchange and the range in momentum space is given by $q_{\rm max}\sim 420$ MeV
fitted to data. We shall see that our results for the correlation function are
very similar, which might be surprising in view of the different ranges
assumed for the interaction. We shall find an explanation for this surprising
result in the fact that the strength of the potentials in Kamiya22  is
taken555 Note that in the approach of Ref. Feijoo21 , the $DD^{*}$ contact
interaction is fixed to the prediction of the hidden gauge model. such as to
fit the experimental scattering lengths reported in Eq.
(LABEL:eq:scattlen_exp), and the correlation functions in the present case are
practically only sensitive to the scattering length.

To finish this section let us write explicitly the correlation function for
the $D^{0}D^{*+}$ and $D^{+}D^{*0}$ coupled-channel problem that we are
discussing,

$\displaystyle
C_{D^{0}D^{*+}}(p_{D^{0}})=1+4\pi\,\theta(q_{\text{max}}-p_{D^{0}})\times$
$\displaystyle\int_{0}^{+\infty}drr^{2}S_{12}(r)\Big{\\{}\big{|}j_{0}(p_{D^{0}}r)+T_{11}(E)\widetilde{G}^{(1)}(r;E)\big{|}^{2}$
$\displaystyle+\big{|}T_{21}(E)\widetilde{G}^{(2)}(r;E)\big{|}^{2}-j^{2}_{0}(p_{D^{0}}r)\Big{\\}}$
(31)

for detecting the $D^{0}D^{*+}$ pair, with
$E=\sqrt{s}=\left[\sqrt{m_{D^{0}}^{2}+p_{D^{0}}^{2}}+\sqrt{m_{D^{*+}}^{2}+p_{D^{0}}^{2}})\right]$,
and

$\displaystyle
C_{D^{+}D^{*0}}(p_{D^{+}})=1+4\pi\,\theta(q_{\text{max}}-p_{D^{+}})\times$
$\displaystyle\int_{0}^{+\infty}drr^{2}S_{12}(r)\Big{\\{}\big{|}j_{0}(p_{D^{+}}r)+T_{22}(E)\widetilde{G}^{(2)}(r;E)\big{|}^{2}$
$\displaystyle+\big{|}T_{12}(E)\widetilde{G}^{(1)}(r;E)\big{|}^{2}-j^{2}_{0}(p_{D^{+}}r)\Big{\\}}$
(32)

for observing $D^{+}D^{*0}$, where now the c.m. energy is given by
$E=\left[\sqrt{m_{D^{+}}^{2}+p_{D^{+}}^{2}}+\sqrt{m_{D^{*0}}^{2}+p_{D^{+}}^{2}})\right]$.

Finally, the quantities $\widetilde{G}^{(i)}$ are simply

$\displaystyle\widetilde{G}^{(i)}(r;E)$
$\displaystyle=\int_{|\vec{q}\,|<q_{\text{max}}}\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{\omega^{(i)}_{1}(q)+\omega^{(i)}_{2}(q)}{2\omega^{(i)}_{1}(q)\,\omega^{i}_{2}(q)}$
$\displaystyle\times\frac{j_{0}(qr)}{s-\left[\omega_{1}^{(i)}(q)+\omega_{2}^{(i)}({q})\right]^{2}+i\eta}$
(33)

and the $T$-matrix elements, $T_{ij}(E)$, are taken from Feijoo21 . We have
taken $(w_{1}=1)$ for the observed channel, as we discussed, and for the other
channel, given the analogy between $D^{0}D^{*+}$ and $D^{+}D^{*0}$, we have
also taken the weight equal to one, $w_{2}=1$.

##  IV Results

Figure 1: (Color online). Correlation function of the $D^{0}D^{*+}$ pair for
different values of the source size. Figure 2: (Color online). Correlation
function of the $D^{+}D^{*0}$ pair for different values of the source size.

In Figs. 1 and 2 we show our results for the correlation function for the
$D^{0}D^{*+}$ and $D^{+}D^{*0}$ channels respectively, and different values of
the source size $R=1,2,3,5$ fm. We can see that for the case $D^{0}D^{*+}$,
the correlation function at very small $p$ and $R=1$ fm is quite large
compared to 1. This is also what is found in Ref. Kamiya22 . As the source
size $R$ is increased, the correlation function becomes closer to unity at
smaller momenta. We also observe a fast convergence of the correlation
function to unity for values of $p$ around 200 MeV. As anticipated, the factor
$\theta(q_{\text{max}}-|\vec{p}\,|)$ in Eq. (24) is inoperative, since
$q_{\text{max}}=420$ MeV in Feijoo21  and we only go in both figures up to 300
MeV. Also in the case of the $D^{+}D^{*0}$, the correlation function obtained
is comparable with the one calculated in Kamiya22 , and features similar to
the ones just described for the correlation function of the $D^{0}D^{*+}$ pair
are also observed in this case. As we shall see below, the $T_{cc}$
correlation functions are almost determined by the scattering lengths and
unitarity. Thus, although the interaction range in coordinate space of
Kamiya22  is much larger than that in the approach of Ref. Feijoo21 , the
$DD^{*}$ isoscalar interaction strength (the sharp ultraviolet cutoff
$q_{\text{max}}$) is adjusted in the former (latter) work to reproduce the
experimental $T_{cc}$ mass. Thus both approaches, the one of Ref. Kamiya22
and that of Ref. Feijoo21 , provide similar scattering lengths and hence
correlation functions. In Sec. V we will also compute the correlation
functions using the $T$-matrix obtained in Ref. Albaladejo22 , finding also
similar results due to the same argument.

Figure 3: (Color online). Correlation function of the $D^{0}D^{*+}$ pair
obtained with the $T$-matrix of Ref. Feijoo21  and different approximations in
terms of the scattering length and the effective range. The source size is
taken $R=1$ fm. Figure 4: (Color online). Correlation function of the
$D^{+}D^{*0}$ pair obtained with the $T$-matrix of Ref. Feijoo21  and
different levels of approximation in terms of the scattering length and the
effective range. The source size is taken $R=1$ fm.

The next step is to see to which magnitudes of the scattering matrix the
correlation function is sensitive. For this purpose we take $R=1$ fm and show
in Figs. 3 and 4 the results using $T=-8\pi\sqrt{s}/(-1/a)$ (blue dashed
lines), $T=-8\pi\sqrt{s}/(-1/a-ip)$ (green long dashed lines) and
$T=-8\pi\sqrt{s}/(-1/a+r_{0}p^{2}/2-ip)$ (brown dotted lines) and compare the
results with those obtained using the full $T$-matrix evaluated in Feijoo21 .
Note that since the scattering length and the effective range are defined only
for the diagonal channels $D^{0}D^{*+}\rightarrow D^{0}D^{*+}$ and
$D^{+}D^{*0}\rightarrow D^{+}D^{*0}$, in Figs. 3 and 4 for the calculation of
the correlation function with the $T$-matrix of Feijoo21  we have taken,
respectively, the weights of the wave functions $\Psi_{D^{+}D^{*0}\rightarrow
D^{0}D^{*+}}$ and $\Psi_{D^{0}D^{*+}\rightarrow D^{+}D^{*0}}$ equal to zero
($w_{2}=0$). This is to say, $T_{21}$ and $T_{12}$ effects are neglected in
Eqs. (31) and (32), respectively. We can see that for the approximation
$T\propto a$ we only get the result of the correlation function at threshold,
as expected, but as $p$ increases we see that the correlation function does
not converge to 1 as it should. However, when we use $T\propto(-1/a-ip)^{-1}$
the agreement with the exact calculation is very good, with differences which
are smaller than ordinary experimental errors, which indicates that the
correlation function in both cases gives us information mostly about the
scattering length. Yet, for the correlation function of the $D^{0}D^{*+}$
channel, the consideration of the effective range $r_{0}$ helps to improve the
agreement with the exact result. The message then is that, provided one finds
experimental input to choose the right parameter of the source, $R$, we could
then get the scattering length for the $D^{0}D^{*+}$ and $D^{+}D^{*0}$
channels, but probably less information about higher order parameters in the
effective range expansion. Determining the scattering length, of course, is a
very valuable information, and it seems a rather general feeling about what
can be accomplished with correlation functions, although for particular cases
one may get a different conclusion.

It is interesting to see where the contribution to the correlation function
comes from. For this purpose in Figs. 5 and 6 we show the results removing in
Eqs. (31) and (32), respectively, the first channel
$|T_{11}\widetilde{G}^{(1)}|^{2}$ and $|T_{22}\widetilde{G}^{(2)}|^{2}$ (blue
dashed lines), the second channel, $|T_{21}\widetilde{G}^{(2)}|^{2}$ and
$|T_{12}\widetilde{G}^{(1)}|^{2}$ (green long dashed lines), and also keeping
only the interference term between the incoming and the outgoing wave of the
observed channel (dotted lines). This is subtracting
$|T_{11}\widetilde{G}^{(1)}|^{2}$ and $|T_{21}\widetilde{G}^{(2)}|^{2}$ from
Eq. (31) and $|T_{22}\widetilde{G}^{(2)}|^{2}$ and
$|T_{12}\widetilde{G}^{(1)}|^{2}$ in Eq. (32). We observe that taking into
account the re-scattering term in the observed channel is essential to go from
negative $C(p)$ to positive at small values of $p$ in the case of
$D^{0}D^{*+}$. The contribution of the re-scattering of the non observed
channel is also important but less than the rescattering of the observed
channel.

Figure 5: (Color online). Different contributions to the correlation function
of the pair $D^{0}D^{*+}$. The source size is taken $R=1$ fm. Figure 6: (Color
online). Different contributions to the correlation function of the pair
$D^{+}D^{*0}$. The source size is taken $R=1$ fm.

Finally, we present here another test. One novelty of our approach is that we
regularize $\widetilde{G}$ with the same cutoff that is used to regularize $G$
when calculating the scattering in Feijoo21 . Although for consistency we have
to use the same cutoff in $G$ and $\widetilde{G}$ as we have shown in the
derivation, it is instructive to see how sensitive is the correlation function
to the parameter $q_{\text{max}}$, which as shown in Eq. (10) provides the
range of the interaction in momentum space. We show in Figs. 7 and 8 what
happens if we increase this cutoff in the evaluation of $\widetilde{G}$ in Eq.
(33), while keeping the $T_{ij}$ matrices of Feijoo21  unchanged. What we see
is that the results are rather insensitive to the value of the cutoff,
particularly in the case of the correlation function of the $D^{0}D^{*+}$
pair, which again shows the limitations of the correlation function to provide
information on this magnitude. Note that we are changing $q_{\text{max}}$ only
in $\widetilde{G}$, but $q_{\text{max}}$ enters the evaluation of the
scattering matrix. This has also a positive reading, since if one wishes to
obtain the $T_{ij}$-matrix from experimental correlation functions, the range
of the interaction needed in the evaluation of $\widetilde{G}$ does not matter
much.

Figure 7: (Color online). Dependence of the correlation function of the pair
$D^{0}D^{*+}$ on the cutoff parameter $q_{\text{max}}$. The source size is
taken $R=1$ fm. Figure 8: (Color online). Dependence of the correlation
function of the pair $D^{+}D^{*0}$ on the cutoff parameter $q_{\text{max}}$.
The source size is taken $R=1$ fm.

It is interesting to see the reason for this relative lack of sensitivity. Let
us take for instance the term
$2j_{0}(p_{D^{0}}r)\mbox{Re}[T_{11}\widetilde{G}^{(1)}]$ in Eq. (31). We have

$\int d^{3}\vec{q}\int drr^{2}4\pi
S_{12}(r)j_{0}(p_{D^{0}}r)j_{0}(qr)\,\cdot\cdot\cdot\cdot$ (34)

We can write

$j_{0}(qr)=\frac{e^{iqr}-e^{-iqr}}{2iqr}$ (35)

and the same for $j_{0}(p_{D^{0}}r)$ and then perform the integration over
$r$. We obtain some polynomial and most importantly the exponential factor

$e^{-Q^{2}R^{2}}$ (36)

with $Q=p_{D^{0}}+q$ and $Q=p_{D^{0}}-q$. With a value of $R=1$ fm this factor
kills values of $q$ of the order of $R^{-1}\approx 200$ MeV in the integration
over $q$. The Bessel function acts as a regulator of the $d^{3}\vec{q}$
integral and if we cut it with a value of $q_{\text{max}}$ reasonably larger
than $200$ MeV this extra regulator becomes inoperative. In other words, the
correlation function will provide little information on the range of the
interaction in momentum space, unless this is reasonably smaller than
$R^{-1}$. This also means that with values of $R$ of the order of $5$ fm any
possible information on the range of the interaction is greatly lost. The
argumentation would proceed similarly with the quadratic terms in the
scattering matrix, where we would have $\int d^{3}\vec{q}\int d^{3}q^{\prime}$
and the same factor of Eq. (36) with $Q=q+q^{\prime},Q=q-q^{\prime}$.
Certainly, the correlation functions can provide information on the $T_{ij}$
matrix which depends on the range of the interaction through $q_{\text{max}}$
in $G$ of Eq. (12), but since we can get similar $T$-matrices by a
simultaneous change of the strength of the potential and $q_{\text{max}}$,
such that $\delta(VG)=0$, once again, it is not easy to obtain information on
this interaction range, as illustrated above in Figs. 3 and 4.

##  V Comparison with other schemes

Figure 9: Correlation functions for the $D^{0}D^{\ast+}$ (left) and
$D^{+}D^{\ast 0}$ (right) pairs. The solid lines, denoted by “Model A”, show
the correlation functions computed as discussed in Sec. IV, using as an input
the amplitudes derived in Ref. Feijoo21 . The dash-dotted line, denoted by
“Model B”, represent the same correlation functions but calculated with the
modifications discussed in Sec. V, using as an input the $T$-matrix from Ref.
Albaladejo22 . Red, blue, and green lines correspond to a source size
$R_{0}=1$, $3$, and $5\,\text{fm}$, respectively. The vertical, magenta,
dotted line on the left plot indicates the $D^{0}D^{\ast+}$ threshold.

In Sec. IV, the correlation functions for the $D^{0}D^{\ast+}$ and
$D^{+}D^{\ast 0}$ pairs (in the $T_{cc}^{+}$ channel) have been computed with
the formalism discussed in Sec. III. There, we have taken as an input the
scattering matrix of Ref. Feijoo21 , which provides a reasonable reproduction
of the experimental mass and width of the $T_{cc}^{+}$ state, and thus
constitutes a good baseline for our study. In this section, and in order to
have an estimation of the theoretical uncertainties in our calculation, we
employ instead the amplitudes calculated in Ref. Albaladejo22 . In this latter
work, the experimental LHCb $D^{0}D^{0}\pi^{+}$ spectrum where the
$T_{cc}^{+}$ peak appears is fitted. The main difference of Ref. Albaladejo22
with respect to the model of Ref. Feijoo21  is that the interaction potential
is not derived from the exchange of vector mesons dictated by the hidden gauge
symmetry approach, but two independent constants appearing in a general-
isospin symmetric potential are fitted to the experimental spectrum.
Furthermore, the loop functions entering the $T$-matrix are not regulated with
a sharp cutoff, but with a Gaussian one, and non-relativistic kinematics are
explicitly used. Thus, in the alternative calculation of this section, we use
the scattering amplitudes obtained in Ref. Albaladejo22  and, for consistency,
we employ non-relativistic kinematics and a Gaussian regulator in the
computation of the quantities $\widetilde{G}_{i}(r;E)$.

The correlation functions obtained with this new scheme are displayed in Fig.
9 with dash-dotted lines, and compared with those reported in Sec. IV, which
are shown with solid lines. As can be seen, there is very good agreement for
the $D^{0}D^{\ast+}$ correlation function for the three sizes considered
($R_{0}=1$ (red), $3$ (blue), and $5\,\text{fm}$ (green)). This agreement is
due to the similar $T_{cc}^{+}$ pole description found in Refs. Feijoo21 ;
Albaladejo22 . For the $D^{+}D^{\ast 0}$ pair there are some discrepancies,
but the trend is very similar in both cases and, more importantly, this
channel is further away from the $T_{cc}^{+}$ pole than the $D^{0}D^{\ast+}$
one, which is the most relevant one to obtain the $T_{cc}^{+}$ properties.

##  VI Conclusions

We have carried out a detailed study of the correlation function $C(p)$ for
the $D^{0}D^{*+}$ and $D^{+}D^{*0}$ channels that build up the
$T_{cc}(3875)^{+}$ state, only a few hundred keV below the $D^{0}D^{*+}$
threshold. We find values of $C(p)$ positive and very large, around 30 for a
source size $R=1$ fm, at the threshold for the $D^{0}D^{*+}$ channel and also
positive and bigger than 1 for the $D^{+}D^{*0}$ channel, but substantially
smaller than for $D^{0}D^{*+}$. Both correlation functions converge rapidly to
unity for values of the relative momentum of the pair, $p$, around 200 MeV. We
have performed several tests to see which are the terms that contribute most
to the correlation function and found that all of them, interference of
incoming wave with the re-scattered wave, and the squared of the re-scattered
one, are important, particularly that of the observed channel. We modified
slightly the formalism of Koonin–Pratt and showed that it is possible to
factorize the scattering amplitude outside the integrals and only the range of
the interaction has an effect in regularizing the loop function entering the
formalism. Even then, we proved that with values of $R=1-5$ fm for the
Gaussian source function, the range of the strong interaction becomes almost
inoperative in the present case, once the proper scattering matrices are used.
This is because the Bessel functions entering the formalism provide a stronger
cut in the loop integration than the one provided by the range of the strong
interaction. We conducted another important test which is that the
substitution of the $T$-matrix by $T\propto(-1/a-ip)^{-1}$ in the effective
range expansion, already provides an excellent approximation to the exact
solution, although for the case of the $D^{0}D^{*+}$ channel, the
consideration of the effective range improves the agreement. This makes us
reach the conclusion that one can only hope to accurately get the scattering
length from these reactions, provided that from phenomenology we can get a
hold on a reasonable value of the parameter of the Gaussian source function.
Even then, the possibility of measuring correlation functions for pairs of
particles that cannot be reached in scattering experiments, opens the door to
find very valuable information on the interaction of hadrons. There is still a
caveat for the case of coupled channels, since we found that the contribution
of the non observed channel to the correlation function is not small. In other
words, for the two channel case the correlations depend on
$T_{11},T_{12},T_{22}$ and we only have two correlation functions, so we
cannot determine the three magnitudes from two observables. Yet, by using the
optical theorem one can relate $|T_{12}|^{2}$ with the imaginary part of
$T_{22}$ and one can exploit this feature. In any case, the value of the
correlation function could serve to test models rather than to build them. In
the present case, the measurement of the correlation functions for the
channels of the $T_{cc}$ would provide additional information to contrast what
has been obtained from the analysis of the data on the $D^{0}D^{0}\pi^{+}$
spectrum observed in the LHCb experiment.

## Acknowledgments

This work was supported by the Spanish Ministerio de Ciencia e Innovación
(MICINN) and European FEDER funds under Contracts No. PID2020-112777GB-I00,
and by Generalitat Valenciana under contract PROMETEO/2020/023. This project
has received funding from the European Union Horizon 2020 research and
innovation programme under the program H2020-INFRAIA-2018-1, grant agreement
No. 824093 of the STRONG-2020 project. M. A. and A. F. are supported through
Generalitat Valencia (GVA) Grants Nos. CIDEGENT/2020/002 and APOSTD-2021-112,
respectively. M. A. and A. F. thank the warm support of ACVJLI.

## References

  * (1) J. A. Oller, E. Oset and A. Ramos, Prog. Part. Nucl. Phys. 45, 157-242 (2000).
  * (2) F. K. Guo, C. Hanhart, U. G. Meißner, Q. Wang, Q. Zhao, and B. S. Zou, Rev. Mod. Phys. 90, no.1, 015004 (2018) [erratum: Rev. Mod. Phys. 94, no. 2, 029901 (2022)].
  * (3) X. K. Dong, F. K. Guo and B. S. Zou, Commun. Theor. Phys. 73, no. 12, 125201 (2021).
  * (4) X. K. Dong, F. K. Guo and B. S. Zou, Few Body Syst. 62, no. 3, 61 (2021).
  * (5) X. K. Dong, F. K. Guo and B. S. Zou, Progr. Phys. 41, 65-93 (2021).
  * (6) R. Aaij et al., LHCb, Phys. Rev. Lett. 115, 072001 (2015), arXiv:1507.3414 [hep-ex].
  * (7) R. Aaij et al., LHCb, Chin. Phys. C 40, 011001 (2016), arXiv:1509.00292 [hep-ex].
  * (8) R. Aaij et al., LHCb, Phys. Rev. Lett. 122, 222001 (2019), arXiv:1904.03947 [hep-ex].
  * (9) R. Aaij et al., (LHCb), Sci. Bull. 66, 1278 (2021), arXiv:2012.10380 [hep-ex].
  * (10) R. Aaij et al., (LHCb), (2022), arXiv:2210.10346 [hep-ex].
  * (11) R. Aaij et a., (LHCb), Phys. Rev. Lett. 125, 242001 (2020).
  * (12) R. Aaij et a., (LHCb), Phys. Rev. D 102, 112003 (2020).
  * (13) R. Aaij et al., (LHCb Collaboration), Nat. Phys. 18, 751 (2022).
  * (14) R. Aaij et al. (LHCb Collaboration), Nat. Commun. 13, 3351 (2022).
  * (15) A. Feijoo, W. H. Liang, and E. Oset, Phys. Rev. D 104, 114015 (2021).
  * (16) X. Z. Ling, M. Z. Liu, L. S. Geng, E. Wang, and J. J. Xie, Phys. Lett. B 826, 136897 (2022).
  * (17) S. Fleming, R. Hodges, and T. Mehen, Phys. Rev. D 104, 116010 (2021).
  * (18) H. Ren, F. Wu, and R. Zhu, Adv. High Energy Phys. 2022, 9103031 (2022).
  * (19) K. Chen, R. Chen, L. Meng, B. Wang, and S. L. Zhu, Eur. Phys. J. C 82, 581 (2022).
  * (20) M. Albaladejo, Phys. Lett. B 829, 137052 (2022)
  * (21) M. L. Du, V. Baru, X. K. Dong, A. Filin, F. K. Guo, C. Hanhart, A. Nefediev, J. Nieves, and Q. Wang, Phys. Rev. D 105, 014024 (2022).
  * (22) V. Baru, X. K. Dong, M. L. Du, A. Filin, F. K. Guo, C. Hanhart, A. Nefediev, J. Nieves and Q. Wang, Phys. Lett. B 833 137290 (2022).
  * (23) N. Santowsky and C. S. Fischer, Eur. Phys. J. C 82, 313 (2022).
  * (24) C. Deng and S. L. Zhu, Phys. Rev. D 105, 054015 (2022).
  * (25) H. W. Ke, X. H. Liu, and X. Q. Li, Eur. Phys. J. C 82, 144 (2022).
  * (26) S. S. Agaev, K. Azizi, and H. Sundu, J. High Energy Phys. 06, 057 (2022).
  * (27) Y. Kamiya, T. Hyodo, and A. Ohnishi, Eur. Phys. J. A 58, 131 (2022).
  * (28) L. Meng, B. Wang, G. J. Wang, and S. L. Zhu, arXiv: 2204.08716.
  * (29) L. M. Abreu, Nucl. Phys. B 985, 115994 (2022).
  * (30) S. Chen, C. Shi, Y. Chen, M. Gong, Z. Liu, W. Sun, and R. Zhang, Phys. Lett. B 833, 137391 (2022).
  * (31) M. Albaladejo and J. Nieves, Eur. Phys. J. C 82, 724 (2022).
  * (32) J. I. Ballot and J. M. Richard, Phys. Lett. B 123, 449 (1983).
  * (33) S. Zouzou, B. Silvestre-Brac, C. Gignoux and J. M. Richard, Z. Phys. C 30, 457 (1986).
  * (34) K. Morita, T. Furumoto, and A. Ohnishi, Phys. Rev. C 91, 024916 (2015).
  * (35) A. Ohnishi, K. Morita, K. Miyahara, and T. Hyodo, Nucl. Phys. A 954, 294 (2016).
  * (36) K. Morita, A. Ohnishi, F. Etminan, and T. Hatsuda, Phys. Rev. C 94, 031901 (2016).
  * (37) T. Hatsuda, K. Morita, A. Ohnishi, and K. Sasaki, Nucl. Phys. A 967, 856 (2017).
  * (38) D. L. Mihaylov et al., Eur. Phys. J. C 78, 394 (2018).
  * (39) J. Haidenbauer, Nucl. Phys. A 981, 1 (2019).
  * (40) K. Morita et al., Phys. Rev. C 101,015201(2020).
  * (41) Y. Kamiya,T. Hyodo,K. Morita,A. Ohnishi, and W. Weise, Phys. Rev. Lett. 124, 132501 (2020).
  * (42) Y. Kamiya, K. Sasaki, T. Fukui, T. Hyodo, K. Morita, K. Ogata, A. Ohnishi, and T. Hatsuda, Phys. Rev. C 105, 014915 (2022).
  * (43) L. Adamczyk et al., STAR. Phys. Rev. Lett. 114, 022301 (2015).
  * (44) S. Acharya et al., ALICE. Phys. Lett. B 774, 64 (2017).
  * (45) J. Adam et al., STAR. Phys. Lett. B 790, 490 (2019).
  * (46) S. Acharya et al., ALICE. Phys. Rev. C 99, 024001 (2019).
  * (47) S. Acharya et al., ALICE. Phys. Rev. Lett. 123, 112002 (2019).
  * (48) S. Acharya et al., ALICE. Phys. Lett. B 797, 134822 (2019).
  * (49) S. Acharya et al., ALICE. Phys. Lett. B 805, 135419 (2020).
  * (50) S. Acharya et al., ALICE. Phys. Rev. Lett. 124, 092301 (2020).
  * (51) S. Acharya et al., ALICE. Nature 588, 232 (2020).
  * (52) S. Acharya et al., ALICE. Phys. Lett. B 822, 136708 (2021).
  * (53) S. Acharya et al., ALICE. Phys. Rev. Lett. 127, 172301 (2021).
  * (54) L. Fabbietti, V. Mantovani Sarti, O. Vázquez Doce, Ann. Rev. Nucl. Part. Sci. 71, 377 (2021).
  * (55) S. Acharya et al., ALICE, arXiv:2201.05352 (2022).
  * (56) Z.-W. Liu, J.-X. Lu, and L.-S. Geng, arXiv:2302.01046 (2023).
  * (57) M. Bando, T. Kugo, and K. Yamawaki, Phys. Rept. 164, 217 (1988).
  * (58) M. Harada and K. Yamawaki, Phys. Rept. 381, 1 (2003).
  * (59) U. G. Meißner, Phys. Rept. 161, 213 (1988).
  * (60) H. Nagahiro, L. Roca, A. Hosaka, and E. Oset, Phys. Rev. D 79, 014015 (2009)
  * (61) M. A. Lisa, S. Pratt, R. Soltz, and U. Wiedmann, Annu. Rev. Nucl. Part. Sci. 55, 357 (2005=.
  * (62) J. Nieves and E. Oset, Phys. Rev. C 47, 1478 (1993).
  * (63) D. Gamermann, J. Nieves, E. Oset, and E. Ruiz Arriola, Phys. Rev. D 81, 014029 (2010).
  * (64) J. Nieves and E. Ruiz Arriola, Nucl. Phys. A 679, 57-117 (2000).
  * (65) J. A. Oller and U. G. Meissner, Phys. Lett. B 500, 263 (2001).
  * (66) S. E. Koonin, Phys. Lett. B 70, 43 (1977).
  * (67) S. Pratt, T. Csörgő, and J. Zimányi, Phys. Rev. C 42, 2646 (1990).
  * (68) W. Bauer, C. Gelbke, and S. Pratt. Annu. Rev. Nucl. Part. Sci. 42, 77 (1992).
  * (69) J. A. Oller and E. Oset, Nucl. Phys. A 620, 438 (1997).