# Mass behavior of hidden-charm pentaquarks with open-strange inspired by
these established $P_{c}$ molecular states

Rui Chen1<EMAIL_ADDRESS>Xiang Liu2,3,4<EMAIL_ADDRESS>1Key
Laboratory of Low-Dimensional Quantum Structures and Quantum Control of
Ministry of Education, Department of Physics and Synergetic Innovation Center
for Quantum Effects and Applications, Hunan Normal University, Changsha
410081, China
2Lanzhou Center for Theoretical Physics, Key Laboratory of Theoretical Physics
of Gansu Province, and Frontier Science Center for Rare Isotopes, Lanzhou
University, Lanzhou 730000, China
3School of Physical Science and Technology, Lanzhou University, Lanzhou
730000, China
4Research Center for Hadron and CSR Physics, Lanzhou University and Institute
of Modern Physics of CAS, Lanzhou 730000, China

###### Abstract

Stimulated by the meson-baryon molecular interpretations of the $P_{c}$ states
($P_{c}(4312)/P_{c}(4440)/P_{c}(4457)$), we systematically study the
interactions between an $S-$wave charm-strange baryon $\Xi_{c}^{(\prime,*)}$
and an anti-charmed meson $\bar{D}^{(*)}$ in a coupled channel analysis.
Effective potentials for the $\Xi_{c}^{(\prime,*)}\bar{D}^{(*)}$ interactions
in a one-boson-exchange model can be related to those in the
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$ systems by using the $SU(3)$ flavor symmetry
and heavy quark symmetry. Our results can predict several promising hidden-
charm molecular pentaquarks with strangeness $|S|=1$, which include the
$\Xi_{c}^{\prime}\bar{D}$ states with $I(J^{P})=0,1(1/2^{-})$, the
$\Xi_{c}^{*}\bar{D}$ states with $0,1(3/2^{-})$, the
$\Xi_{c}^{\prime}\bar{D}^{*}$ states with $0(1/2^{-})$ and $0,1(3/2^{-})$, and
the $\Xi_{c}^{*}\bar{D}^{*}$ states with $0(1/2^{-},3/2^{-},5/2^{-})$.

##  I Introduction

In 2015, the LHCb Collaboration reported two hidden-charm pentaquarks namely
$P_{c}(4380)$ and $P_{c}(4450)$ in the $\Lambda_{b}^{0}\to J/\psi pK^{-}$
decay process Aaij:2015tga . Their observations immediately inspired theorists
propose several different interpretations of the $P_{c}(4380)$ and
$P_{c}(4450)$ states, i.e., the molecular state assignments Chen:2015loa ;
Chen:2015moa ; Karliner:2015ina ; Roca:2015dva ; Mironov:2015ica ; He:2015cea
; Meissner:2015mza ; Burns:2015dwa ; Shimizu:2016rrd ; Chen:2016heh ;
Eides:2015dtr ; Huang:2015uda ; Chen:2016otp ; Yang:2015bmv ; He:2016pfa ;
Yamaguchi:2016ote , the diquark-diquark-antiquark configuration Maiani:2015vwa
; Li:2015gta ; Ghosh:2015ksa ; Anisovich:2015zqa ; Wang:2015ava , the diquark-
triquark configuration Lebed:2015tna ; Zhu:2015bba , the re-scattering effect
Guo:2015umn ; Liu:2015fea ; Mikhasenko:2015vca , and so on (see reviews
Chen:2016qju ; Liu:2019zoy ; Brambilla:2019esw ; Guo:2017jvc ;
Esposito:2016noz ; Hosaka:2016pey  for details). Because the $P_{c}(4380)$ and
$P_{c}(4450)$ are very close to the mass threshold of a charmed baryon and an
anti-charmed meson, the hadronic molecular state assignments to the $P_{c}$
states are the most popular proposal, which already had been predicted
previously Yang:2011wz ; Wu:2010jy ; Wang:2011rga ; Karliner:2015ina ;
Li:2014gra ; Wu:2010vk .

In 2019, the LHCb Collaboration updated the observations of the $P_{c}$ states
in the same process with data collected in run 1 and run 2 Aaij:2019vzc . They
discovered three narrow structures ($P_{c}(4312)^{+}$, $P_{c}(4440)^{+}$, and
$P_{c}(4457)^{+}$) in the $J/\psi p$ invariant mass spectrum, where the
$P_{c}(4440)^{+}$ and $P_{c}(4457)^{+}$ correspond to the fine structures of
the former $P_{c}(4450)$ Aaij:2015tga . The refinement observations of the
$P_{c}$ states may provide a strong evidence of the hidden-charm molecular
pentaquarks Chen:2019asm ; Liu:2019tjn ; He:2019ify ; Meng:2019ilv ;
Burns:2019iih ; Anwar:2018bpu , although there were the other discussions on
these three $P_{c}$ states Weng:2019ynv ; Ali:2019npk ; Giron:2019bcs ;
Cheng:2019obk ; Mutuk:2019snd . In particular, the $P_{c}(4312)$,
$P_{c}(4440)$ and $P_{c}(4457)$ can be assigned as the loosely bound
$\Sigma_{c}\bar{D}$ state with $(I=1/2,J^{P}=1/2^{-})$, the
$\Sigma_{c}\bar{D}^{*}$ state with $(I=1/2,J^{P}=1/2^{-})$ and the
$\Sigma_{c}\bar{D}^{*}$ state with $(I=1/2,J^{P}=3/2^{-})$, respectively,
based on the one-boson-exchange potentials with considering the coupled-
channel effect Chen:2019asm .

Stimulated by the hidden-charm molecular assignments to the $P_{c}$ states,
many groups further search for the possible hidden-charm molecular pentaquarks
with strangeness $|S|=1$ Chen:2016ryt ; Wang:2019nvm ; Wang:2015wsa ;
Feijoo:2015kts ; Lu:2016roh ; Xiao:2019gjd ; Zhang:2020cdi ; Shen:2020gpw ;
Ferretti:2020ewe ; Zhu:2021lhd ; Xiao:2021rgp , $|S|=2$ Wang:2020bjt , and
$|S|=3$ Wang:2021hql . For example, in Ref. Chen:2016ryt , we study the single
$\Lambda_{c}\bar{D}_{s}^{*}$, $\Sigma_{c}\bar{D}_{s}^{*}$,
$\Sigma_{c}^{*}\bar{D}_{s}^{*}$, $\Xi_{c}\bar{D}^{*}$,
$\Xi_{c}^{\prime}\bar{D}^{*}$, and $\Xi_{c}^{*}\bar{D}^{*}$ interactions by
considering the one-eta-exchange and/or one-pion-exchange processes, and
predict several possible strange hidden-charm molecular pentaquarks.

Experimental, very recently, the LHCb Collaboration further analyzed the
$\Xi_{b}^{-}\to J/\psi\Lambda K^{-}$ decay process, and reported an evidence
of the strange hidden-charm pentaquark $P_{cs}(4459)$ in the $J/\psi\Lambda$
invariant mass spectrum LHCb:2020jpq . Its resonance parameters are

$\displaystyle M=4458.8\pm 2.9_{-1.2}^{+4.7}\text{MeV},\quad\Gamma=17.3\pm
6.5_{-5.7}^{+8.0}\text{MeV},$

respectively. Its spin-parity is not determined, yet. Since its mass is just
below the $\Xi_{c}\bar{D}^{*}$ threshold with around 19 MeV, the
$P_{cs}(4459)$ was explained as a strange hidden-charm $\Xi_{c}\bar{D}^{*}$
molecule Chen:2020kco ; Chen:2021tip ; Chen:2020uif ; Peng:2020hql ;
Wang:2020eep .

All these $P_{c}$ and $P_{cs}$ states are a little below the mass thresholds
of a pair of a charmed/or charm-strange baryon and an anti-charmed meson. If
they are the hidden-charm molecules, in this line, one can expect the
existence of more hidden-charm molecular pentaquarks based on the heavy quark
symmetry and the $SU(3)$ flavor symmetry. In this work, we will perform a
systematic investigation on the $\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}$
interactions, and we still adopt the one-boson-exchange (OBE) model, including
the $\pi$, $\sigma$, $\eta$, $\rho$, and $\omega$ exchanges. The corresponding
OBE effective potentials for the $\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}$
systems can be deduced from the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ interactions
by using the heavy quark symmetry and the $SU(3)$ flavor symmetry. Here, we
also consider the coupled channel effect in our calculations. It is very
important in generating the $P_{c}$ and $P_{cs}$ states as hidden-charm
molecular pentaquarks Chen:2019asm ; Chen:2020kco ; Chen:2021tip . By
revisiting the $\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}$ interactions, one can
provide valuable information to search for the possible hidden-charm molecular
pentaquarks with strangeness $|S|=1$. Our comprehensive investigation can also
help us to probe the inner structures or underly mechanism of the $P_{c}$ and
$P_{cs}$ states.

This paper is organized as follows. After introduction, we deduce the OBE
effective potentials in Sec. II. In Sec. III, we present the corresponding
numerical results. The paper ends with a summery in Sec. IV.

##  II OBE effective potentials

According to the chiral symmetry, the effective Lagrangians of the
interactions between light quarks $(u,d,s)$ and the light mesons ($\sigma$,
$\pi$, $\eta$, $\rho$ and $\omega$) are constructed as

$\displaystyle\mathcal{L}_{q}$ $\displaystyle=$ $\displaystyle-
g_{\pi}\bar{\psi}\gamma^{\mu}\gamma^{5}\partial_{\mu}\left(\pi^{i}\tau^{i}+\eta\right)\psi-
g_{\sigma}\bar{\psi}\sigma\psi$ $\displaystyle-
g_{\rho}\bar{\psi}\gamma^{\mu}\left(\rho_{\mu}^{i}\tau^{i}+\omega_{\mu}\right)\psi-
f_{\rho}\bar{\psi}\sigma^{\mu\nu}\partial_{\mu}\left(\rho_{\nu}^{i}\tau^{i}+\omega_{\nu}\right)\psi,$

where $\tau$ is the flavor spatial operator. Therefore, the OBE effective
potentials depend on the spin and isospin of the discussed systems, and have
the form of

$\displaystyle V_{q_{1}q_{2}}$ $\displaystyle=$
$\displaystyle\mathcal{V}_{\sigma,\eta,\omega}(\bm{\sigma}_{1},\bm{\sigma}_{2})+\bm{\tau_{1}}\cdot\bm{\tau_{2}}\mathcal{V}_{\pi,\rho}(\bm{\sigma}_{1},\bm{\sigma}_{2}),$
(2.2)

where $q_{1}$ and $q_{2}$ are the interacting light quarks, respectively.
$\bm{\sigma}_{1}$ and $\bm{\sigma}_{2}$ stand for the spin operators. Here, we
can simply divide the OBE effective potentials into isospin-related part and
isospin-unrelated part, which correspond to the $\pi/\rho$ exchanges and
$\sigma/\eta/\omega$ effective potentials, respectively. For the total
effective potentials, one should sum over the interactions between all the
light quarks in the discussed systems.

Compared to the $\Sigma_{c}^{(*)}\bar{D}^{*}$ systems, the OBE effective
potentials in the isospin-unrelated part for the
$\Xi_{c}^{(^{\prime},*)}\bar{D}^{*}$ systems are very similar as collected in
Table 1. With the help of the $SU(3)_{F}$ symmetry, one can further obtain
relations for the OBE effective potentials in the isospin-related part between
the $\Sigma_{c}^{(*)}\bar{D}^{*}$ and $\Xi_{c}^{(^{\prime},*)}\bar{D}^{*}$
systems. The total isospin for the $\Sigma_{c}^{(*)}\bar{D}^{*}$ systems is
either $I=1/2$ or $I=3/2$, whereas $I=0$ or $I=1$ for the
$\Xi_{c}^{(^{\prime},*)}\bar{D}^{*}$ systems. Here, we firstly expand their
isospin wave functions
$|I_{q_{1}};I_{q_{2}},I_{q_{3}}(I_{q_{2}q_{3}});I\rangle$ in terms of the
$|I_{q_{1}},I_{q_{2}}(I_{q_{1}q_{2}});I_{q_{3}};I\rangle$ basis,

$\displaystyle|\bar{D}^{(*)}\Sigma_{c}^{(*)}(I=\frac{1}{2})\rangle$
$\displaystyle=$
$\displaystyle\left|\frac{1}{2};\frac{1}{2},\frac{1}{2}(1);\frac{1}{2}\right\rangle$
$\displaystyle=$
$\displaystyle\frac{\sqrt{3}}{2}\left|\frac{1}{2},\frac{1}{2}(0);\frac{1}{2};\frac{1}{2}\right\rangle+\frac{1}{2}\left|\frac{1}{2},\frac{1}{2}(1);\frac{1}{2};\frac{1}{2}\right\rangle,$
$\displaystyle|\bar{D}^{(*)}\Sigma_{c}^{(*)}(I=\frac{3}{2})\rangle$
$\displaystyle=$
$\displaystyle\left|\frac{1}{2};\frac{1}{2},\frac{1}{2}(1);\frac{3}{2}\right\rangle=\left|\frac{1}{2},\frac{1}{2}(1);\frac{1}{2};\frac{3}{2}\right\rangle,$
(2.4) $\displaystyle|\bar{D}^{(*)}\Xi_{c}^{(^{\prime},*)}(I)\rangle$
$\displaystyle=$
$\displaystyle\left|\frac{1}{2};\frac{1}{2},0(\frac{1}{2});I\right\rangle=\left|\frac{1}{2},\frac{1}{2}(0);0;I\right\rangle,$
(2.5)

here, we use

$\displaystyle|I_{q_{1}};I_{q_{2}}I_{q_{3}}(I_{q_{2}q_{3}});I\rangle$
$\displaystyle=$
$\displaystyle\sum_{I_{12}}(-1)^{I_{1}+I_{2}+I_{3}+I}\sqrt{(2I_{12}+1)(2I_{23}+1)}$
(2.8) $\displaystyle\times\left\\{\begin{array}[]{ccc}I_{1}&I_{2}&I_{12}\\\
I_{3}&I&I_{23}\end{array}\right\\}|I_{1}I_{2}(I_{12});I_{3};I\rangle.$

Once sandwiching the isospin operator $\bm{\tau_{1}}\cdot\bm{\tau_{2}}$ by the
$|I_{q_{1}},I_{q_{2}}(I_{q_{1}q_{2}});I_{q_{3}};I\rangle$ basis, we can obtain
a serial of relations of the OBE effective potentials in the isospin-related
part, i.e.,

$\displaystyle\mathcal{V}_{\Sigma_{c}^{(*)}\bar{D}^{(*)}}^{I=1/2}$
$\displaystyle=$
$\displaystyle\frac{3}{2}\mathcal{U}^{0}+\frac{1}{2}\mathcal{U}^{1},$ (2.9)
$\displaystyle\mathcal{V}_{\Sigma_{c}^{(*)}\bar{D}^{(*)}}^{I=3/2}$
$\displaystyle=$ $\displaystyle 2\mathcal{U}^{1},$ (2.10)
$\displaystyle\mathcal{V}_{\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}}^{I=0}$
$\displaystyle=$ $\displaystyle\mathcal{U}^{0},$ (2.11)
$\displaystyle\mathcal{V}_{\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}}^{I=1}$
$\displaystyle=$ $\displaystyle\mathcal{U}^{1},$ (2.12)

with

$\displaystyle\mathcal{U}^{0}$ $\displaystyle=$ $\displaystyle\left\langle
I^{\prime}_{q_{1}q_{2}(q_{3})}=0|\bm{\tau_{1}}\cdot\bm{\tau_{2}}|I_{q_{1}q_{2}(q_{3})}=0\right\rangle,$
(2.13) $\displaystyle\mathcal{U}^{1}$ $\displaystyle=$
$\displaystyle\left\langle
I^{\prime}_{q_{1}q_{2}(q_{3})}=1|\bm{\tau_{1}}\cdot\bm{\tau_{2}}|I_{q_{1}q_{2}(q_{3})}=1\right\rangle.$
(2.14)

Table 1: The one-boson-exchange effective potentials for the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ and $\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}$ systems. Here, $\sigma=(u\bar{u}+d\bar{d}+s\bar{s})/\sqrt{3}$, $\eta=(u\bar{u}+d\bar{d}-2s\bar{s})/\sqrt{6}$, $\omega=(u\bar{u}+d\bar{d})\sqrt{2}$, $\pi^{-}(\rho^{-})=\bar{u}d$, $\pi^{0}(\rho^{0})=(u\bar{u}-d\bar{d})/\sqrt{2}$, and $\pi^{+}(\rho^{+})=u\bar{d}$. | $\sigma$ | $\eta$ | $\omega$ | $\pi$ | $\rho$
---|---|---|---|---|---
$\Sigma_{c}^{(*)}\bar{D}^{(*)}[I=1/2]$ | $\mathcal{U}_{\sigma}$ | $\mathcal{U}_{\eta}$ | $2\mathcal{U}_{\omega}$ | $\frac{3}{2}\mathcal{U}_{\pi}^{0}+\frac{1}{2}\mathcal{U}_{\pi}^{1}$ | $\frac{3}{2}\mathcal{U}_{\rho}^{0}+\frac{1}{2}\mathcal{U}_{\rho}^{1}$
$\Sigma_{c}^{(*)}\bar{D}^{(*)}[I=3/2]$ | $\mathcal{U}_{\sigma}$ | $\mathcal{U}_{\eta}$ | $2\mathcal{U}_{\omega}$ | $2\mathcal{U}_{\pi}^{1}$ | $2\mathcal{U}_{\rho}^{1}$
$\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}[I=0]$ | $\mathcal{U}_{\sigma}$ | $-\frac{1}{2}\mathcal{U}_{\eta}$ | $\mathcal{U}_{\omega}$ | $\mathcal{U}_{\pi}^{0}$ | $\mathcal{U}_{\rho}^{0}$
$\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}[I=1]$ | $\mathcal{U}_{\sigma}$ | $-\frac{1}{2}\mathcal{U}_{\eta}$ | $\mathcal{U}_{\omega}$ | $\mathcal{U}_{\pi}^{1}$ | $\mathcal{U}_{\rho}^{1}$

In Table 1, we summarize the OBE effective potentials for the
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$ and $\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}$
systems. In Ref. Chen:2019asm , we have already prepared the concrete OBE
effective potentials for the coupled $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ systems,
and find the interactions from the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ systems
with $I=1/2$ are strong attractive, for the isoquartet systems, the OBE model
provides the weak attractive or repulsive interactions. According to the
relations in Table 1, we can give a qualitative conclusion that the $S-$wave
isoscalar $\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}$ systems may be the possible
strange hidden-charm molecular candidates.

In the heavy quark symmetry, charmed baryons are divided into two multiplets
according to the $SU(3)$ flavor symmetry of light quark cluster, $3_{F}\otimes
3_{F}=\bar{3}_{F}\oplus 6_{F}$, the $\Lambda_{c}$ and $\Xi_{c}$ are in the
$\bar{3}_{F}$ multiplet, whereas, the $\Sigma_{c}^{(*)}$ and
$\Xi_{c}^{(\prime,*)}$ are in the $6_{F}$ multiplet. Thus, the
$\Xi_{c}\bar{D}^{(*)}$ systems aren’t related to the
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$ systems in terms of the $SU(3)$ flavor
symmetry and heavy quark symmetry. In the following, we will deduce the OBE
effective potentials for the $\Xi_{c}\bar{D}^{(*)}$ systems. The general
procedures of the derivations of the OBE effective potentials include three
steps. After constructed the effective Lagrangians, one can firstly write down
the scattering amplitudes $\mathcal{M}_{\text{OBE}}(h_{1}h_{2}\to h_{3}h_{4})$
for the discussed processes. Then, the OBE effective potentials
$\mathcal{V}(\bm{q})$ can be related to the corresponding scattering
amplitudes by using the Breit approximation,
$\mathcal{V}(\bm{q})=-\mathcal{M}(h_{1}h_{2}\to
h_{3}h_{4})/\sqrt{\prod_{i}2M_{i}\prod_{f}2M_{f}}$ with $M_{i}$ and $M_{f}$
being the masses of the initial states ($h_{1}$, $h_{2}$) and final states
($h_{3}$, $h_{4}$), respectively. At last, one can obtain the OBE effective
potentials in the coordinate space $\mathcal{V}(\bm{r})$ after performing the
Fourier transformation, i.e.,

$\displaystyle\mathcal{V}_{E}^{h_{1}h_{2}\to
h_{3}h_{4}}(\bm{r})=\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}e^{i\bm{q}\cdot\bm{r}}\mathcal{V}_{E}^{h_{1}h_{2}\to
h_{3}h_{4}}(\bm{q})\mathcal{F}^{2}(q^{2},m_{E}^{2}).$

Here, for compensating the off-shell effect of the exchanged bosons, a
monopole form factor
$\mathcal{F}(q^{2},m_{E}^{2})=(\Lambda^{2}-m_{E}^{2})/(\Lambda^{2}-q^{2})$ is
introduced at every interactive vertex, where $m_{E}$ and $q$ are the mass and
four-momentum of the exchanged meson, respectively. $\Lambda$ is the cutoff. A
reasonable cutoff value is around 1.00 GeV Tornqvist:1993ng ; Tornqvist:1993vu
. Especially, we also take this empirical value in reproducing the masses of
the $P_{c}(4312)$, $P_{c}(4440)$, and $P_{c}(4457)$Chen:2019asm .

The relevant effective Lagrangians are constructed in the heavy quark limit
and chiral symmetry Yan:1992gz ; Wise:1992hn ; Burdman:1992gh ;
Casalbuoni:1996pg ; Falk:1992cx ; Liu:2011xc , i.e.,

$\displaystyle\mathcal{L}_{H}$ $\displaystyle=$ $\displaystyle
g_{S}\langle\bar{H}^{(\overline{Q})}_{a}\sigma
H^{(\overline{Q})}_{b}\rangle+ig\langle\bar{H}^{(\overline{Q})}_{a}\gamma_{\mu}{\mathcal{A}}_{ab}^{\mu}\gamma_{5}H^{(\overline{Q})}_{b}\rangle$
(2.15)
$\displaystyle-i\beta\langle\bar{H}^{(\overline{Q})}_{a}v_{\mu}\left(\mathcal{V}^{\mu}-\rho^{\mu}\right)_{ab}H^{(\overline{Q})}_{b}\rangle$
$\displaystyle+i\lambda\langle\bar{H}^{(\overline{Q})}_{a}\sigma_{\mu\nu}F^{\mu\nu}(\rho)H^{(\overline{Q})}_{b}\rangle,$
$\displaystyle\mathcal{L}_{\mathcal{B}_{\bar{3}}}$ $\displaystyle=$
$\displaystyle
l_{B}\langle\bar{\mathcal{B}}_{\bar{3}}\sigma\mathcal{B}_{\bar{3}}\rangle+i\beta_{B}\langle\bar{\mathcal{B}}_{\bar{3}}v^{\mu}(\mathcal{V}_{\mu}-\rho_{\mu})\mathcal{B}_{\bar{3}}\rangle,$
$\displaystyle\mathcal{L}_{\mathcal{B}^{(^{\prime},*)}_{6}}$ $\displaystyle=$
$\displaystyle
l_{S}\langle\bar{\mathcal{S}}_{\mu}\sigma\mathcal{S}^{\mu}\rangle-\frac{3}{2}g_{1}\varepsilon^{\mu\nu\lambda\kappa}v_{\kappa}\langle\bar{\mathcal{S}}_{\mu}{\mathcal{A}}_{\nu}\mathcal{S}_{\lambda}\rangle$
$\displaystyle+i\beta_{S}\langle\bar{\mathcal{S}}_{\mu}v_{\alpha}\left(\mathcal{V}^{\alpha}-\rho^{\alpha}\right)\mathcal{S}^{\mu}\rangle+\lambda_{S}\langle\bar{\mathcal{S}}_{\mu}F^{\mu\nu}(\rho)\mathcal{S}_{\nu}\rangle,$
$\displaystyle\mathcal{L}_{\mathcal{B}_{\bar{3}}\mathcal{B}^{(^{\prime},*)}_{6}}$
$\displaystyle=$ $\displaystyle
ig_{4}\langle\bar{\mathcal{S}^{\mu}}{\mathcal{A}}_{\mu}\mathcal{B}_{\bar{3}}\rangle+i\lambda_{I}\varepsilon^{\mu\nu\lambda\kappa}v_{\mu}\langle\bar{\mathcal{S}}_{\nu}F_{\lambda\kappa}\mathcal{B}_{\bar{3}}\rangle+h.c..$

In Eqs. (2.15)-(LABEL:lag2), the multiplet fields $H^{(\bar{Q})}$ and
$\mathcal{S}$ are linear combinations of the $S-$wave charmed mesons and
charmed baryons in the $6_{F}$ flavor representation, respectively.
$H^{(\bar{Q})}=[\tilde{\mathcal{P}}^{*\mu}\gamma_{\mu}-\tilde{\mathcal{P}}\gamma_{5}]\frac{1-\hbox
to0.0pt{/\hss}v}{2}$ with
$\tilde{\mathcal{P}}=\left(\bar{D}^{0},\,D^{-}\right)^{T}$ and
$\tilde{\mathcal{P}}^{*}=\left(\bar{D}^{*0},\,D^{*-}\right)^{T}$.
$\mathcal{S}_{\mu}=-\sqrt{\frac{1}{3}}(\gamma_{\mu}+v_{\mu})\gamma^{5}\mathcal{B}_{6}^{(\prime)}+\mathcal{B}_{6\mu}^{*}$.
$A_{\mu}$ and $\mathcal{V}_{\mu}$ correspond to the axial current and vector
current, respectively,
$A_{\mu}=\frac{1}{2}(\xi^{{\dagger}}\partial_{\mu}\xi-\xi\partial_{\mu}\xi^{{\dagger}})=\frac{i}{f_{\pi}}\partial_{\mu}\mathbb{P}+\ldots$,
and
$\mathcal{V}_{\mu}=\frac{1}{2}(\xi^{{\dagger}}\partial_{\mu}\xi-\xi\partial_{\mu}\xi^{{\dagger}})=\frac{i}{2f_{\pi}^{2}}\left[\mathbb{P},\partial_{\mu}\mathbb{P}\right]+\ldots$
with $\xi=\text{exp}(i\mathbb{P}/f_{\pi})$ and the pion decay constant
$f_{\pi}=132$ MeV. $\rho_{ba}^{\mu}=ig_{V}\mathbb{V}_{ba}^{\mu}/\sqrt{2}$,
$F^{\mu\nu}(\rho)=\partial^{\mu}\rho^{\nu}-\partial^{\nu}\rho^{\mu}+\left[\rho^{\mu},\rho^{\nu}\right]$.
$\mathbb{P}$ and $\mathbb{V}$ stand for the isoscalar and vector matrixes,
respectively. Matrices for the $\mathcal{B}_{\bar{3}}$,
$\mathcal{B}_{6}^{(^{\prime},*)}$, $\mathbb{P}$, and $\mathbb{V}$ are
expressed as

$\displaystyle\mathcal{B}_{\bar{3}}$ $\displaystyle=$
$\displaystyle{\left(\begin{array}[]{ccc}0&\Lambda_{c}^{+}&\Xi_{c}^{+}\\\
-\Lambda_{c}^{+}&0&\Xi_{c}^{0}\\\
-\Xi_{c}^{+}&-\Xi_{c}^{0}&0\end{array}\right),}$
$\displaystyle\mathcal{B}_{6}^{(^{\prime},*)}$ $\displaystyle=$
$\displaystyle{\left(\begin{array}[]{ccc}\Sigma_{c}^{(*)++}&\frac{1}{\sqrt{2}}\Sigma_{c}^{(*)+}&\frac{1}{\sqrt{2}}\Xi_{c}^{(^{\prime},*)+}\\\
\frac{1}{\sqrt{2}}\Sigma_{c}^{(*)+}&\Sigma_{c}^{(*)0}&\frac{1}{\sqrt{2}}\Xi_{c}^{(^{\prime},*)0}\\\
\frac{1}{\sqrt{2}}\Xi_{c}^{(^{\prime},*)+}&\frac{1}{\sqrt{2}}\Xi_{c}^{(^{\prime},*)0}&\Omega_{c}^{(*)0}\\\
\end{array}\right)},$ $\displaystyle\mathbb{P}$ $\displaystyle=$
$\displaystyle{\left(\begin{array}[]{ccc}\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}&\pi^{+}&K^{+}\\\
\pi^{-}&-\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}&K^{0}\\\
K^{-}&\bar{K}^{0}&-\frac{2}{\sqrt{6}}\eta\end{array}\right)},$
$\displaystyle{V}$ $\displaystyle=$
$\displaystyle{\left(\begin{array}[]{ccc}\frac{\rho^{0}}{\sqrt{2}}+\frac{\omega}{\sqrt{2}}&\rho^{+}&K^{*+}\\\
\rho^{-}&-\frac{\rho^{0}}{\sqrt{2}}+\frac{\omega}{\sqrt{2}}&K^{*0}\\\
K^{*-}&\bar{K}^{*0}&\phi\end{array}\right)},$

respectively.

In this work, we adopt the same values of the coupling constants in Ref.
Liu:2011xc ; Chen:2019asm ; Chen:2020kco , where $g=0.59$, $g_{1}=0.94$, and
$g_{4}=1.06$ are extracted from the decay widths of $\Gamma(D^{*}\to D\pi)$
and $\Gamma(\Sigma_{c}^{(*)}\to\Lambda_{c}\pi)$ Isola:2003fh ; Liu:2011xc ;
pdg , respectively. In the charmed mesons sector, $g_{S}=0.76$ Bardeen:2003kt
, $\beta=0.9$ Isola:2003fh , $\lambda=0.56$ GeV-1 Isola:2003fh ,
$g_{V}=m_{\rho}/f_{\pi}=5.9$. For the remaining coupling constants between the
heavy baryons and the light meson ($\sigma$, $\rho$, $\omega$), their values
are estimated from the nucleon-nucleon interactions Liu:2011xc ,
$l_{B}=-3.65$, $\beta_{B}g_{V}=6.0$, $l_{S}=6.2$, $\beta_{S}g_{V}=-12.0$,
$\lambda_{S}g_{V}=-19.2$ $\text{GeV}^{-1}$, $\lambda_{I}g_{V}=6.8$
$\text{GeV}^{-1}$. With these preparations, we can deduce the detailed
expressions of the OBE effective potentials, i.e.,

$\displaystyle\mathcal{V}^{\Xi_{c}\bar{D}\to\Xi_{c}\bar{D}}(r)$
$\displaystyle=$ $\displaystyle
2\mathbb{A}Y_{\Lambda,m_{\sigma}}+\frac{\mathbb{B}}{4}\left(\mathcal{G}(I)Y_{\Lambda,m_{\rho}}+Y_{\Lambda,m_{\omega}}\right),\quad$
(2.21)
$\displaystyle\mathcal{V}^{\Xi_{c}^{\prime}\bar{D}\to\Xi_{c}\bar{D}}(r)$
$\displaystyle=$ $\displaystyle 0,$ (2.22)
$\displaystyle\mathcal{V}^{\Xi_{c}\bar{D}^{*}\to\Xi_{c}\bar{D}}(r)$
$\displaystyle=$ $\displaystyle 0,$ (2.23)
$\displaystyle\mathcal{V}^{\Xi_{c}^{*}\bar{D}\to\Xi_{c}\bar{D}}(r)$
$\displaystyle=$ $\displaystyle 0,$ (2.24)
$\displaystyle\mathcal{V}^{\Xi_{c}^{\prime}\bar{D}^{*}\to\Xi_{c}\bar{D}}(r)$
$\displaystyle=$
$\displaystyle\frac{\mathbb{C}}{6\sqrt{6}}\left(\mathcal{G}(I)\mathcal{Z}^{13}_{\Lambda_{1},m_{\pi
1}}+\mathcal{Z}^{13}_{\Lambda_{1},m_{\eta 1}}\right)$ (2.25)
$\displaystyle-\frac{\mathbb{D}}{3\sqrt{6}}\left(\mathcal{G}(I)\mathcal{Z}^{\prime
13}_{\Lambda_{1},m_{\rho 1}}+\mathcal{Z}^{\prime 13}_{\Lambda_{1},m_{\omega
1}}\right),$
$\displaystyle\mathcal{V}^{\Xi_{c}^{*}\bar{D}^{*}\to\Xi_{c}\bar{D}}(r)$
$\displaystyle=$
$\displaystyle-\frac{\mathbb{C}}{6\sqrt{2}}\left(\mathcal{G}(I)\mathcal{Z}^{14}_{\Lambda_{2},m_{\pi
2}}+\mathcal{Z}^{14}_{\Lambda_{2},m_{\eta 2}}\right)$ (2.26)
$\displaystyle-\frac{\mathbb{D}}{6\sqrt{2}}\left(\mathcal{G}(I)\mathcal{Z}^{\prime
14}_{\Lambda_{2},m_{\rho 2}}+\mathcal{Z}^{\prime 14}_{\Lambda_{2},m_{\omega
2}}\right),$

for the $\Xi_{c}^{(\prime,*)}\bar{D}^{(*)}\to\Xi_{c}\bar{D}$ processes, where
$\mathbb{A}=l_{B}g_{s}$, $\mathbb{B}=\beta\beta_{B}g_{v}^{2}$,
$\mathbb{C}=gg_{4}/f_{\pi}^{2}$, $\mathbb{D}=\lambda\lambda_{I}g_{v}^{2}$,
$\mathcal{G}(I)$ is the isospin factor, which is taken as $1$ for the
isospin-$1$ system, and $-3$ for the isospin-$0$ system. The functions
$Y_{\Lambda,m}$, $\mathcal{Z}^{ij}_{\Lambda,m_{a}}$, and $\mathcal{Z}^{\prime
ij}_{\Lambda,m_{a}}$ denote

$\displaystyle Y_{\Lambda,m}$ $\displaystyle=$ $\displaystyle\frac{1}{4\pi
r}(e^{-mr}-e^{-\Lambda r})-\frac{\Lambda^{2}-m^{2}}{8\pi\Lambda}e^{-\Lambda
r},$ (2.27) $\displaystyle\mathcal{Z}^{ij}_{\Lambda,m_{a}}$ $\displaystyle=$
$\displaystyle\left(\mathcal{E}_{ij}\nabla^{2}+\mathcal{F}_{ij}r\frac{\partial}{\partial
r}\frac{1}{r}\frac{\partial}{\partial r}\right)Y_{\Lambda,m_{a}},$ (2.28)
$\displaystyle\mathcal{Z}^{\prime ij}_{\Lambda,m_{a}}$ $\displaystyle=$
$\displaystyle\left(2\mathcal{E}_{ij}\nabla^{2}-\mathcal{F}_{ij}r\frac{\partial}{\partial
r}\frac{1}{r}\frac{\partial}{\partial r}\right)Y_{\Lambda,m_{a}},$ (2.29)

The variables in Eqs. (2.21)-(2.26) are defined as
$\Lambda_{i}^{2}=\Lambda^{2}-q_{i}^{2}$, $m_{{i}}^{2}=m^{2}-q_{i}^{2}$, with
$i=1,2$. $q_{1}=0.64$ MeV and $q_{2}=38.14$ MeV. The spin-spin interaction and
tensor force operators read as
$\mathcal{E}_{13}=\vec{\sigma}\cdot\vec{\epsilon}_{4}^{{\dagger}}$,
$\mathcal{E}_{14}=\sum_{1/2,a;1,b}^{3/2,a+b}\chi_{3,a}^{{\dagger}}\vec{\epsilon}_{3,b}^{{\dagger}}\cdot\vec{\epsilon}_{4}^{{\dagger}}\chi_{1}$,
$\mathcal{F}_{13}=S\left(\hat{r},\vec{\sigma},\vec{\epsilon}_{4}^{{\dagger}}\right)$,
and
$\mathcal{F}_{14}=\sum_{1/2,a;1,b}^{3/2,a+b}\chi_{3,a}^{{\dagger}}S\left(\hat{r},\vec{\epsilon}_{3,b}^{{\dagger}},\vec{\epsilon}_{4}^{{\dagger}}\right)\chi_{1}$.
In addition, one can refer to the concrete subpotentials for the
$\Xi_{c}^{(\prime,*)}\bar{D}^{(*)}\to\Xi_{c}\bar{D}^{*}$ processes in Ref.
Chen:2020kco .

##  III Numerical results

Before producing numerical calculations, we would like to make several remarks
on the bound state solutions for the reasonable loosely bound molecular state
after performing the coupled channel analysis Chen:2017xat ,

  1. 1.

The reasonable cutoff value in the monopole form factor is around 1.00 GeV
according to the experience of the nucleon-nucleon interaction
Tornqvist:1993ng ; Tornqvist:1993vu .

  2. 2.

The binding energy is around several to a few tens MeV.

  3. 3.

The root-mean-square (RMS) radii are around a few fm or larger as the size of
the loosely bound molecule should be much larger than the size of the
component hadrons.

  4. 4.

For an $S-$wave molecular state composed by two mesons, the asymptotic form of
its wave function can be expressed as $\psi(r)\sim e^{-\sqrt{2\mu E}r}/r$,
where $\mu$ and $E$ stand for the reduced mass and the binding energy,
respectively. In the coupled channel analysis, the $\mu$ is the reduced mass
for the dominant channel, and $E$ is measured from the dominant channel,
$E=M_{\text{lowest}}-M_{\text{dominant}}+E_{\text{binding}}$. When we use the
approximated wave function, we can obtain the relation between the molecular
size and its binding energy, $R\sim 1/\sqrt{2\mu E}$ Chen:2017xat ;
Close:2003sg . By using this relation, the system with the lowest mass
threshold is the dominant channel.

With the above preparations, we firstly solve the coupled channel
Shr$\ddot{\text{o}}$dinger equation to find the bound state solutions
(including the binding energies $E$, the RMS radii $r_{RMS}$, and the
probabilities for all the discussed channels $p_{i}$) for the $S-$wave coupled
$\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}$ systems with all the possible quantum
numbers. The cutoff value is taken from 0.8 GeV to 5.0 GeV. And then, we check
whether the obtained bound state solutions match the above criterion of the
reasonable loosely bound molecular states (the binding energy around a few to
10 MeV and the RMS radius around or larger than 1 fm). Finally, we compare
their cutoff value in these reasonable loosely bound states. In our analysis,
we conclude that the state with the obtained reasonable loosely bound state
solutions in the cutoff range $\Lambda$ around 1.00 GeV is the prime loosely
bound molecular candidates. With the increasing of the cutoff value, the
possible of the existence of the loosely bound molecule becomes lower. For the
reasonable loosely bound states with the cutoff $1.00<\Lambda<2.00$ GeV, they
may be the possible molecular candidate. If the cutoff $\Lambda$ is larger
than 3.00 GeV, we conclude that they cannot be a good molecular candidate.

Figure 1: The $\Lambda$ dependence of the binding energy $E$, root-mean-square
radii, and the probabilities for all discussed channels. The blue solid lines
and the red dashed lines in the left figures correspond to the binding
energies curves and root-mean-square radii curves, respectively.

As shown in Figure 1, we present the $\Lambda$ dependence of the bound state
properties ($E$, $r_{RMS}$, and $p_{i}$) for the coupled
$\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
states with $0,1(1/2^{-})$, the coupled
$\Xi_{c}^{*}\bar{D}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ states
with $0,1(3/2^{-})$, and the coupled
$\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ states with $0(1/2^{-})$
and $0,1(3/2^{-})$. For the coupled
$\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
state with $0(1/2^{-})$, when the cutoff is taken around 0.8 GeV, the binding
energy is around $-10$ MeV, the RMS radius is around 1.0 fm, and the dominant
channel is the $S-$wave $\Xi_{c}^{\prime}\bar{D}$ with its probability around
70%. If we adopt the former remarks, these bound state solutions are
consistent with the reasonable loosely bound state properties. Therefore, we
can conclude that this state can be the possible strange hidden-charm
molecular pentaquark. Since the probabilities for the remaining channels are
almost 30 percent, the coupled channel effects play an important role to
generate this coupled bound state. For the isovector coupled
$\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
state with $1/2^{-}$, the binding energy appears at cutoff $\Lambda$ around
1.10 GeV, the RMS radius is over or around 1.0 fm, and the $S-$wave
$\Xi_{c}^{\prime}\bar{D}$ component is dominant, the corresponding probability
is around 90 %. The reasonable cutoff, the reasonable RMS radius and the small
binding energy indicate that the oupled
$\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
state with $1(1/2^{-})$ can be the other possible strange hidden-charm
molecular candidate.

For the coupled
$\Xi_{c}^{*}\bar{D}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ states
with $0,1(3/2^{-})$, and the coupled
$\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ states with
$0,1(3/2^{-})$, their bound states solutions are similar to the corresponding
results for the coupled
$\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
states with $0,1(1/2^{-})$. In the range of $1.00<\Lambda<2.00$ GeV, we can
obtain the loosely binding energies and the reasonable RMS radii. And the
dominant channels are the systems with lowest mass threshold in these
discussed coupled systems. Meanwhile, we also find the interactions from the
isoscalar states are stronger than those in the isovector states as the
corresponding cutoff values are a little smaller than those in the isovector
bound states with the same binding energy. If we still adopt the criterions of
the loosely bound molecular states, the coupled
$\Xi_{c}^{*}\bar{D}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ states
with $0,1(3/2^{-})$, and the coupled
$\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ states with
$0,1(3/2^{-})$ can be the possible strange hidden-charm molecular candidates.

As shown in Figure 1, we also obtain the loosely bound state solutions for the
coupled $\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ states with
$0(1/2^{-})$ when the cutoff $\Lambda$ is taken around 1.00 GeV. The dominant
channel is the $S-$wave $\Xi_{c}^{\prime}\bar{D}^{*}$ component. Since its
probability is over 90 % in the cutoff range, the coupled channel effects play
a minor role in forming this bound state. Therefore, this bound state can be
recommended as the possible strange hidden-charm molecular pentaquark.

Figure 2: The $\Lambda$ dependence of the binding energy $E$ and RMS radius
for the single $\Xi_{c}^{*}\bar{D}^{*}$ systems with
$I(J^{P})=0(1/2^{-},3/2^{-},5/2^{-})$ and $1(1/2^{-},3/2^{-})$. The blue solid
lines and the red dashed lines correspond to the binding energies curves and
root-mean-square radii curves, respectively.

Meanwhile, we also perform a single $\Xi_{c}^{*}\bar{D}^{*}$ channel analysis
to search for possible molecular states. When we vary the cutoff value in the
range of 0.80$\sim$5.00 GeV, we find that

  1. 1.

There don’t exist bound state solutions for the isovector
$\Xi_{c}^{*}\bar{D}^{*}$ system with $5/2^{-}$.

  2. 2.

For the isovector $\Xi_{c}^{*}\bar{D}^{*}$ system with $3/2^{-}$, it cannot be
a good strange hidden-charm molecular candidate as its cutoff $\Lambda$ is
larger than 3.00 GeV.

  3. 3.

For the $\Xi_{c}^{*}\bar{D}^{*}$ with $1(1/2^{-})$ and $0(5/2^{-})$, they may
the possible strange hidden-charm molecular candidates as their binding
energies appear at the cutoff $\Lambda$ around 2.00 GeV.

  4. 4.

If we still take the cutoff value $\Lambda$ around 1.00 GeV as a reasonable
parameter, the isoscalar $\Xi_{c}^{*}\bar{D}^{*}$ molecules with
$J^{P}=(1/2^{-},3/2^{-})$ can be prime strange hidden-charm molecular
candidates.

  5. 5.

As shown in Figure 2, the OBE interactions become stronger with the increasing
of the cutoff value for one bound state. When we adopt this rough property in
the isoscalar single $\Xi_{c}^{*}\bar{D}^{*}$ bound states, one can conclude
that the interactions with the higher spin are a little weaker attractive
according to the cutoff relation
$\Lambda(1/2^{-})<\Lambda(3/2^{-})<\Lambda(5/2^{-})$ with the same binding
energy. This is also consist with predictions in Ref. Wang:2019nvm .

Table 2: The bound state solutions of the coupled $\Xi_{c}^{(\prime,*)}\bar{D}^{(*)}$ systems with $I(J^{P})=0(1/2^{-})$ and $1(1/2^{-})$. The cutoff $\Lambda$, the root-mean-square $r_{RMS}$, and the mass of the bound state $M$ are in the units of GeV, fm, and MeV, respectively. $p_{i}(\%)$ denotes the probability of the $i-$th channel for the investigated system. The largest probability of the quantum number configuration for a bound state are remarked by bold typeface. And $\tilde{E}$, in the unit of MeV, is the mass gap between the bound state mass and the threshold of the dominant channel. $I(J^{P})$ | ${0({1}/{2}^{-})}$ | ${1({1}/{2}^{-})}$ | ${1({1}/{2}^{-})}$
---|---|---|---
$\Lambda$ |  | 0.99 | 1.00 | 1.01 | 1.27 | 1.28 | 1.29 | 2.03 | 2.04 | 2.05
$r_{RMS}$ |  | 1.29 | 0.60 | 0.49 | 0.40 | 0.36 | 0.34 | 0.89 | 0.40 | 0.34
$M$ |  | 4334.64 | 4327.69 | 4319.82 | 4327.53 | 4315.33 | 4302.35 | 4583.26 | 4574.51 | 4565.02
$\tilde{E}$ |  |  $-$109.37 |  $-$116.32 |  $-$124.19 |  $-$116.48 |  $-$128.68 |  $-$141.66 |  $-$71.2 |  $-$79.95 |  $-$89.44
$p_{i}(\%)$ | $\Xi_{c}\bar{D}$ | 30.36 | 16.51 | 12.37 | 2.12 | 1.48 | 1.19 | $\ldots$ | $\ldots$ | $\ldots$
| $\Xi_{c}^{\prime}\bar{D}$ | 32.56 | 38.99 | 40.73 | 62.36 | 62.33 | 62.10 | $\ldots$ | $\ldots$ | $\ldots$
| $\Xi_{c}\bar{D}^{*}$ | 21.75 | 26.38 | 27.89 | 35.30 | 35.97 | 36.50 | $\ldots$ | $\ldots$ | $\ldots$
| $\Xi_{c}^{\prime}\bar{D}^{*}$ | 13.07 | 15.58 | 16.40 | 0.17 | 0.18 | 0.19 | 16.33 | 8.57 | 6.91
| $\Xi_{c}^{*}\bar{D}^{*}$ | 2.26 | 2.54 | 2.61 | 0.04 | 0.03 | 0.02 | 83.67 | 91.43 | 93.09

In our calculations, we also obtain the other kind of bound state solutions,
which is very different with those for a reasonable loosely bound hadronic
molecule. In this case, their RMS radii are around 0.5 fm or more less,
therefore they cannot be reasonable hadronic molecular states but the tightly
bound states. According to the relation $R\sim
1/\sqrt{2\tilde{\mu}\tilde{E}}$, the system with a higher mass is always the
dominant channel for the tightly bound state.

As shown in Table 2, we obtain three tightly bound state solutions, the
coupled
$\Xi_{c}\bar{D}/\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
states with $I(J^{P})=0(1/2^{-})$ and $1(1/2^{-})$ and the coupled
$\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ state with $1(1/2^{-})$,
where the other binding energies
$\tilde{E}=E+(M_{\text{lowest}}-M_{\text{dominant}})$ are around 100 MeV, the
corresponding dominant channels are $\Xi_{c}^{\prime}\bar{D}$ and
$\Xi_{c}^{*}\bar{D}^{*}$, respectively. Compared to the coupled
$\Xi_{c}\bar{D}/\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
states with $I(J^{P})=0(1/2^{-})$ and $1(1/2^{-})$, the cutoff value in the
coupled $\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ state with
$1(1/2^{-})$ is a little away from the empirical value $\Lambda$ around 1.00
GeV.

For the
$\Xi_{c}\bar{D}/\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
interactions, if we recall that such interactions are obtained by adding the
lowest channel from the
$\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
interactions, it is not strange to find the analogous results in the
$\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
interactions. Thus, the
$\Xi_{c}\bar{D}/\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
tightly bound states with $0,1(1/2^{-})$ here are not the independent states,
but correspond to the $\Xi_{c}^{\prime}\bar{D}$ loosely bound molecular
pentaquarks with $0,1(1/2^{-})$.

##  IV conclusion and discussion

As a hot issue in the hadron physics, whether the new hadron states are the
hadronic molecules is still open to discuss. In 2019, the discovery of the
three $P_{c}$ states could provide a strong evidence of the existence of the
hidden-charm molecular pentaquarks Aaij:2019vzc . In Ref. Chen:2019asm , the
$P_{c}(4312)$, $P_{c}(4440)$, and $P_{c}(4457)$ states can be assigned as the
$\Sigma_{c}\bar{D}$ molecular state with $1/2(1/2^{-})$, the
$\Sigma_{c}\bar{D}^{*}$ states with $1/2(1/2^{-})$ and $1/2(3/2^{-})$,
respectively. This is not the end of the story. Very recently, the LHCb
Collaboration reported a new evidence of the strange hidden-charm pentaquark
$P_{cs}(4459)$ in the $\Xi_{b}^{-}\to J/\psi\Lambda K^{-}$ process. Theorists
also proposed that the $P_{cs}(4459)$ can be the strange hidden-charm
molecular state composed of the $S-$wave $\Xi_{c}\bar{D}^{*}$ state
Chen:2020kco ; Chen:2021tip ; Chen:2020uif ; Peng:2020hql ; Wang:2020eep .

Figure 3: A summary of the mass spectrum of the $P_{c}$ states in the meson-
baryon molecular scenario Chen:2019asm  and the predicted possible hidden-
charm molecular pentaquarks candidates with $|S|=1$. Here, we roughly estimate
the mass positions of the predicted hidden-charm molecular pentaquarks
according to the values of the cutoff $\Lambda$. The red and green lines label
the predicted molecular candidates with $I=0$ and $I=1$, respectively.

In this work, we preform the coupled channel analysis on the interactions
between a charm-strange baryon $\Xi_{c}^{(^{\prime},*)}$ and an anti-charmed
meson $\bar{D}^{(*)}$ in the framework of the OBE model. By using the $SU(3)$
flavor symmetry and the heavy quark symmetry, we can obtain the OBE effective
potentials relations between the $\Xi_{c}^{(\prime,*)}\bar{D}^{(*)}$
interactions and the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ interactions, where we
have already deduced the concrete OBE effective potentials in Ref.
Chen:2019asm . For the $\Xi_{c}\bar{D}^{(*)}$ systems, their $SU(3)$ flavor
partners are the $\Lambda_{c}\bar{D}^{(*)}$ systems, which is not prepared in
our former work. In this work, we derive the corresponding OBE effective
potentials in the general procedures.

Our results can predict several possible hidden-charm molecular pentaquarks
with strangeness $|S|=1$. As shown in the Figure 3, there can exist seven
possible isoscalar hidden-charm molecular pentaquarks with $|S|=1$, they are
mainly composed of the $\Xi_{c}^{\prime}\bar{D}$ state with
$I(J^{P})=0(1/2^{-})$, the $\Xi_{c}\bar{D}^{*}$ state with $0(3/2^{-})$, the
$\Xi_{c}^{*}\bar{D}$ state with $0(3/2^{-})$, the
$\Xi_{c}^{\prime}\bar{D}^{*}$ states with $0(1/2^{-},3/2^{-})$, and the
$\Xi_{c}^{*}\bar{D}^{*}$ states with $0(1/2^{-},3/2^{-})$, respectively. The
$\Xi_{c}^{*}\bar{D}^{*}$ state with $0(5/2^{-})$ may be also the possible
strange hidden-charm molecular candidate. Meanwhile, we find there may exist
three isovector hidden-charm molecular pentaquarks with $|S|=1$, like the
$\Xi_{c}^{\prime}\bar{D}$ state with $1(1/2^{-})$, the $\Xi_{c}^{*}\bar{D}$
state with $1(3/2^{-})$, and the $\Xi_{c}^{\prime}\bar{D}^{*}$ state with
$1(3/2^{-})$. Our results also indicate that the coupled channel effect play a
very important role in forming these hidden-charm molecular candidates with
strangeness $|S|=1$, especially for the
$\Xi_{c}^{\prime}\bar{D}/\Xi_{c}\bar{D}^{*}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
coupled states with $I(J^{P})=0,1(1/2^{-})$, the
$\Xi_{c}^{*}\bar{D}/\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$
coupled state with $0(3/2^{-})$, and the
$\Xi_{c}^{\prime}\bar{D}^{*}/\Xi_{c}^{*}\bar{D}^{*}$ coupled states with
$0,1(3/2^{-})$. We are looking forward that the future experiments can search
for possible hidden-charm pentaquarks with strangeness around these predicted
mass thresholds.

Compared to the $P_{c}$ states assigned to the meson-baryon molecules, as
shown in Figure 3, there are four more possible hidden-charm molecular
pentaquarks with strangeness $|S|=1$, which include one isoscalar
$\Xi_{c}\bar{D}^{*}$ state with $3/2^{-}$ and three isovector bound states. In
our previous work Chen:2020kco ; Chen:2020uif , the $P_{cs}(4459)$ can be
explained as the isoscalar $\Xi_{c}\bar{D}^{*}$ state with $3/2^{-}$, and the
$\Xi_{c}^{*}\bar{D}$ channel is also very important. As we seen, the mass
difference between the $\Lambda_{c}\bar{D}^{*}$ and $\Sigma_{c}\bar{D}^{*}$
systems is much larger than that between the $\Xi_{c}\bar{D}^{*}$ and
$\Xi_{c}^{*}\bar{D}$ systems, which can weaken the contribution from the
coupled channel effect, this may explain the reason why there cannot exist the
possible hidden-charm molecular pentaquarks composed of the
$\Lambda_{c}\bar{D}^{*}$ state with $1/2(3/2^{-})$.

As shown in Table 1, the OBE interactions for the isovector
$\Xi_{c}^{(^{\prime},*)}\bar{D}^{(*)}$ systems are weaker repulsive or
stronger attractive than those from the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$
systems with $I=3/2$. This is the main reason that we may predict three
isovector strange hidden-charm molecular pentaquarks instead of the hidden-
charm molecular pentaquarks with $I=3/2$ Chen:2019asm .

## ACKNOWLEDGMENTS

Rui Chen is very grateful to Shi-Lin Zhu for helpful discussions and
constructive suggestions. This work is supported by the China National Funds
for Distinguished Young Scientists under Grant No. 11825503 and by the
National Program for Support of Top-notch Young Professionals. R. C. is
supported by the National Postdoctoral Program for Innovative Talent.

## References

  * (1) R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 115, 072001 (2015)
  * (2) R. Chen, X. Liu, X. Q. Li and S. L. Zhu, Phys. Rev. Lett. 115, 132002 (2015)
  * (3) H. X. Chen, W. Chen, X. Liu, T. G. Steele and S. L. Zhu, Phys. Rev. Lett. 115, 172001 (2015)
  * (4) M. Karliner and J. L. Rosner, Phys. Rev. Lett. 115, no. 12, 122001 (2015)
  * (5) L. Roca, J. Nieves and E. Oset, Phys. Rev. D 92, 094003 (2015)
  * (6) A. Mironov and A. Morozov, JETP Lett. 102, 271 (2015)
  * (7) J. He, Phys. Lett. B 753, 547-551 (2016)
  * (8) U. G. Meissner and J. A. Oller, Phys. Lett. B 751, 59 (2015)
  * (9) T. J. Burns, Eur. Phys. J. A 51, no.11, 152 (2015)
  * (10) Y. Shimizu, D. Suenaga and M. Harada, Phys. Rev. D 93, 114003 (2016)
  * (11) R. Chen, X. Liu and S. L. Zhu, Nucl. Phys. A 954, 406 (2016)
  * (12) M. I. Eides, V. Y. Petrov and M. V. Polyakov, Phys. Rev. D 93, no. 5, 054039 (2016)
  * (13) H. Huang, C. Deng, J. Ping and F. Wang, Eur. Phys. J. C 76, no.11, 624 (2016) [arXiv:1510.04648 [hep-ph]].
  * (14) H. X. Chen, E. L. Cui, W. Chen, X. Liu, T. G. Steele and S. L. Zhu, Eur. Phys. J. C 76, no.10, 572 (2016)
  * (15) G. Yang and J. Ping, Phys. Rev. D 95, no.1, 014010 (2017)
  * (16) J. He, Phys. Rev. D 95, no.7, 074004 (2017)
  * (17) Y. Yamaguchi and E. Santopinto, Phys. Rev. D 96, no.1, 014018 (2017)
  * (18) L. Maiani, A. D. Polosa and V. Riquer, Phys. Lett. B 749, 289 (2015)
  * (19) G. N. Li, X. G. He and M. He, JHEP 12, 128 (2015)
  * (20) R. Ghosh, A. Bhattacharya and B. Chakrabarti, Phys. Part. Nucl. Lett. 14, no.4, 550-552 (2017)
  * (21) Z. G. Wang and T. Huang, Eur. Phys. J. C 76, no.1, 43 (2016)
  * (22) V. V. Anisovich, M. A. Matveev, J. Nyiri, A. V. Sarantsev and A. N. Semenova, Int. J. Mod. Phys. A 30, 1550190 (2015)
  * (23) R. F. Lebed, Phys. Lett. B 749, 454 (2015) [arXiv:1507.05867 [hep-ph]].
  * (24) R. Zhu and C. F. Qiao, Phys. Lett. B 756, 259-264 (2016) [arXiv:1510.08693 [hep-ph]].
  * (25) F. K. Guo, U. G. Meissner, W. Wang and Z. Yang, Phys. Rev. D 92, 071502 (2015)
  * (26) X. H. Liu, Q. Wang and Q. Zhao, Phys. Lett. B 757, 231-236 (2016)
  * (27) M. Mikhasenko, arXiv:1507.06552 [hep-ph].
  * (28) H. X. Chen, W. Chen, X. Liu and S. L. Zhu, Phys. Rept. 639, 1 (2016)
  * (29) Y. R. Liu, H. X. Chen, W. Chen, X. Liu and S. L. Zhu, Prog. Part. Nucl. Phys. 107, 237-320 (2019)
  * (30) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan, Phys. Rept. 873, 1-154 (2020)
  * (31) 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)
  * (32) A. Esposito, A. Pilloni and A. D. Polosa, Phys. Rept. 668, 1-97 (2017)
  * (33) A. Hosaka, T. Iijima, K. Miyabayashi, Y. Sakai and S. Yasui, PTEP 2016, no.6, 062C01 (2016)
  * (34) Z. C. Yang, Z. F. Sun, J. He, X. Liu and S. L. Zhu, Chin. Phys. C 36, 6 (2012)
  * (35) J. J. Wu, R. Molina, E. Oset and B. S. Zou, Phys. Rev. Lett. 105, 232001 (2010)
  * (36) W. L. Wang, F. Huang, Z. Y. Zhang and B. S. Zou, Phys. Rev. C 84, 015203 (2011)
  * (37) X. Q. Li and X. Liu, Eur. Phys. J. C 74, no.12, 3198 (2014)
  * (38) J. J. Wu, R. Molina, E. Oset and B. S. Zou, Phys. Rev. C 84, 015202 (2011)
  * (39) R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 122, no. 22, 222001 (2019)
  * (40) R. Chen, Z. F. Sun, X. Liu and S. L. Zhu, Phys. Rev. D 100, no. 1, 011502 (2019)
  * (41) M. Z. Liu, Y. W. Pan, F. Z. Peng, M. Sanchez Sanchez, L. S. Geng, A. Hosaka and M. Pavon Valderrama, Phys. Rev. Lett. 122, no. 24, 242001 (2019)
  * (42) J. He, Eur. Phys. J. C 79, no. 5, 393 (2019)
  * (43) L. Meng, B. Wang, G. J. Wang and S. L. Zhu, Phys. Rev. D 100, no. 1, 014031 (2019)
  * (44) T. J. Burns and E. S. Swanson, Phys. Rev. D 100, no.11, 114033 (2019)
  * (45) J. Ferretti, E. Santopinto, M. Naeem Anwar and M. A. Bedolla, Phys. Lett. B 789, 562-567 (2019)
  * (46) X. Z. Weng, X. L. Chen, W. Z. Deng and S. L. Zhu, Phys. Rev. D 100, no.1, 016014 (2019)
  * (47) A. Ali and A. Y. Parkhomenko, Phys. Lett. B 793, 365-371 (2019)
  * (48) J. F. Giron, R. F. Lebed and C. T. Peterson, JHEP 05, 061 (2019)
  * (49) J. B. Cheng and Y. R. Liu, Phys. Rev. D 100, no.5, 054002 (2019)
  * (50) H. Mutuk, Chin. Phys. C 43, no.9, 093103 (2019)
  * (51) R. Chen, J. He and X. Liu, Chin. Phys. C 41, no.10, 103105 (2017)
  * (52) B. Wang, L. Meng and S. L. Zhu, Phys. Rev. D 101, no.3, 034018 (2020)
  * (53) Z. G. Wang, Eur. Phys. J. C 76, no.3, 142 (2016)
  * (54) A. Feijoo, V. K. Magas, A. Ramos and E. Oset, Eur. Phys. J. C 76, no.8, 446 (2016)
  * (55) J. X. Lu, E. Wang, J. J. Xie, L. S. Geng and E. Oset, Phys. Rev. D 93, 094009 (2016)
  * (56) C. W. Xiao, J. Nieves and E. Oset, Phys. Lett. B 799, 135051 (2019)
  * (57) Q. Zhang, B. R. He and J. L. Ping, [arXiv:2006.01042 [hep-ph]].
  * (58) C. W. Shen, H. J. Jing, F. K. Guo and J. J. Wu, Symmetry 12, no.10, 1611 (2020)
  * (59) J. Ferretti and E. Santopinto, JHEP 04, 119 (2020)
  * (60) J. T. Zhu, L. Q. Song and J. He, Phys. Rev. D 103, no.7, 074007 (2021)
  * (61) C. W. Xiao, J. J. Wu and B. S. Zou, Phys. Rev. D 103, no.5, 054016 (2021)
  * (62) F. L. Wang, R. Chen and X. Liu, Phys. Rev. D 103, no.3, 034014 (2021)
  * (63) F. L. Wang, X. D. Yang, R. Chen and X. Liu, Phys. Rev. D 103, no.5, 054025 (2021)
  * (64) R. Aaij et al. [LHCb], Sci. Bull. 66, 1278-1287 (2021)
  * (65) R. Chen, Phys. Rev. D 103, no.5, 054007 (2021)
  * (66) R. Chen, Eur. Phys. J. C 81, no.2, 122 (2021)
  * (67) H. X. Chen, W. Chen, X. Liu and X. H. Liu, Eur. Phys. J. C 81, no.5, 409 (2021)
  * (68) F. Z. Peng, M. J. Yan, M. Sánchez Sánchez and M. P. Valderrama, Eur. Phys. J. C 81, 666 (2021)
  * (69) Z. G. Wang, Int. J. Mod. Phys. A 36, no.10, 2150071 (2021)
  * (70) N. A. Tornqvist, Z. Phys. C 61, 525 (1994)
  * (71) N. A. Tornqvist, Nuovo Cim. A 107, 2471 (1994)
  * (72) T. M. Yan, H. Y. Cheng, C. Y. Cheung, G. L. Lin, Y. C. Lin and H. L. Yu, Phys. Rev. D 46, 1148 (1992) [Erratum-ibid. D 55, 5851 (1997)].
  * (73) M. B. Wise, Phys. Rev. D 45, 2188 (1992).
  * (74) G. Burdman and J. F. Donoghue, Phys. Lett. B 280, 287 (1992).
  * (75) R. Casalbuoni, A. Deandrea, N. Di Bartolomeo, R. Gatto, F. Feruglio and G. Nardulli, Phys. Rept. 281, 145 (1997)
  * (76) A. F. Falk and M. E. Luke, Phys. Lett. B 292, 119 (1992)
  * (77) Y. R. Liu and M. Oka, Phys. Rev. D 85, 014015 (2012)
  * (78) P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020)
  * (79) C. Isola, M. Ladisa, G. Nardulli and P. Santorelli, Phys. Rev. D 68, 114001 (2003).
  * (80) W. A. Bardeen, E. J. Eichten, and C. T. Hill, Phys. Rev. D 68, 054024 (2003).
  * (81) R. Chen, A. Hosaka, and X. Liu, Phys. Rev. D 97, 036016 (2018).
  * (82) F. E. Close and P. R. Page, Phys. Lett. B 578, 119-123 (2004)


Advance Access Publication Date: Day Month Year Reflection

Brusselmans et al.

[$\ast$]Corresponding author<EMAIL_ADDRESS>
0Year

# On the importance of assessing topological convergence in Bayesian
phylogenetic inference

Marius Brusselmans  Luiz Max Carvalho  Samuel L. Hong  Jiansi Gao  Frederick
A. Matsen IV0000-0003-0607-6025  Andrew Rambaut  Philippe Lemey  Marc A.
Suchard  Gytis Dudas  Guy Baele0000-0002-1915-7732  Department of
Microbiology, Immunology and Transplantation, Rega Institute, KU Leuven,
Leuven, Belgium  School of Applied Mathematics, Getulio Vargas Foundation,
Praia de Botafogo, 190, 22250-900, Rio de Janeiro, Brazil  Computational
Biology Program, Fred Hutchinson Cancer Center, Seattle, WA 98109, USA  Howard
Hughes Medical Institute, Computational Biology Program, Fred Hutchinson
Cancer Research Center, Seattle, Washington, USA  Department of Genome
Sciences, University of Washington, Seattle, Washington, USA  Department of
Statistics, University of Washington, Seattle, Washington, USA  Institute of
Ecology and Evolution, University of Edinburgh, Edinburgh, EH9, 3FL, UK
Department of Biostatistics, Fielding School of Public Health, University of
California, Los Angeles, CA 90095, USA  Institute of Biotechnology, Life
Sciences Center, Vilnius University, Vilnius, Lithuania

(2024; 2024; Date)

###### Abstract

Modern phylogenetics research is often performed within a Bayesian framework,
using sampling algorithms such as Markov chain Monte Carlo (MCMC) to
approximate the posterior distribution. These algorithms require careful
evaluation of the quality of the generated samples. Within the field of
phylogenetics, one frequently adopted diagnostic approach is to evaluate the
effective sample size (ESS) and to investigate trace graphs of the sampled
parameters. A major limitation of these approaches is that they are developed
for continuous parameters and therefore incompatible with a crucial parameter
in these inferences: the tree topology. Several recent advancements have aimed
at extending these diagnostics to topological space. In this short reflection
paper, we present a case study illustrating how these topological diagnostics
can contain information not found in standard diagnostics, and how decisions
regarding which of these diagnostics to compute can impact inferences
regarding MCMC convergence and mixing. Given the major importance of detecting
convergence and mixing issues in Bayesian phylogenetic analyses, the lack of a
unified approach to this problem warrants further action, especially now that
additional tools are becoming available to researchers.

###### keywords:

effective sample size, topologies, Bayesian inference, phylogenetics,
phylodynamics, convergence, mixing

##  1 Introduction

### Background

When performing Bayesian phylogenetic inference using Markov chain Monte Carlo
(MCMC) algorithms, the standard practice is to visually inspect trace plots
and to compute the corresponding effective sample size (ESS) of the sampled
parameters. This can be done using a variety of software packages such as
Tracer (Rambaut et al., 2018), Beastiary (Wirth and Duchene, 2022) or CODA
(Plummer et al., 2020). However, such software packages only produce
diagnostics for simple, often univariate model parameters, which the topology
of the phylogenetic tree is not. If these diagnostics suggest the MCMC
convergence and mixing of the simple model parameters are satisfactory, it is
usually assumed that this will be the case for the topology as well. This is
potentially problematic, as the tree topology is often of key interest in
phylogenetic and phylodynamic studies, and obtaining a correct (consensus)
phylogeny is essential when performing outbreak investigation and monitoring
ongoing epidemics (see e.g. Attwood et al. (2022)).

Recent research has focused on convergence and mixing diagnostics applicable
to the tree topology, including studies by Lanfear et al. (2016), Guimarães
Fabreti and Höhna (2021), and Magee et al. (2023). The former two studies
focus on finding ways to apply the principles of traces and effective sample
sizes to the topology as a whole, while the latter study considers the
presence of each split in the tree as an individual parameter to be evaluated
using classical diagnostics. We here explore the topological methods
introduced by Lanfear et al. (2016) and Guimarães Fabreti and Höhna (2021), as
well as those from Magee et al. (2023) that were stated in the original
publication to be preferred over others on the data from an Ebola virus (EBOV)
study by Dudas et al. (2017). Of note, we also include the multidimensional
scaling (MDS) ESS metric – a slightly more conservative tree ESS measure
(Magee et al., 2023) – that also enables to project the high-dimensional
phylogenies to a small number of dimensions suitable for visualisation
(Kruskal, 1964a, b).

We find that the evaluation of convergence and mixing of samples in
topological space can reveal issues not typically detected by standard
diagnostics for (convergence and mixing of) continuous parameters. We selected
this EBOV study because of its large size - it being the most densely sampled
viral outbreak prior to SARS-CoV-2 - and the rich complexity of the model that
was applied, making it a prime case study of a challenging phylogenetic
analysis that could be susceptible to hidden convergence and mixing issues.

### The 2016 Ebola virus data

We analyse data from a genomic epidemiology study of the largest EBOV outbreak
to date, which investigated the impact of several potential predictors, such
as climate and demographic information, on EBOV spread in West Africa from
2014 to 2015. Other aspects, such as the impact of international travel
restrictions and reasons for the containment of the epidemic to mostly West
Africa, are also addressed in that study and we refer the interested reader to
the original publication for further details (Dudas et al., 2017). In summary,
Dudas et al. (2017) performed Bayesian phylogenetic inference using MCMC on a
total of 1610 EBOV genome sequences sampled between March 2014 and October
2015 using an HKY + $\Gamma_{4}$ nucleotide substitution model, a relaxed
molecular clock model with an underlying lognormal distribution and a flexible
non-parametric coalescent model (Gill et al., 2013) as the tree prior.

We downloaded the 1 000 posterior sample trees, which were sampled every 10
000 iterations, and log-files from the original publication from:
https://github.com/ebov/space-time/tree/master/Analyses/Phylogenetic. Of note,
the burn-in was already discarded from the posterior sample trees file shared
by Dudas et al. (2017).

### Phylogenetic distance metrics

A central concept in our exploration is that of phylogenetic distance, a
quantitative measure of similarity between two phylogenetic trees. Such a
distance can be computed in several ways, using what we will refer to as
phylogenetic distance metrics. These distances are zero for two identical
trees and are expected to increase as trees grow more dissimilar. The metrics
considered in this paper are: the Robinson-Foulds distance (Robinson and
Foulds, 1981) and its weighted counterpart (Robinson and Foulds, 1979), the
path difference (Steel and Penny, 1993), the branch score (Kuhner and
Felsenstein, 1994), and the Kendall-Colijn distance (Kendall and Colijn, 2016)
with $\lambda$=0 (i.e., disregarding the branch lengths).

### Topological convergence diagnostics

We here briefly discuss the topological convergence diagnostics considered in
this manuscript.

The topology trace plot, as defined by Lanfear et al. (2016), works very
similarly to standard trace plots used for continuous parameters. The value of
the trace being graphed is the phylogenetic distance from each sampled tree to
a chosen reference tree. The choice of reference tree is somewhat arbitrary;
choosing a tree that is part of the chain will cause a “slump” in the trace
plot as the distance from this tree to itself is inevitably equal to zero. It
is good practice to try several reference trees and compare the resulting
trace plots. In this manuscript, we use the first tree as a reference and
exclude it from the graphs to avoid scaling issues.

The pseudo-ESS (Lanfear et al., 2016) is the ESS of the vector of phylogenetic
distances from an arbitrarily chosen focal tree to all other trees in the
sample (note the analogy with the topology trace plot). Because there is
randomness involved in the choice of the focal tree, the computation is
repeated using each tree in the sample as a focal tree. The lowest value and
median value are then reported.

The approximate ESS (Lanfear et al., 2016) adopts the common ESS calculation
of dividing the actual sample size by the autocorrelation time, but it uses a
topological version of autocorrelation time estimated by determining the
thinning interval at which the average phylogenetic distance between
subsequently sampled trees ceases to increase as the thinning interval
increases.

The Fréchet correlation ESS (Magee et al., 2023) is defined analogously to the
standard one-dimensional continuous ESS, but substitutes Pearson
autocorrelation for an alternative definition of autocorrelation between trees
using Fréchet (co)variances, which make use of the relationship between
covariance and Euclidian distance – here substituted with whatever
phylogenetic distance is being used.

The split frequency ESS (Magee et al., 2023) is computed by treating each tree
as a vector of binary split indicators (split is present/absent). Fréchet
(co)variances can then be computed using the Euclidean norm (as the ”trees”
are now reduced to vectors of binary indicators), which allow for the
computation of an ESS. Note that this is the only method that does not use any
kind of phylogenetic distance metric.

The multidimensional scaling ESS (Magee et al., 2023) (MDS ESS) is computed by
performing multidimensional scaling of the matrix of squared pairwise
distances between trees. The first dimension of the resulting multidimensional
scaling matrix is then used to compute an ESS on.

##  2 Methods

We computed the topological ESS estimators described in the previous section
using the treess package version 1.0.1 (Magee et al., 2023) for the CRAN R
software package v4.3.0 (R-team, 2022), the phylogenetic distances required
for these ESS estimators using the R packages phangorn v2.11.1 (Schliep et
al., 2022) and TreeDist v2.6.1 (Martin, 2023), the per-split ESS’s with the R
package convenience (Guimarães Fabreti and Höhna, 2021), and the approximate
Subtree-Prune-Regraft (aSPR) distances using RSPR version 1.3.1 (Whidden et
al., 2013). We constructed the tanglegrams using Baltic version 0.2.2,
available at: https://github.com/evogytis/baltic. Finally, we also made use of
Tracer v1.7.2 (Rambaut et al., 2018) and the TreeAnnotator v1.10.4 tool
associated with the BEAST 1.10.4 software package for summarising maximum
clade credibility (MCC) trees (Suchard et al., 2018), as well as the R package
ggtree v3.8.2 (Xu et al., 2022) for visualisation of these trees.

##  3 Topological trace plots and ESS values

Phylogenetic distances from first tree

Figure 1: Topology trace plots for the Ebola analysis using five different
phylogenetic distance metrics. The X-axis shows the index of the sampled tree
(thinned sample), while the Y-axis shows the phylogenetic distance from the
first tree. Red vertical lines denote iteration 333 and 666. Note how certain
distance metrics show clear jumps between the separate runs, while for others
the traces seem totally homogeneous.

Figure 1 shows the topology trace plots for the MCMC sample using five
different phylogenetic distance metrics. What is immediately apparent from the
thinned trace plot is the fact that, using certain distance metrics, the
posterior sample of trees seems to be divided into three distinct parts at
iterations 333 and 666. This suggests that three independent replicate
analyses were combined to obtain a posterior samples of trees in the study of
Dudas et al. (2017), which has been confirmed by the respective authors.
Bayesian phylogenetic inference on large data sets indeed commonly employs the
practice of concatenating the samples of several independent chains to both
reduce computation time (by increasing the ESS values of continuous
parameters) and assess convergence towards the same posterior. Three of the
topology trace plots in Figure 1 therefore suggest a discrepancy between the
posterior space explored by the three chains. Whether this indicates failure
to converge to the same posterior, a case of extremely slow / poor mixing, or
something else entirely is not clear. Intriguingly, this behaviour is not
detected at all using certain other metrics, more specifically the path
difference and Kendall-Colijn distance (with $\lambda$=0).

Subsample | 1-1000 | 1-333 | 334-666 | 667-1000
---|---|---|---|---
Estimator (RF) |  |  |  |
Approximate ESS | 362 | 181 | 150 | 143
Split frequency ESS | 446 | 279 | 273 | 276
Fréchet correlation ESS | 83 | 229 | 209 | 211
Median pseudo-ESS | 44 | 333 | 334 | 333
Minimum pseudo-ESS | 15 | 276 | 285 | 271
MDS ESS | 2 | 7 | 3 | 5
Estimator (wRF) |  |  |  |
Approximate ESS | 864 | 299 | 292 | 296
Split frequency ESS | 446 | 279 | 273 | 276
Fréchet correlation ESS | 46 | 238 | 189 | 212
Median pseudo-ESS | 12 | 282 | 280 | 282
Minimum pseudo-ESS | 3 | 189 | 189 | 177
MDS ESS | 2 | 20 | 8 | 6
Estimator (PD) |  |  |  |
Approximate ESS | 1000 | 333 | 334 | 333
Split frequency ESS | 446 | 279 | 273 | 276
Fréchet correlation ESS | 280 | 292 | 279 | 299
Median pseudo-ESS | 1000 | 333 | 334 | 333
Minimum pseudo-ESS | 631 | 205 | 166 | 239
MDS ESS | 1000 | 333 | 334 | 333
Estimator (BS) |  |  |  |
Approximate ESS | 710 | 290 | 278 | 287
Split frequency ESS | 446 | 279 | 273 | 276
Fréchet correlation ESS | 22 | 222 | 139 | 189
Median pseudo-ESS | 4 | 261 | 230 | 256
Minimum pseudo-ESS | 2 | 187 | 108 | 89
MDS ESS | 2 | 88 | 12 | 6
Estimator (KC) |  |  |  |
Approximate ESS | 1000 | 333 | 334 | 333
Split frequency ESS | 446 | 302 | 296 | 226
Fréchet correlation ESS | 593 | 279 | 273 | 275
Median pseudo-ESS | 1000 | 333 | 334 | 333
Minimum pseudo-ESS | 518 | 261 | 148 | 181
MDS ESS | 861 | 333 | 280 | 272
Table 1: ESS estimates from the treess package for the full sample and the
three sub-samples. Values below the often used (although somewhat arbitrary)
threshold of 200 are made bold. The phylogenetic distance metrics are here
labeled using the following codes: RF = Robinson-Foulds, wRF = weighted
Robinson-Foulds, PD = Path Difference, BS = Branch Score, KC = Kendall-Colijn.
Like the topology trace plots, ESS’s computed using KC or PD do not show the
same discrepancies as ESS’s computed using other metrics. The different
estimators also do not all behave in the same way.

Applying the standard treess procedure for the different topological ESS
estimators and using the different phylogenetic distance metrics produces a
wide range of results which are listed in table 1. The ESS of the complete
sample is not always greater than the sum of the ESS’s of the subsamples,
which can be expected if the subsamples explore differents parts of tree
space. Another distinction similar to the one in the topology trace plots can
be noted: the path difference and Kendall-Colijn distance do not seem to pick
up the same patterns as the other metrics, and generally seem overly
optimistic about mixing based on the ESS values.

##  4 ESS values of individual tree splits

The approach by Guimarães Fabreti and Höhna (2021) does not aim to compute an
ESS for the entire topology. Instead, each split in the tree is considered
individually as a binary parameter (1=present, 0=absent) for which an ESS is
computed using standard methodology. Figure 2 shows the cumulative density of
ESS values of the $6859$ splits observed in trees generated by the MCMC
algorithm. The vast majority of these (94%) have an ESS above the often used
cutoff value of $200$. Almost half (46%) have an ESS of 1000, equal to the
actual sample size. Although not pictured here, the three subsamples showed a
nearly identically shaped distribution of individual split ESS’s. From this
perspective, the mixing in topological space is quite optimistic.

Figure 2: Cumulative density of ESS’s of individual splits in the sampled
trees. ESS values 200 and 1000 are indicated and connected to their
corresponding cumulative density values with dotted lines. Less than 6% of the
splits have an ESS below the often used threshold of 200, suggesting
satisfactory convergence in topological space. These ESS values were computed
using the convenience R package.

##  5 Trace plots and ESS values of continuous parameters and tree statistics

Figure 3 shows trace graphs of the age of the root node and the tree length.
The height of the root - a continuous parameter often evaluated using standard
methodology - shows a well mixed chain with a satisfactory ESS. The trace of
the tree length, which is more rarely assessed as it is not typically a
statistic of interest, shows three somewhat distinguishable subsamples, as
well as a very poor ESS. This is more in-line with the topological
diagnostics.

Figure 3: Trace plots for the root node height (”Tree height”) and tree length
(”Tree length”). Tree length refers to the sum of all tree branches, and is
thus a statistic closely related to the topology of the tree. The root height,
a commonly evaluated parameter, shows satisfactory mixing, while the tree
length, which is rarely evaluated, shows poor mixing.

##  6 MCC trees

The observed discrepancies in the explored posterior topological space seen in
the sample raises the question of what these differences actually reflect in
terms of phylogenetic inference. In order to obtain some insights into this,
we can compare the topology of the Maximum-Clade-Credibility (MCC) trees
between the three subsamples.

Figure 4 shows a tanglegram comparing the MCC trees of the three subsamples
and the entire sample by linking the corresponding taxa to each other. Taxa
are coloured by country of origin. The figure also contains the MCC tree of
samples 334 - 1000.The ESS values generated by tracer when only considering
the last two subsections are also noticeably better than for the full sample:
547 and 272 for root height and tree length respectively. The trees all have
the same overall shape, but several clades end up in different locations
depending on the subsample. An obvious example of this is the clade of taxa
largely sampled in Guinea with relatively recent sampling dates that switched
branches with a closely clustered clade of Sierra-Leone based taxa. There also
appears to be substantial uncertainty about the locations of nodes near the
root of the tree.

Figure 4: Tanglegram of the MCC trees of the three subsamples, the combination
of subsamples 2 and 3, and the whole sample. The tips of the trees are
connected to each other by lines, coloured by the country of origin.
Disagreement between the subsamples regarding the location of several clades
is apparent by the fact that the lines connecting the tips of these clades are
not parallel.

##  7 Approximate SPR distances

Another commonly used phylogenetic distance metric is the subtree prune-and-
regraft (SPR) distance. This metric is often preferred by phylogeneticists
because it reflects the minimum number of SPR operations required to go from
one tree to another, which is more interpretable than other distance metrics.
These and similar operations are often used to search phylogenetic tree space.

Calculation of SPR distances is computationally intensive for large trees such
as the Ebola trees considered in this paper. In our case, calculation of a
single SPR distance could last anywhere from 10 minutes to more than a day.
This is substantially longer than any other metric discussed in the paper, and
makes it unusable for our purposes. Therefore, we used the approximate SPR
distances (aSPR) instead ( de Oliveira Martins et al. (2008)), whose
computation time - similarly to the five other distance metrics - was in the
seconds. Figure 5 shows a heat map of all the pairwise aSPR distances between
sampled trees. Distances above 1300 are coloured grey. The three sections of
the combined sample can be clearly distinguished by the block structure, as
the distances between them tend to be substantially larger than the
differences within them, corroborating the findings in figure 1. The
differences between the three chains in terms of SPR moves are substantial and
indicate that they in fact explore different modes in the posterior
distribution of trees.

Figure 5: Pairwise approximate SPR (aSPR) distances between trees. Distances
above 1300 are coloured grey. Discrepancies between the three subsamples are
particularly pronounced, as the number of aSPR moves between trees is
substantially larger between the subsamples than within them.

##  8 Conclusion

In this case study, we applied an ensemble of diagnostics in topological space
on a set of large trees, obtained through Bayesian phylogenetic inference
under a set of complex models and through combining the output of three
independent replicate analyses. Importantly, we find that the different
approaches to evaluating topological convergence can lead to drastically
different conclusions, a finding that to the best of our understanding has not
been observed to this extent before. It is likely that a combination of the
complexity of the model, the size of the dataset, and the fact that we looked
at samples that result from three different / independent analysis replicates
– as opposed to only just one – allowed us to observe these phenomena that
went unnoticed before. This finding stresses the importance of assessing
topological convergence to the posterior and not merely continuous parameter
and (joint) density convergence, which is the current benchmark in nearly all
Bayesian phylogenetic and phylodynamic studies.

The various diagnostics we employed have revealed discrepancies between the
posterior spaces explored by the MCMC algorithm that would typically have gone
unnoticed using standard diagnostics for continuous parameters. Whether these
discrepancies affect inferences on estimates of parameters of interest
downstream in the analysis is not yet clear, and warrants further research.
The impact of decisions such as which ESS estimator and topological distance
metrics to use is not yet fully understood, and it is likely that these
different approaches capture different aspects of convergence and mixing. SPR
distances between the trees did suggest that the trees explored by the three
subsamples differ in a topologically substantial way, which is line with the
findings from the topological ESS estimators and the topology trace plots
using (weighted) Robinson-Foulds distances and the branch score, as opposed to
the other metrics we employed.

Graphing a topological trace plot using the (weighted) Robinson-Foulds
distance or the branch score seemed to be the most informative approach for
this data set, as these metrics were capable of discriminating between the
three subsamples and suggested a discrepancy in terms of explored topological
space. Analogously, the topological ESS’s estimators computed with those
metrics tended to be more pessimistic - or realistic - about the information
contained in the sample. These methods can rather easily implemented in
existing software packages, which could improve current practice regarding
convergence assessment of Bayesian phylogenetic and phylodynamic analyses.

Finally, our findings raise important questions as to how the output of
replicate Bayesian phylogenetic and phylodynamic analyses should be combined
when discrepancies in topological space are detected. These replicate analyses
produce approximate samples from different regions of the posterior
distribution. A key open question stemming from our work is how to combine
these samples in such a way that preserves the mode structure, i.e., such that
samples from higher modes are present in the combined sample more often than
samples from the lower modes. Simply combining the samples with equal weights
- which is currently standard practice - might not produce samples from the
correct posterior, but techniques such as importance sampling could be
employed to weight samples proportionally to their density.

##  9 Competing interests

No competing interest is declared.

##  10 Author contributions statement

M.B. and G.B. initialised the study. M.B. performed the analyses and wrote the
manuscript. G.D. created the tanglegrams and computed the aSPR distances.
L.M.C., J.G., F.A.M., A.R., M.A.S., G.D., S.L., P.L., and G.B. provided
valuable guidance and feedback. All authors reviewed the manuscript.

##  11 Acknowledgments

The authors thank Andrew Magee for their valuable technical support. S.L.H.
and G.B. acknowledge support from the Research Foundation - Flanders (“Fonds
voor Wetenschappelijk Onderzoek - Vlaanderen,” G0E1420N). G.B. acknowledges
support from the Internal Funds KU Leuven (Grant No. C14/18/094), from the
Research Foundation - Flanders (“Fonds voor Wetenschappelijk Onderzoek -
Vlaanderen,” G098321N) and from the European Union Horizon 2023 RIA project
LEAPS (grant agreement no. 101094685). P.L., M.A.S. and A.R. acknowledge
support from the Wellcome Trust (Collaborators Award 206298/Z/17/Z, ARTIC
network), the European Research Council (grant agreement no. 725422 –
ReservoirDOCS) and the National Institutes of Health (NIH) (R01 AI153044).
F.A.M. and M.A.S. acknowledge further support from the NIH through R01
AI162611. P.L. acknowledges support from the Research Foundation, Flanders
(“Fonds voor Wetenschappelijk Onderzoek - Vlaanderen,” G066215N, G0D5117N and
G0B9317N) and from the European Union Horizon 2020 project MOOD (grant
agreement no. 874850). M.B. and G.B. acknowledge support from the DURABLE
EU4Health project 02/2023-01/2027, which is co-funded by the European Union
(call EU4H-2021-PJ4) under Grant Agreement No. 101102733. Views and opinions
expressed are however those of the author(s) only and do not necessarily
reflect those of the European Union or the European Health and Digital
Executive Agency. Neither the European Union nor the granting authority can be
held responsible for them.

## References

  * Attwood et al. [2022] S. Attwood, S. Hill, D. Aanensen, T. Connor, and O. Pybus.  Phylogenetic and phylodynamic approaches to understanding and combating the early sars-cov-2 pandemic.  _Nature Reviews Genetics_ , 23:1–16, 04 2022.  10.1038/s41576-022-00483-8.
  * de Oliveira Martins et al. [2008] L. de Oliveira Martins, E. Leal, and H. Kishino.  Phylogenetic detection of recombination with a Bayesian prior on the distance between trees.  _PloS one_ , 3(7):e2651–e2651, 2008.
  * Dudas et al. [2017] G. Dudas, L. M. Carvalho, T. Bedford, A. J. Tatem, G. Baele, N. R. Faria, D. J. Park, J. T. Ladner, A. Arias, D. Asogun, F. Bielejec, S. L. Caddy, M. Cotten, J. D’Ambrozio, S. Dellicour, A. D. Caro, J. W. Diclaro, S. Duraffour, M. J. Elmore, L. S. Fakoli, O. Faye, M. L. Gilbert, S. M. Gevao, S. Gire, A. Gladden-Young, A. Gnirke, A. Goba, D. S. Grant, B. L. Haagmans, J. A. Hiscox, U. Jah, J. R. Kugelman, D. Liu, J. Lu, C. M. Malboeuf, S. Mate, D. A. Matthews, C. B. Matranga, L. W. Meredith, J. Qu, J. Quick, S. D. Pas, M. V. T. Phan, G. Pollakis, C. B. Reusken, M. Sanchez-Lockhart, S. F. Schaffner, J. S. Schieffelin, R. S. Sealfon, E. Simon-Loriere, S. L. Smits, K. Stoecker, L. Thorne, E. A. Tobin, M. A. Vandi, S. J. Watson, K. West, S. Whitmer, M. R. Wiley, S. M. Winnicki, S. Wohl, R. Wölfel, N. L. Yozwiak, K. G. Andersen, S. O. Blyden, F. Bolay, M. W. Carroll, B. Dahn, B. Diallo, P. Formenty, C. Fraser, G. F. Gao, R. F. Garry, I. Goodfellow, S. Günther, C. T. Happi, E. C. Holmes, B. Kargbo, S. Keïta, P. Kellam, M. P. G. Koopmans, J. H. Kuhn, N. J. Loman, N. Magassouba, D. Naidoo, S. T. Nichol, T. Nyenswah, G. Palacios, O. G. Pybus, P. C. Sabeti, A. Sall, U. Ströher, I. Wurie, M. A. Suchard, P. Lemey, and A. Rambaut.  Virus genomes reveal factors that spread and sustained the Ebola epidemic.  _Nature_ , 544(7650):309–315, Apr 2017.  10.1038/nature22040.  URL https://doi.org/10.1038%2Fnature22040.
  * Epskamp et al. [2023] D. Epskamp, G. Costantini, J. Haslbeck, A. Isvoranu, A. Cramer, L. Waldorp, V. Schmittmann, and D. Borsboom.  qgraph: Graph plotting methods, psychometric data visualization and graphical model estimation.  _CRAN repository_ , 2023.
  * Gill et al. [2013] M. Gill, P. Lemey, N. Faria, A. Rambaut, B. Shapiro, and M. A. Suchard.  Improving Bayesian population dynamics inference: A coalescent-based model for multiple loci.  _Molecular Biological Evolution_ , 30(3):713–724, 2013.
  * Guimarães Fabreti and Höhna [2021] L. Guimarães Fabreti and S. Höhna.  Convergence assessment for Bayesian phylogenetic analysis using MCMC simulation.  _Methods in Ecology and Evolution_ , 13(1):77–90, 2021.
  * Kendall and Colijn [2016] M. Kendall and C. Colijn.  Mapping phylogenetic trees to reveal distinct patterns of evolution.  _Molecular Biology and Evolution_ , 33(10):2735–2743, 2016.
  * Kruskal [1964a] J. B. Kruskal.  Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis.  _Psychometrika_ , 29:1–27, 1964a.
  * Kruskal [1964b] J. B. Kruskal.  Nonmetric multidimensional scaling: A numerical method.  _Psychometrika_ , 29:115–129, 1964b.
  * Kuhner and Felsenstein [1994] M. Kuhner and J. Felsenstein.  A simulation comparison of phylogeny algorithms under equal and unequal evolutionary rates.  _Molecular Biology and Evolution_ , 11(3):459–468, 1994.
  * Lanfear et al. [2016] R. Lanfear, X. Hua, and D. L. Warren.  Estimating the effective sample size of tree topologies from Bayesian phylogenetic analyses.  _Genome Biology and Evolution_ , 8(8):2319–2332, 2016.
  * Magee et al. [2023] A. Magee, M. Karcher, F. Matsen IV, and V. Minin.  How trustworthy is your tree? Bayesian phylogenetic effective sample size through the lens of Monte Carlo error.  _Bayesian Analysis_ , 1(1):1–29, 2023.
  * Martin [2023] S. R. Martin.  Treedist: Calculate and map distances between phylogenetic trees.  2023\.
  * Plummer et al. [2020] M. Plummer, N. Best, K. Cowles, K. Vines, D. Sarkay, D. Bates, R. Almond, and A. Magnusson.  Package ’coda’.  _CRAN_ , 2020.
  * R-team [2022] R-team.  R: A language and environment for statistical computing.  2022\.
  * Rambaut et al. [2018] A. Rambaut, A. Drummond, D. Xie, G. Baele, and M. Suchard.  Posterior summarisation in Bayesian phylogenetics using Tracer 1.7 (available at http://beast.community/tracer).  _Systematic Biology_ , 67(5):901–904, 2018.
  * Robinson and Foulds [1979] D. Robinson and L. Foulds.  Comparison of weighted labelled trees.  _In: Horadam, A.F., Wallis, W.D. (eds) Combinatorial Mathematics VI_ , 748:119–126, 1979.
  * Robinson and Foulds [1981] D. Robinson and L. Foulds.  Comparison of phylogenetic trees.  _Mathematical Biosciences_ , 53(1-2):131–147, 1981.
  * Schliep et al. [2022] K. Schliep, E. Paradis, L. de Oliveira Martins, A. Potts, and T. White.  Package ’phangorn’.  _CRAN_ , 2022.
  * Steel and Penny [1993] M. A. Steel and D. Penny.  Distributions of tree comparison metrics - some new results.  _Syst Biol._ , 42:126–141, 1993.
  * Suchard et al. [2018] M. A. Suchard, P. Lemey, G. Baele, A. J. Drummond, and A. Rambaut.  Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10.  _Virus Evolution_ , 4(1):vey016, 2018.
  * Whidden et al. [2013] C. Whidden, R. Beiko, and N. Zeh.  Fixed-parameter algorithms for maximum agreement forests.  _SIAM Journal on Computing_ , 42(4):1431–1466, 2013.
  * Wirth and Duchene [2022] W. Wirth and S. Duchene.  Real-Time and Remote MCMC Trace Inspection with Beastiary.  _Molecular Biology and Evolution_ , 39(5):msac095, 05 2022.  ISSN 1537-1719.  10.1093/molbev/msac095.  URL https://doi.org/10.1093/molbev/msac095.
  * Xu et al. [2022] S. Xu, L. Li, X. Luo, M. Chen, W. Tang, L. Zhan, Z. Dai, T. Lam, Y. Guan, and G. Yu.  Ggtree: A serialized data object for visualization of a phylogenetic tree and annotation data.  _iMeta_ , 1(e56), 2022.

##  12 Supplementary Materials

Figure S1: Tanglegram of the MCC trees of the three subsamples, the latter two
subsamples combined, and the whole sample. Nodes with posterior support $<$0.1
are collapsed. When collapsing the nodes with low posterior support ($<$10%),
the differences in topology between the subsamples seem less pronounced.
However, adding the first subsamples to the latter two (MCC tree ‘23’)
strongly impacts the (structure of the) overall MCC tree (‘total’).

((a)) Subsample 1-333

((b)) Subsample 334-666

((c)) Subsample 667-1000

((d)) Subsample 334-1000

((e)) Total sample

Figure S2: MCC trees of the different subsamples with taxa coloured by origin
country, analogous to the tanglegrams of Figures 4 and S1. Figure S3: Network
graph of approximate SPR (aSPR) distances between sampled trees. Each node
represents a sampled tree and distances between nodes approximate the aSPR
distance between the trees. Nodes are coloured by subsection of the sample.
The edges correspond to trees that were subsequently sampled. This graph was
generated using a force-directed layout based on the Fruchterman-Reingold
algorithm, as implemented in the R package qgraph version 1.9.8. (Epskamp et
al. [2023]). The exact location of the nodes has no real meaning, only their
relative locations. As expected from the topological trace plots, the three
subsamples are also distinctly recognisable in this network graph. Note that
the first tree in each subsample is somewhat distinct from the rest of the
trees in that subsample, potentially pointing to a somewhat too low burn-in
while processing the sampled trees from the Bayesian phylogenetic analysis.