# Diabatic Representation of Exotic Hadrons in the Dynamical Diquark Model

Richard F. Lebed<EMAIL_ADDRESS>Steven R. Martinez<EMAIL_ADDRESS>Department of Physics, Arizona State University, Tempe, AZ 85287, USA

(July, 2022)

###### Abstract

We apply the diabatic formalism, an extension of the adiabatic approximation
inherent to the Born-Oppenheimer (BO) approach of atomic physics, to the
problem of mixing between exotic multiquark hadrons and their nearby di-hadron
thresholds. The unperturbed BO eigenstates are obtained using the dynamical
diquark model, while the diabatic calculation introduces a mixing potential
between these states and the threshold states. We solve the resulting coupled
Schrödinger equations numerically for hidden-charm tetraquarks of both open
and closed strangeness to obtain physical mass eigenvalues, and explore the
di-hadron state content and spatial extent of the eigenstates. As an explicit
example, $X(3872)$ emerges with a dominant $D^{0}\bar{D}^{*0}$ component, but
also contains a considerable diquark-antidiquark component that can contribute
significantly to its radiative decay widths, and this component also generates
a full multiplet of other diquark-based exotic hadrons to be compared with
experiment.

Exotic hadrons, diquarks

##  I Introduction

With the approach of the two-decade mark for experimental evidence of heavy-
quark exotic hadrons Choi  _et al._ (2003), the field remains at a remarkable
point: Well over 50 candidates have been observed with a high degree of
statistical significance, and yet no single theoretical paradigm
satisfactorily accounts for all of them Lebed  _et al._ (2017); Chen  _et al._
(2016); Hosaka  _et al._ (2016); Esposito  _et al._ (2017); Guo  _et al._
(2018); Ali  _et al._ (2017); Olsen  _et al._ (2018); Karliner  _et al._
(2018); Yuan (2018); Liu  _et al._ (2019); Brambilla  _et al._ (2020); Chen
_et al._ (2022). Even so, several of these states lie exceptionally close to
di-hadron thresholds, the most extraordinary example being the first one
discovered, $X(3872)$:

$m_{X(3872)}\,–\,m_{D^{0}}\,–\,m_{D^{*0}}=-0.04\pm 0.09\ {\rm MeV}\,,$ (1)

using the average mass value for each hadron in Eq. (1) tabulated by the
Particle Data Group (PDG) Zyla  _et al._ (2020).

No other hadron lies in such close proximity to a decay threshold, suggesting
a unique importance for the $D^{0}{\bar{D}}^{*0}$ component (the inclusion of
charge conjugates being understood throughout) within the state. At minimum, a
state so close to threshold exhibits a number of “universality” properties
that depend only upon the large $D^{0}{\bar{D}}^{*0}$ scattering length
Braaten and Kusunoki (2004). More typically, $X(3872)$ has frequently been
considered as a hadronic molecule analogous to the deuteron $d$, which is
largely (but not entirely) bound by $\pi$ exchanges between its nucleon
components. However, the 2.23 MeV $d$ binding energy is many times larger than
Eq. (1), and corresponds to a typical hadronic size comparable to its root-
mean-square (rms) charge radius $\sqrt{\\!\left<r_{d}^{2}\right>}\\!=2.13$ fm
Tiesinga  _et al._ (2021). By the same token, $X(3872)$ would be expected to
be many femtometers across, larger than most nuclei, and its physical
observables would be utterly dominated by long-distance
$D^{0}$-${\bar{D}}^{*0}$ interactions.

Nevertheless, $X(3872)$ exhibits certain properties suggesting the
significance of its short-distance wave-function components. Its known decays
to conventional charmonium ($J/\psi$ and $\chi_{c1}$) account for more than
10% of $\Gamma_{X(3872)}$, and one of its observed radiative decay modes,
$\gamma\psi(2S)$, has a branching fraction of at least a few percent Zyla  _et
al._ (2020). Since charmonium rms charge radii are predicted from quark-
potential models to be significantly less than 1 fm, the large expected
separation of the charm quarks in the $D^{0}$-${\bar{D}}^{*0}$ pair would
naively predict much smaller branching ratios for these processes. But
$X(3872)$ also shares the $J^{PC}$ quantum numbers $1^{++}$ of the yet-unseen
conventional charmonium state $\chi_{c1}(2P)$ that quark-potential models
predict (e.g., Ref. Barnes  _et al._ (2005)) to lie several tens of MeV above
$m_{X(3872)}$. As a result, $X(3872)$ has long been suggested to contain a
substantial core of conventional $\chi_{c1}(2P)$ Suzuki (2005), in part as a
mechanism to explain its surprising decay patterns.

Nevertheless, conventional charmonium is not the only short-distance component
available to $X(3872)$. The valence quark content $c\bar{c}u\bar{u}$ of
$D^{0}{\bar{D}}^{*0}$ allows for an alternative binding mechanism: that of a
diquark-antidiquark $(cu)_{\bar{\bm{3}}}(\bar{c}\bar{u})_{\bm{3}}$ pair, each
one bound through the attractive color channel
${\bm{3}}\otimes{\bm{3}}\\!\to\\!\bar{\bm{3}}$.

A number of other exotic candidates lie quite close to (within a few MeV of)
di-hadron thresholds, notably, $Z_{c}(3900)$ [$D{\bar{D}}^{*}$], $Z_{c}(4020)$
[$D^{*}{\bar{D}}^{*}$], $P_{c}(4312)$ [$\Sigma_{c}{\bar{D}}$],
$P_{c}(4450)$/$P_{c}(4457)$ [$\Sigma_{c}{\bar{D}}^{*}$], $Z_{b}(10610)$
[$B{\bar{B}}^{*}$], $Z_{b}(10650)$ [$B^{*}{\bar{B}}^{*}$], and others. In
fact, all of these $Z$ meson states lie slightly above the corresponding
thresholds, arguing against a traditional bound-state molecular picture. And
yet, one cannot deny the significance of the proximity of the thresholds in
these cases, suggesting a special importance of those particular hadron pairs
for the exotic state. On the other hand, some exotic candidates [e.g., the $Y$
states or $Z_{c}(4430)$] lack an obvious nearby threshold. A complete
theoretical framework accommodating all of the heavy-quark exotic candidates
must therefore recognize the physical significance of such nearby thresholds
on their formation, mass, and decay modes.

One formalism that has predicted a specific spectrum of multiquark heavy-quark
exotic hadrons is the dynamical diquark model Brodsky  _et al._ (2014); Lebed
(2017), in which the formation of a diquark-antidiquark
($\delta^{\vphantom{1}}$-${\bar{\delta}}$) exotic meson requires the color-
nonsinglet diquark quasiparticles $\delta^{\vphantom{1}},{\bar{\delta}}$ not
to dissociate instantly into a di-meson pair, but rather to persist as
components of a single multiquark state connected by a color flux tube. Since
each diquark is a color triplet containing a heavy quark, the same potentials
that are computed on the lattice to describe heavy quarkonium and its hybrids
can be imported into Schrödinger equations, which are solved numerically to
obtain the spectrum of $\delta^{\vphantom{1}}{\bar{\delta}}$ eigenstates. This
procedure has been performed both for the multiplet average masses Giron  _et
al._ (2019), and for the detailed spectrum once spin and isospin fine-
structure effects are included Giron  _et al._ (2020); Giron and Lebed (2020a,
b); Giron  _et al._ (2021). Pentaquarks are handled similarly Giron and Lebed
(2021), by replacing the antidiquark ${\bar{\delta}}$ with a color-triplet
triquark ${\bar{\theta}}\\!\equiv\\![\bar{Q}_{\bar{\bf
3}}(q_{1}q_{2})_{\bar{\bf 3}}]_{\bf 3}$ Lebed (2015).

This treatment of obtaining eigenstates of heavy (quasi-)particles at
separation $r$ connected by a static potential manifests the well-known Born-
Oppenheimer (BO) approximation Born and Oppenheimer (1927) from atomic
physics. Intrinsic to the approximation is the assumption that the light
degrees of freedom (d.o.f.) of the state adjust instantly to changes in the
configuration $\Gamma$ of the heavy sources (an adiabatic approximation), and
that the eigenstates of the potential $V_{\Gamma}(r)$ thus derived from each
such configuration $\Gamma$ change gradually with $V_{\Gamma}(r)$. However,
when the value of $V_{\Gamma}(r)$ crosses the energy of a di-hadron threshold,
the physical eigenstates undergo a rapid level crossing between a
predominantly $\delta^{\vphantom{1}}{\bar{\delta}}$ state and a predominantly
di-hadron state. The mixing of configurations induces diabatic changes to the
system, and solving for its new eigenstates requires a generalization beyond
the strict BO limiting case.

Such a diabatic formalism has been extensively developed in atomic physics,
and in recent years it has been standardized into textbook form Baer (2006).
This approach provides a specific, nonperturbative method for incorporating
the mixing of coupled-channel contributions from di-hadron states into the
calculations, and we briefly review the relevant formalism below. The diabatic
formalism was first applied to heavy-quark systems quite recently in Ref.
Bruschini and González (2020); in that case, the analysis examined mixing of
di-hadron thresholds with conventional quarkonium, rather than with the
4-quark $\delta^{\vphantom{1}}{\bar{\delta}}$ states of the dynamical diquark
model that are used here.

The initial study of Ref. Bruschini and González (2020) specifically considers
states lying below or just above a di-hadron threshold. In this work we adopt
the same restriction, specifically to study the lightest hidden-charm exotics
of valence quark content $c\bar{c}q\bar{q}^{\prime}$ (where $q^{(\prime)}$ are
$u$ or $d$ quarks), $c\bar{c}s\bar{s}$, and $c\bar{c}q\bar{s}$. “Lightest” in
this sense means the members of the (positive-parity) ground-state BO
multiplet $\Sigma^{+}_{g}(1S)$; observed candidates with (presumed)
corresponding flavor contents include $X(3872)$, $Y(4140)$, and
$Z_{cs}(4000)$, respectively. Of course, the same methods can be applied as
well to orbitally excited multiplets such as $\Sigma^{+}_{g}(1P)$ [containing,
e.g., $Y(4220)$], hidden-bottom [e.g., $Z_{b}(10610)$] states, and fully
charmed $c\bar{c}c\bar{c}$ [e.g., $X(6900)$] states, and can also be
generalized to resonant states in order to study mass shifts and strong-decay
widths Bruschini and González (2021a).

Here we focus on identifying the $\delta^{\vphantom{1}}$-${\bar{\delta}}$ and
meson-meson content (including distinct contributing partial waves) of the
lightest hidden-charm states, and extract interesting features such as the
expectation values $\left<r\right>$ and $\langle r^{2}\rangle^{1/2}$ of the
$\delta^{\vphantom{1}}$-${\bar{\delta}}$ separation $r$ in each mass
eigenstate. We show, for example, that the wave function of $X(3872)$ is
indeed dominated by $D^{0}{\bar{D}}^{*0}$, but not overwhelmingly so, and find
that the spectrum of states described by the original uncoupled dynamical
diquark model is not greatly disrupted by the existence of di-hadron
thresholds.

This paper is organized as follows. In Sec. II we review the state notation
for $\delta^{\vphantom{1}}{\bar{\delta}}$ systems, focusing for now only on
the $\Sigma^{+}_{g}(1S)$ multiplet. Section III reviews the diabatic mixing
formalism needed for the current set of calculations. Our numerical results
appear in Sec. IV, and in in Sec. V we summarize and indicate the next
directions for future calculations.

##  II States of the dynamical diquark Model

The complete spectra of ${c\bar{c}q\bar{q}^{\prime}},c\bar{c}s\bar{s}$,
$c\bar{c}q\bar{s}$, and $c\bar{c}c\bar{c}$ states as
$\delta^{\vphantom{1}}{\bar{\delta}}$ eigenstates of the dynamical diquark
model are presented in Refs. Lebed (2017), Giron and Lebed (2020b), Giron  _et
al._ (2021), and Giron and Lebed (2020c), respectively. In this paper, all
relevant states are accommodated by the lowest ($\Sigma^{+}_{g}$) Born-
Oppenheimer (BO) potential, which consists of the (light) gluon field in its
ground state connecting the heavy diquark
[$\delta^{\vphantom{1}}\\!\equiv(Qq)_{\bar{\bf 3}}$]-antidiquark
[${\bar{\delta}}\\!\equiv\\!(\bar{Q}\bar{q}^{\prime})_{\mathbf{3}}$] or
diquark-triquark [${\bar{\theta}}\\!\equiv\\!(\bar{Q}_{\bar{\bf
3}}(q_{1}q_{2})_{\bar{\bf 3}})]_{\bf 3}$ quasiparticles. In all cases,
$\delta^{\vphantom{1}},{\bar{\delta}},{\bar{\theta}}$ are assumed to transform
as color triplets (or antitriplets), and each quasiparticle contains no
internal orbital angular momentum.

In this work, we consider ${Q\bar{Q}}q\bar{q}^{\prime}$ states, in which
$q,\bar{q}^{\prime}$ may assume any of the flavors $\\{u,d,s\\}$. The
classification scheme, regardless of the combination, begins with 6 states,
here grouped by $J^{PC}$ quantum numbers. This spectrum,

$\displaystyle J^{PC}=0^{++}:$ $\displaystyle
X_{0}\equiv\left|0_{\delta}^{\vphantom{1}},0_{\bar{\delta}}\right>_{0}\,,\ \
X_{0}^{\prime}\equiv\left|1_{\delta}^{\vphantom{1}},1_{\bar{\delta}}\right>_{0}\,,$
$\displaystyle J^{PC}=1^{++}:$ $\displaystyle
X_{1}\equiv\frac{1}{\sqrt{2}}\left(\left|1_{\delta}^{\vphantom{1}},0_{\bar{\delta}}\right>_{1}\\!+\left|0_{\delta}^{\vphantom{1}},1_{\bar{\delta}}\right>_{1}\right)\,,$
$\displaystyle J^{PC}=1^{+-}:$ $\displaystyle\,Z\
\equiv\frac{1}{\sqrt{2}}\left(\left|1_{\delta}^{\vphantom{1}},0_{\bar{\delta}}\right>_{1}\\!-\left|0_{\delta}^{\vphantom{1}},1_{\bar{\delta}}\right>_{1}\right)\,,$
$\displaystyle\,Z^{\prime}\equiv\left|1_{\delta}^{\vphantom{1}},1_{\bar{\delta}}\right>_{1}\,,$
$\displaystyle J^{PC}=2^{++}:$ $\displaystyle
X_{2}\equiv\left|1_{\delta}^{\vphantom{1}},1_{\bar{\delta}}\right>_{2}\,,$ (2)

which specifies the full multiplet of $\Sigma^{+}_{g}$ $S$-wave states, is
written with the total $\delta^{\vphantom{1}}({\bar{\delta}})$ spin denoted by
$s^{\vphantom{\dagger}}_{\delta^{\vphantom{1}}}(s_{{\bar{\delta}}})$, and with
the overall state total spin signified by an outer subscript. When needed,
$9j$ angular momentum recoupling coefficients may be used to transform these
states to another spin basis. For example, the transformation coefficients to
the basis of good total heavy-quark (${Q\bar{Q}}$) and light-quark
(${q\bar{q}^{\prime}}$) spin read

$\displaystyle\left<(s_{q}\,s_{\bar{q}})s_{q\bar{q}^{\prime}},(s_{Q}\,s_{\bar{Q}})s_{Q\bar{Q}},S\,\right|\left.(s_{q}\,s_{Q})s_{\delta}^{\vphantom{1}},(s_{\bar{q}}\,s_{\bar{Q}})s_{\bar{\delta}},S\right>$
$\displaystyle=\left([s_{q\bar{q}^{\prime}}][s_{Q\bar{Q}}][s_{\delta}^{\vphantom{1}}][s_{\bar{\delta}}]\right)^{1/2}\begin{Bmatrix}s_{q}&s_{\bar{q}}&s_{q\bar{q}^{\prime}}\\\
s_{Q}&s_{\bar{Q}}&s_{Q\bar{Q}}\\\
s_{\delta}^{\vphantom{1}}&s_{\bar{\delta}}&S\end{Bmatrix}\,,\ \ $ (3)

with $[s]\\!\equiv\\!2s+1$ denoting the multiplicity of a spin-$s$ state.
Using Eqs. (2) and (3), one may write

$\displaystyle J^{PC}=0^{++}:$ $\displaystyle
X_{0}=\frac{1}{2}\left|0_{q\bar{q}^{\prime}},0_{Q\bar{Q}}\right>_{0}+\frac{\sqrt{3}}{2}\left|1_{q\bar{q}^{\prime}},1_{Q\bar{Q}}\right>_{0}\,,$
$\displaystyle
X_{0}^{\prime}=\frac{\sqrt{3}}{2}\left|0_{q\bar{q}^{\prime}},0_{Q\bar{Q}}\right>_{0}-\frac{1}{2}\left|1_{q\bar{q}^{\prime}},1_{Q\bar{Q}}\right>_{0}\,,$
$\displaystyle J^{PC}=1^{++}:$ $\displaystyle
X_{1}=\left|1_{q\bar{q}^{\prime}},1_{Q\bar{Q}}\right>_{1}\,,$ $\displaystyle
J^{PC}=1^{+-}:$
$\displaystyle\,Z\;=\frac{1}{\sqrt{2}}\left(\left|1_{q\bar{q}^{\prime}},0_{Q\bar{Q}}\right>_{1}\\!-\left|0_{q\bar{q}^{\prime}},1_{Q\bar{Q}}\right>_{1}\right)\,,$
$\displaystyle\,Z^{\prime}=\frac{1}{\sqrt{2}}\left(\left|1_{q\bar{q}^{\prime}},0_{Q\bar{Q}}\right>_{1}\\!+\left|0_{q\bar{q}^{\prime}},1_{Q\bar{Q}}\right>_{1}\right)\,,$
$\displaystyle J^{PC}=2^{++}:$ $\displaystyle
X_{2}=\left|1_{q\bar{q}^{\prime}},1_{Q\bar{Q}}\right>_{2}\,.$ (4)

Further specifying the chosen combination of $\\{u,d,s\\}$ light-quark flavors
enlarges this set. For example, considering combinations of $\\{u,d\\}$ alone
expands the set to 12 states: 6 each with $I\\!=\\!0$ and $I\\!=\\!1$, but
which nonetheless maintain spin structures in the forms of Eqs. (2) or (4).
For the purposes of this work, we identify states solely based upon total
$J^{PC}$, effectively ignoring fine structure due to isospin. Additionally, we
separately examine unique flavor combinations of light quarks:
$c\bar{c}s\bar{s}$, $c\bar{c}q\bar{s}$, and ${c\bar{c}q\bar{q}^{\prime}}$,
where henceforth $q,q^{\prime}\in\\{u,d\\}$. The only sources of SU(3)flavor
dependence in these calculations arise through distinct explicit diquark and
meson masses.

All states considered within this work are accommodated within the ground-
state BO multiplet $\Sigma^{+}_{g}(1S)$, but it is worth noting that Ref.
Lebed (2017) provides a classification of states in higher multiplets such as
$\Sigma^{+}_{g}(nP)$, as well as those with excited-glue BO potentials such as
$\Pi_{u}^{+}$.

##  III Diabatic Mixing Formalism

In this work, we begin with the same construction as in the original dynamical
diquark model Giron  _et al._ (2019). That is, one separates the light d.o.f.
from the heavy d.o.f. by writing the Hamiltonian as

$H=K_{\rm heavy}+H_{\rm light}=\frac{\mathbf{p}^{2}}{2\mu_{\rm heavy}}+H_{\rm
light},$ (5)

such that the Schrödinger equation now reads

$\left(\frac{\mathbf{p}^{2}}{2\mu_{\rm heavy}}+H_{\rm
light}-E\right)|\psi\rangle=0.$ (6)

Under the current analysis, “light field” refers to either just the glue
fields (in the case of a $\delta^{\vphantom{1}}$-${\bar{\delta}}$ state) or
both glue and exchanged light-quark fields (in the case of the meson-meson
states).

We now implement the Ansatz that the states defined in Sec. II may appreciably
mix with nearby meson-meson thresholds sharing the same $J^{PC}$ quantum
numbers. Thus, one must determine and solve the multi-channel Schrödinger
equation connecting the $\delta^{\vphantom{1}}{\bar{\delta}}$ states to such
threshold states. We closely follow the work of Ref. Bruschini and González
(2020), which carries out this process using conventional quarkonium rather
than $\delta^{\vphantom{1}}{\bar{\delta}}$ states. Applying the diabatic
expansion to the eigenstates of said Schrödinger equation yields Bruschini and
González (2020)

$|\psi\rangle=\sum_{i}\int
d\mathbf{r}^{\prime}\tilde{\psi}_{i}(\mathbf{r}^{\prime},\mathbf{r}_{0})\>|\mathbf{r}^{\prime}\rangle\>|\xi_{i}(\mathbf{r}_{0})\rangle,$
(7)

where $\mathbf{r}^{\prime}$ denotes the separation of the heavy sources,
$\mathbf{r}_{0}$ is a freely set fiducial parameter, and $|\xi_{i}\rangle$ are
eigenstates of the light-field Hamiltonian. Inserting Eq. (7) into Eq. (6) and
applying $\langle\xi_{j}(\mathbf{r}_{0})|$ on the left-hand side produces

$\sum_{i}\left[-\frac{\hbar^{2}}{2\mu_{i}}\delta^{\vphantom{1}}_{ij}\nabla^{2}+V_{ji}(\mathbf{r,r_{0}})-E\delta^{\vphantom{1}}_{ji}\right]\tilde{\psi}_{i}(\mathbf{r,r_{0}})=0,$
(8)

with the diabatic potential matrix defined as

$V_{ji}(\mathbf{r,r_{0}})\equiv\langle\xi_{j}(\mathbf{r}_{0})|H_{\rm
light}|\xi_{i}(\mathbf{r}_{0})\rangle.$ (9)

We identify the $i\\!=\\!0$ term with $\delta^{\vphantom{1}}{\bar{\delta}}$
states, and $i\\!>\\!0$ terms with meson-meson states.

This result may be written more compactly in matrix notation as

$[\text{K}+\text{V}(\mathbf{r})]\mathbf{\Psi}(\mathbf{r})=\emph{E}\,\mathbf{\Psi}(\mathbf{r})\,,$
(10)

with the parameter $\mathbf{r}_{0}$ implicit. Neglecting interactions between
distinct meson-meson components, as is done in analogous lattice-QCD studies
Bulava  _et al._ (2019), the potential matrix then becomes

$\text{V}=\begin{pmatrix}V_{\delta^{\vphantom{1}}{\bar{\delta}}}(\mathbf{r})&V_{\rm
mix}^{(1)}(\mathbf{r})&\cdots&V_{\rm
mix\vphantom{\overline{M}_{2}}}^{(N)}(\mathbf{r})\\\ V_{\rm
mix}^{(1)}(\mathbf{r})&V_{M_{1}\overline{M}_{2}}^{(1)}(\mathbf{r})&&\\\
\vdots&&\ddots\\\ V_{\rm
mix\vphantom{\overline{M}_{2}}}^{(N)}(\mathbf{r})&&&V_{M_{1}\overline{M}_{2}}^{(N)}(\mathbf{r})\\\
\end{pmatrix},$ (11)

with the kinetic-energy operator expressed as

$\text{K}=\begin{pmatrix}-\frac{\hbar^{2}}{2\mu_{\delta^{\vphantom{1}}{\bar{\delta}}}}&&&\\\
&-\frac{\hbar^{2}}{2\mu^{(1)}_{M_{1}\overline{M}_{2}}}&&\\\ &&\ddots&\\\
&&&-\frac{\hbar^{2}}{2\mu^{(N)}_{M_{1}\overline{M}_{2}}}\\\
\end{pmatrix}\nabla^{2},$ (12)

where omitted elements are zeroes. Inspection of the $i>0$ diagonal elements
leads to the identification of $V_{M_{1}\overline{M}_{2}}^{(i)}(\mathbf{r})$
as simply being the energy associated with that of the pure free $i^{\rm th}$
meson-meson state. That is,

$V_{M_{1}\overline{M}_{2}}^{(i)}(\mathbf{r})=T_{M_{1}\overline{M}_{2}}\,,$
(13)

where

$T_{M_{1}\overline{M}_{2}}\equiv
m^{\vphantom{\dagger}}_{M_{1}}\\!+m_{\overline{M}_{2}}\,.$ (14)

For $S$-wave $\delta^{\vphantom{1}}{\bar{\delta}}$ states, we identify
$V_{\delta^{\vphantom{1}}{\bar{\delta}}}(\mathbf{r})$ as the uncoupled
$\Sigma^{+}_{g}$ potential Lebed (2017), with parameters that are calculated
on the lattice Berwein  _et al._ (2015). Thus, we parametrize

$V_{\delta^{\vphantom{1}}{\bar{\delta}}}(r)=-\frac{\alpha}{r}+\sigma
r+V_{0}+m_{\delta^{\vphantom{1}}}+m_{{\bar{\delta}}}\,,$ (15)

where $\alpha,\sigma,$ and $V_{0}$ are
$0.053~{}\text{GeV$\cdot$fm},1.097~{}\text{GeV/fm},$ and
$-0.180~{}\text{GeV}$, respectively. For each flavor sector, one also requires
values for the $\delta^{\vphantom{1}},{\bar{\delta}}$ masses, and a
corresponding list of di-hadron thresholds with matching $J^{PC}$ quantum
numbers must be identified. One then needs only to determine the appropriate
form for the mixing potentials. Reference Bruschini and González (2020) argues
for a Gaussian form, which we also adopt here. Explicitly,

$|V_{\rm
mix}^{(i)}(r)|=\frac{\Delta}{2}\exp\\!\left\\{-\frac{1}{2}\frac{\left[V^{\vphantom{\dagger}}_{\delta^{\vphantom{1}}{\bar{\delta}}}(r)-T_{M_{1}\overline{M}_{2}}^{(i)}\right]^{2}}{(\rho\sigma)^{2}}\right\\},$
(16)

where $\sigma$ is the same string-tension parameter as in the Cornell-like
potential of Eq. (15). $\Delta$ is a free parameter with units of energy
indicating the strength of the mixing, and $\rho$ is the radial scale for the
level crossing Bruschini and González (2020).

It is a useful first exercise to consider this procedure for a single
threshold. The diabatic-potential matrix then reads

$\text{V}(r)=\begin{pmatrix}V^{\vphantom{\dagger}}_{\delta^{\vphantom{1}}{\bar{\delta}}}(r)&V^{\vphantom{\dagger}}_{\rm
mix}(r)\\\ V^{\vphantom{\dagger}}_{\rm
mix}(r)&T_{M_{1}\overline{M}_{2}}\end{pmatrix},$ (17)

where we have replaced $\mathbf{r}$ with its magnitude, since all relevant
potentials at this stage depend solely upon the diquark separation distance.
The corresponding single-threshold kinetic-energy operator becomes

$\text{K}=\begin{pmatrix}-\frac{\hbar^{2}}{2\mu_{\delta^{\vphantom{1}}{\bar{\delta}}}}&\\\
&-\frac{\hbar^{2}}{2\mu_{M_{1}\overline{M}_{2}}}\end{pmatrix}\nabla^{2}.$ (18)

The diabatic-potential matrix eigenvalues [denoted $V_{-}(r)$ and $V_{+}(r)$]
and corresponding eigenvectors [$|\xi_{-}(r)\rangle$ and $|\xi_{+}(r)\rangle$]
may be directly related to the $\delta^{\vphantom{1}}{\bar{\delta}}$ and di-
meson threshold light-field eigenstates via a generic transformation matrix:

$\text{R}(r)=\begin{pmatrix}\ \ \cos\theta(r)&\sin\theta(r)\\\
-\sin\theta(r)&\cos\theta(r)\end{pmatrix},$ (19)

such that

$\text{R}(r)\text{V}(r)\text{R}^{\dagger}(r)=\text{diag}\left\\{V_{-}(r),V_{+}(r)\right\\}\,.$
(20)

Here, $\theta$ is identified as the mixing angle between the
$\delta^{\vphantom{1}}{\bar{\delta}}$ and meson-meson terms. Since
$\text{R}(r)$ diagonalizes $\text{V}(r)$, one can deduce an expression for
$\theta$ in terms of its matrix elements Bruschini and González (2020):

$\theta(r)=\frac{1}{2}\arctan\left(\frac{2V_{\rm
mix}(r)}{T_{M_{1}\overline{M}_{2}}-V^{\vphantom{\dagger}}_{\delta^{\vphantom{1}}{\bar{\delta}}}(r)}\right).$
(21)

As an example, we plot $\theta$ for
$\delta^{\vphantom{1}}({\bar{\delta}})=cq\,(\overline{cq}^{\prime})$ and
$M_{1}\overline{M}_{2}=D\bar{D}^{*}$ in Fig. 1 using phenomenologically viable
values of $\Delta$ and $\rho$ [see Eq. (22)], as obtained below.

Table 1: Calculated eigenvalues and component-state admixtures for the ${c\bar{c}q\bar{q}^{\prime}}$ sector obtained from solving Eq. (10) for specific $J^{PC}$ numbers. Also presented are the expectation values $\langle r\rangle$ of the radial coordinate $r$ (corresponding to $\delta$-${\bar{\delta}}$ separation) as well as $\langle r^{2}\rangle^{1/2}$ for each state. Suppressed entries indicate contributions that are individually $<\\!1\%$. $J^{PC}$ | $E\ ({\rm MeV})$ | $\delta^{\vphantom{1}}{\bar{\delta}}$ | $D\bar{D}^{*}$ | $D_{s}\bar{D}_{s}$ | $D^{*}\bar{D}^{*}$ | $D_{s}^{*}\bar{D}_{s}^{*}$ | $\langle r\rangle({\rm fm})$ | $\langle r^{2}\rangle^{1/2}({\rm fm})$
---|---|---|---|---|---|---|---|---
$0^{++}$ | 3905.4 | 63.0% |  | 27.4% | 8.4% | 1.2% | 0.596 | 0.605
$1^{++}$ | 3871.5 | 8.6% | 91.4% |  |  |  | 4.974 | 5.459
$2^{++}$ | 3922.3 | 83.1% |  | 1.5% | 13.9% | 1.5% | 0.443 | 0.497

##  IV Results

###  IV.1 ${c\bar{c}q\bar{q}^{\prime}}$ Sector

In order to obtain values for the $\Delta$ and $\rho$ parameters in the mixing
potentials, we fit them such that the lowest $1^{++}$ eigenvalue in the
${c\bar{c}q\bar{q}^{\prime}}$ sector is near the measured value
$m_{X(3872)}\\!=\\!3871.7$ MeV. Of course, a range of $\Delta,\rho$ value
pairs achieves this result. Following the approach of Ref. Bruschini and
González (2020), we select a pair that produces the correct asymptotic
behavior for the single-threshold mixing angle $\theta(r)$, as modeled in Fig.
1. Specifically, armed with unquenched lattice-QCD calculations Prkacin  _et
al._ (2006), we require that $\theta\\!\to\\!\pi/2$ rapidly for
$r\\!>\\!r_{c}$, the value at which the $\delta^{\vphantom{1}}{\bar{\delta}}$
Cornell-like potential Eq. (15) crosses the $M_{1}\overline{M}_{2}$ threshold
Eq. (14). According to Eq. (21), the $\Delta,\rho$ pair must also produce the
nonzero value $\theta\\!=\\!\pi/4$ at $r\\!=\\!r_{c}$, indicating a nontrivial
mixing of the $\delta^{\vphantom{1}}{\bar{\delta}}$ and meson-meson
components. The specific values used here are

$\Delta_{\delta^{\vphantom{1}}{\bar{\delta}}}=0.300\>\text{GeV}\,,\;\rho_{\delta^{\vphantom{1}}{\bar{\delta}}}=0.185\>\text{fm}\,.$
(22)

Figure 1: The mixing angle $\theta(r)$ [Eq. (21), using Eqs. (16) and (22)]
between $\delta^{\vphantom{1}}{\bar{\delta}}=c\bar{c}qq^{\prime}$ and
$D\bar{D}^{*}$ in the $1^{++}$ channel, as a function of heavy-source radial
separation $r$. The dashed line shows the critical radius $r_{c}$, i.e., the
separation at which $V_{\delta{\bar{\delta}}}(r)=T_{D\bar{D}^{*}}$, using Eqs.
(14)–(15) and (23).

Again, these values are not unique, but rather serve as proof of principle for
the approach. The diquark mass is then the only remaining free parameter, with
its value taken as reported in Ref. Giron  _et al._ (2021):

$m_{\delta^{\vphantom{1}}=(cq)}=m_{{\bar{\delta}}=(\bar{c}\bar{q})}=1.9271\;\text{GeV}.$
(23)

One may then numerically solve Eq. (10) as a coupled set of equations. We
follow the procedure described in Sec. IV.F of Ref. Bruschini and González
(2020), solving for the lowest bound states. The resulting eigenvalues and
their corresponding state mixtures are collected in Table 1. Note the omission
of a $J^{PC}=1^{+-}$ entry in Table 1 as compared with Eqs. (2) or (4); while
the diabatic formalism allows for the formation of $C\\!=\\!-1$ eigenstates,
the di-hadron thresholds alone provide no mechanism to lift the degeneracy
with the $1^{++}$ eigenstate.111In the full dynamical diquark model, this
degeneracy is lifted by operators that couple to individual diquark spins
Giron  _et al._ (2020); Giron and Lebed (2020b). Notably, we find that the
$1^{++}$ eigenstate, which we associate with $X(3872)$, consists primarily of
$D\bar{D}^{*}$ content, consistent with expectations given its proximity to
that threshold. In contrast, we find that the other eigenstates ($0^{++}$ and
$2^{++}$) are primarily of $\delta^{\vphantom{1}}{\bar{\delta}}$ content.

Using the results of this formalism, it is also straightforward to calculate
certain transition rates, which of course provide numerous predictions for
comparison with experimental results. Here, we focus specifically on the
radiative transition of the $J^{PC}=1^{++}$ eigenstate to $J/\psi$ and
$\psi(2S)$. For E1 and M1 transitions of states within the dynamical diquark
model, substantial work has already been performed in Ref. Gens  _et al._
(2021). There, a standard equation for E1 partial widths for the process
$i\\!\to\\!\gamma f$ is adapted to the case of exotic-to-exotic transitions,
and it may be recast for the present case as

$\Gamma_{\rm E1}\left(n^{2s_{Q\bar{Q}}+1}(J_{q\bar{q}})_{J}\rightarrow
n^{\prime{2s+1}}(L^{\prime})_{J^{\prime}}+\gamma\right)=\frac{4}{3}C_{fi}\,\delta_{s^{\vphantom{\prime}}_{Q\bar{Q}}s^{\prime}_{Q\bar{Q}}}Q^{2}_{\delta}\alpha\,|\langle\psi_{f}|r|\psi_{i}\rangle|^{2}E^{3}_{\gamma}\frac{E_{f}^{Q\bar{Q}}}{M_{i}^{Q\bar{Q}q\bar{q}}}\,,$
(24)

where

$C_{fi}\equiv{\rm
max}\left(J_{q\bar{q}},L^{\prime}\right)(2J^{\prime}+1)\begin{Bmatrix}L^{\prime}&J^{\prime}&s_{Q\bar{Q}}\\\
J&J_{q\bar{q}}&1\end{Bmatrix}^{2}\,.$ (25)

In this expression, $\alpha$ is the QED fine-structure constant, $E_{\gamma}$
is the photon energy (measured in the rest frame of the decaying exotic
initial state $Q\bar{Q}q\bar{q}$ of mass $M_{i}^{Q\bar{Q}q\bar{q}}$, in which
the conventional quarkonium final state $Q\bar{Q}$ recoils with energy
$E_{f}^{Q\bar{Q}}$), and $Q_{\delta^{\vphantom{1}}}\\!=\\!Q_{Q}\\!+Q_{q}$ is
the diquark charge in units of proton charge (meaning that the diquark is
treated as a single quasiparticle coupling to the photon). Furthermore,
$s^{(\prime)}_{Q\bar{Q}}$ denotes the total heavy-quark spin eigenvalue for
the initial (final) state, $n^{(\prime)}$ indicates the principal quantum
number for the initial (final) state, $J_{q\bar{q}}$ ($L^{\prime}$) represents
the total angular momentum carried by the light d.o.f. in the initial (final)
state, and $J$ ($J^{\prime}$) is the total angular momentum eigenvalue of the
initial (final) state, respectively.

Before proceeding, we hasten to emphasize that Eq. (24) refers to E1
transitions between two states with the same two fundamental charged
constituents, which together generate an electric dipole moment. This
expression is reasonable (and indeed is commonly used) for conventional
quarkonium, albeit often including relativistic corrections that are outlined
in Ref. Gens  _et al._ (2021). This reference applied Eq. (24) only to the
case of transitions between two $\delta^{\vphantom{1}}{\bar{\delta}}$ states,
again treating the diquarks as single compact quasiparticles. Applying it to
$\delta^{\vphantom{1}}{\bar{\delta}}\\!\to Q\bar{Q}$ transitions is much more
questionable since, among other possible objections, such a transition
requires the annihilation of the $q\bar{q}$ pair. At minimum, Eq. (24) must at
least be modified to accommodate the Okubo-Zweig-Iizuka (OZI) suppressed
amplitude $q\bar{q}\\!\to\\!g$ in this process. While a proper treatment of
the issue lies outside the scope of this work, we may at least note that the
magnitude $\epsilon$ of the OZI-suppression is expected to be significantly
$<\\!1$, and so $\Gamma_{\rm E1}$ of Eq. (24) is diminished by a factor
$\epsilon^{2}$. One may expect $\epsilon$ to depend upon both the radial
excitation number $n$ of the initial state, as well as the $q\bar{q}$ spin
state.

The most sensitive term to calculate in Eq. (24) is the overlap
$|\langle\psi_{f}|r|\psi_{i}\rangle|^{2}$ of the initial and final radial wave
functions weighted by $r$, the characteristic spatial separation of the heavy
sources. Recall that the full exotic eigenstate in this analysis is taken to
be a mixture of both $\delta^{\vphantom{1}}{\bar{\delta}}$ and
$M_{1}\overline{M}_{2}$ components. However, the expected relative separation
of the mesons in the $D^{0}\bar{D}^{*0}$ component of $X(3872)$ is expected to
be much greater than that of a pure $\delta^{\vphantom{1}}{\bar{\delta}}$
state. Explicitly, using the results of Table 1, we find that the
$D^{0}\bar{D}^{*0}$ component alone should have $\langle r\rangle\\!=\\!5.092$
fm, in agreement with the crude estimate provided by the Bohr-radius analogue
$r\\!\sim\\!\frac{1}{\sqrt{2\mu E_{\rm bind.}}}$. Conversely, the value of
$\langle r\rangle$ for the pure
$\delta^{\vphantom{1}}({\bar{\delta}})=cq\,(\bar{c}\bar{q}^{\prime})$
component is only 0.361 fm. Since the decays $X(3872)\\!\to\\!\gamma\psi$
requires the annihilation of the light $q\bar{q}$ pair in $X(3872)$, one
naively expects that the $\delta^{\vphantom{1}}{\bar{\delta}}$ component
should dominate in this process. Let us test this expectation.

Once the fractional $\delta^{\vphantom{1}}{\bar{\delta}}$ content in the
$1^{++}$ eigenstate is known, one may calculate the E1 partial width from Eq.
(24). Using the potential of Eq. (15) and numerical inputs for its parameters
given there, and adopting $m_{c}\\!=\\!1840$ MeV Bruschini and González (2020)
in place of $m_{\delta^{\vphantom{1}}}$, we compute the radial wave functions
of $J/\psi$, $\psi(2S)$, and $\chi_{c1}(2P)$ treated as pure conventional
charmonium states. We also solve for the $(1S)$
$\delta^{\vphantom{1}}{\bar{\delta}}=cq\,(\bar{c}\bar{q})$ wave function,
which was an essential ingredient in Ref. Gens  _et al._ (2021). Introducing
these values into Eq. (24), we find

$\displaystyle\Gamma_{\text{E}1}[X(3872)\rightarrow\gamma J/\psi]$
$\displaystyle=$ $\displaystyle 469\,\epsilon^{2}\,\text{keV},$
$\displaystyle\Gamma_{\text{E}1}[X(3872)\rightarrow\gamma\,\psi(2S)]$
$\displaystyle=$ $\displaystyle 1.56\,\epsilon^{2}\,\text{keV},$ (26)

while

$\displaystyle\Gamma_{\text{E}1}[\chi_{c1}(2P)\rightarrow\gamma J/\psi]$
$\displaystyle=$ $\displaystyle{\rm 20.8}\,\text{keV},$
$\displaystyle\Gamma_{\text{E}1}[\chi_{c1}(2P)\rightarrow\gamma\,\psi(2S)]$
$\displaystyle=$ $\displaystyle 78.0\,\text{keV}.$ (27)

We note that the results of Eqs. (IV.1) are comparable to the predictions of
71 keV and 95 keV, respectively, from Ref. Swanson (2004) [which also uses Eq.
(24), but with somewhat different numerical values for the matrix elements].

On the other hand, Eq. (IV.1) predicts (setting $\epsilon\\!\to\\!1$) wildly
different $X(3872)$ radiative widths than the current PDG averages Zyla  _et
al._ (2020), $10.1\\!\pm\\!4.7$ keV for $J/\psi$ and $54\pm 25$ keV for
$\psi(2S)$. The PDG separately gives an average for their ratio:

$R\equiv\frac{\Gamma[X(3872)\to\gamma\psi(2S)]}{\Gamma[X(3872)\to\gamma
J/\psi]}=2.6\pm 0.6\,,$ (28)

while the corresponding ratio from Eqs. (IV.1), taking the two values of
$\epsilon$ equal, is only $3.3\times 10^{-3}$. Interestingly, this number
approximately equals an early estimate for $R$ based upon a molecular-model
calculation using vector-meson dominance Swanson (2004). However, a modern
molecular-model calculation employing an effective Lagrangian supports values
of $R$ in a typical range of order several tenths Guo  _et al._ (2015), and
these predictions are numerically stable when a small admixture of the
conventional-charmonium state $\chi_{c1}(2P)$ is included Cincioglu and
Ozpineci (2019).

The origin of such a small predicted value for $R$ in the dynamical diquark
model can be traced to the fact that the $\delta^{\vphantom{1}}{\bar{\delta}}$
component of $X(3872)$ is a ground-state ($n\\!=\\!1$) radial mode like
$J/\psi$, while $\chi_{c1}(2P)$ is a first-excited ($n\\!=\\!2$) radial mode
like $\psi(2S)$. Indeed, our calculated idealized
$\delta^{\vphantom{1}}{\bar{\delta}}$ wave function has almost exactly the
same size as our calculated $c\bar{c}$ $J/\psi$ wave function, leading to
almost complete overlap and a large transition to $J/\psi$. In a more
realistic treatment including a finite diquark $(cq)$ size Giron  _et al._
(2020), one expects this overlap enhancement to be muted. Nevertheless, one
still expects the phase-space advantage of $J/\psi$ over $\psi(2S)$ to
generate a rather small $R$ value in this model. If a value of $R$
substantially larger than 1 persists in the data, we conclude that one must
have some underlying preferential coupling to $\psi(2S)$ over $J/\psi$, either
from a significant $\chi_{c1}(2P)$ component in $X(3872)$, from an enhanced
ratio of $\psi(2S)$ to $J/\psi$ effective-theory couplings as suggested in
Ref. Guo  _et al._ (2015), or from an enhanced wave-function overlap of
$\delta^{\vphantom{1}}{\bar{\delta}}$ to the spatially larger $\psi(2S)$
state.

The overall size of the prediction for $\Gamma[X(3872)\\!\to\\!\gamma J/\psi]$
from the first of Eqs. (IV.1) also points to the necessity of $\epsilon$ being
no larger than about $1/6$, which is not an unreasonable expectation for an
OZI-suppressed amplitude.

Lastly, in a model for $X(3872)$ with not just
$\delta^{\vphantom{1}}{\bar{\delta}}$ but also $D^{0}\bar{D}^{*0}$ components
[and possibly $\chi_{c1}(2P)$ as well], the full value for its radiative decay
widths can include contributions from more than one source. In that case, it
becomes crucial to determine the relative phase of each contribution, since
interference may be critical to obtaining a physically accurate result. This
effect in a mixed $D^{0}\bar{D}^{*0}$-$\chi_{c1}(2P)$ model was considered in
Ref. Cincioglu and Ozpineci (2019).

###  IV.2 ${c\bar{c}s\bar{s}}$ and $c\bar{c}q\bar{s}$ Sectors

In the ${c\bar{c}s\bar{s}}$ sector, one encounters an additional free
parameter, the $(cs)$ diquark mass. We fix this mass so that the $0^{++}$
energy eigenvalue matches the $X(3915)$ mass Zyla  _et al._ (2020). This
choice produces the result

$m_{\delta^{\vphantom{1}}=(cs)}=m_{{\bar{\delta}}=(\bar{c}\bar{s})}=1.9450\;\text{GeV}.$
(29)

Note that this value only slightly exceeds $m_{\delta^{\vphantom{1}}=(cq)}$
given in Eq. (23), despite containing a heavier $s$ quark; in particular, it
is substantially smaller than the value (2.080 GeV) obtained from the analysis
of Ref. Giron  _et al._ (2021), because that work (unlike here) incorporates
spin-splitting fine structure for the $\Sigma^{+}_{g}(1S)$
$\delta^{\vphantom{1}}{\bar{\delta}}$ multiplet. The remaining eigenstates may
then be calculated, and are presented in Table 2. Remarkably, the $1^{++}$ and
$2^{++}$ states exhibit rather little $D^{(*)}_{(s)}\bar{D}^{*}_{(s)}$
content, while the $0^{++}$ state shows significant mixing with
$D_{s}\overline{D}_{s}$.

Table 2: The same as in Table 1, for the ${c\bar{c}s\bar{s}}$ sector. $J^{PC}$ | $E\ ({\rm MeV})$ | $\delta^{\vphantom{1}}{\bar{\delta}}$ | $D_{s}\bar{D}_{s}$ | $D^{*}\\!\bar{D}^{*}$ | $D_{s}\bar{D}_{s}^{*}$ | $D_{s}^{*}\bar{D}_{s}^{*}$ | $\langle r\rangle({\rm fm})$ | $\langle r^{2}\rangle^{1/2}({\rm fm})$
---|---|---|---|---|---|---|---|---
$0^{++}$ | 3922.0 | 50.6% | 39.9% | 8.3% |  | 1.2% | 0.627 | 0.757
$1^{++}$ | 3978.7 | 86.9% |  | 1.5% | 11.0% |  | 0.421 | 0.470
$2^{++}$ | 3998.7 | 95.2% |  |  | 1.5% | 3.3% | 0.394 | 0.437
Table 3: The same as in Tables 1 & 2, for the $c\bar{c}q\bar{s}$ sector. $J^{P}$ | $E\ ({\rm MeV})$ | $\delta^{\vphantom{1}}{\bar{\delta}}$ | $D^{*}\\!\bar{D}_{s}$ | $D\bar{D}_{s}^{*}$ | $D^{*}\\!\bar{D}_{s}^{*}$ | $\langle r\rangle({\rm fm})$ | $\langle r^{2}\rangle^{1/2}({\rm fm})$
---|---|---|---|---|---|---|---
$0^{+}$ | 3983.5 | 93.2% |  |  | 6.8% | 0.397 | 0.440
$1^{+}$ | 3914.5 | 66.7% | 16.5% | 15.9% |  | 0.482 | 0.548
$2^{+}$ | 3962.1 | 90.4% | 1.8% | 1.8% | 6.0% | 0.411 | 0.460

The calculations for $c\bar{c}q\bar{s}$ states follow entirely from the
previously determined parameters, and are presented separately in Table 3.
This sector in particular contains a unique threshold structure, in that the
$D^{*}\\!\bar{D}_{s}$ and $D\bar{D}^{*}_{s}$ thresholds differ by less than 4
MeV, for each fixed value of total charge for the meson pair. Our calculations
show that thresholds extremely close to each other tend to evenly share state
content (assuming that the allowed $\ell$ quantum numbers are the same).

In all three sectors, whenever
$\delta^{\vphantom{1}}{\bar{\delta}}\,(\ell\\!=\\!0)$ content dominates, we
furthermore find a preferential coupling to $S$-wave meson-meson combinations
over $D$-wave combinations, even if the latter threshold is closer to the mass
eigenvalue. This result is exhibited by the
${c\bar{c}q\bar{q}^{\prime}}~{}2^{++}$, ${c\bar{c}s\bar{s}}~{}1^{++}$, and
$c\bar{c}q\bar{s}~{}2^{+}$ states. The ${c\bar{c}q\bar{q}^{\prime}}~{}2^{++}$
state is the most extreme of these, preferring to couple more strongly to
$D^{*}_{s}\bar{D}^{*}_{s}$ over $D_{s}\bar{D}_{s}$, despite the $\sim$280 MeV
difference. This behavior is further exhibited by the results for individual
thresholds; for example, in the $2^{++}{c\bar{c}q\bar{q}^{\prime}}$ state, the
distinct $D^{*}\bar{D}^{*}$ channels with $\ell=0$ and $\ell=2$ contribute
$11.6\%$ and $1.2\%$, respectively. The expected suppression associated with
the near-orthogonal overlap of an $S$-wave
$\delta^{\vphantom{1}}{\bar{\delta}}$ state and a $D$-wave meson-meson state
provides a natural explanation.

We also note that the expectation values $\left<r\right>$ (or $\langle
r^{2}\rangle^{1/2}$) for each state increase with increasing meson-meson
content, consistent with predictions of $\left<r\right>$ for hadron molecular
states Liu  _et al._ (2008). As discussed above, the
$1^{++}{c\bar{c}q\bar{q}^{\prime}}$ state exhibits a several-fm value for
$\left<r\right>$, as is expected for a pure $D\bar{D}^{*}$ molecular state
with a binding energy of tenths of an MeV Liu  _et al._ (2008).

##  V Conclusions

The introduction into a diquark-based model—in this paper, the variant known
as the dynamical diquark model—of effects caused by the proximity of some of
its eigenstates to di-hadron thresholds creates a model for multiquark exotic
hadrons that combines the best features of both diquark and hadron-molecular
models. One universal model can then incorporate exotic candidates that lie
quite close to such thresholds, but which still maintain a close physical
connection to other exotic candidates with no obvious di-hadron
interpretation.

Coupled-channel calculations between bound states and thresholds have of
course been carried out many times in the past, originally in the context of
atomic and molecular physics, but more recently for exotic-hadron candidates
as well. The so-called diabatic approach used here earns distinction as the
natural (and rigorous) generalization of the Born-Oppenheimer approximation
used to compute the unperturbed spectrum of heavy-heavy systems such as
quarkonium and its hybrid excitations, and in this work for the first time,
the diquark-antidiquark states of the dynamical diquark model.

The precise functional form for the mixing potential between unperturbed
states and thresholds used here is phenomenological in origin, but it is
motivated by lattice-QCD results involving the breaking of the color flux tube
between the heavy sources. Clearly, advances in lattice simulations can be
incorporated into improved future calculations.

We find, using parameters from the most recent analysis in the dynamical
diquark model, that the famous $X(3872)$ can naturally contain a dominant
component of $D^{0}\bar{D}^{*0}$ and yet originate as the unique isosinglet
$1^{++}$ member of the lightest multiplet of hidden-charm diquark-antidiquark
states. Masses of all the remaining members of the multiplet are then
predicted, along with their di-hadron content generated by nearby thresholds.
Since the original model has been extended to study hidden-charm, hidden-
strangeness and open-strange states as well, we also present calculations for
those flavor sectors. General results include the observation that $S$-wave
thresholds always produce larger effects than $D$-wave thresholds, even if the
latter are substantially closer to the mass eigenvalue, and that coincident
di-hadron thresholds with different arrangements of the light flavors in the
hadrons have comparable effects.

We also use our results on state content to calculate the radiative decay
widths for $X(3872)\\!\to\\!\gamma J/\psi$ and $\gamma\psi(2S)$, which
provides crucial information on the short-distance components of $X(3872)$,
since its di-hadron component is spatially much larger than that of its
diquark-antidiquark component. Interestingly, we obtain a $\gamma\psi(2S)$
width much smaller than the current measured value, due not just to the fact
that the compact component amounts to only about 10% of the full state, but
also that the overlap of the $1S$ diquark-antidiquark and $2S$ $c\bar{c}$
states is small. We conclude that obtaining the full radiative width may
require including a comparable contribution from the diffuse di-hadron
component, or even from a $\chi_{c1}(2P)$ component (which can also be
included in the diabatic formalism).

These initial calculations, while quite encouraging, remain quite incomplete.
First, no spin- or flavor-dependent effects (besides explicit differences in
the diquark mass) have been incorporated into the unperturbed side of this
calculation. The Hamiltonian for the full dynamical diquark model contains
spin-spin and spin-isospin operators that, for example, distinguish masses of
$X(3872)$, $Z_{c}(3900)$, and $Z_{c}(4020)$. In this calculation, the only
mass splittings arise from the explicit differences of $D^{(*)}_{(s)}$ masses.
One major future research thrust is the inclusion of explicit fine-structure
effects into the initial unperturbed states.

A second possible improvement involves the treatment of the threshold
contribution. Di-hadron molecular models have been alluded to several times in
this paper, but the actual treatment of the di-hadron state uncoupled from the
diquark-antidiquark state, according to Eq. (14), is simply that of two free
hadrons. Should it be desirable to regard the pair as forming a true di-hadron
molecule, then directly replacing Eq. (14) with a mildly attractive potential
would be straightforward in this formalism.

Finally, the diabatic formalism as presented here strictly only applies to
states either below or not too far above significant di-hadron thresholds.
Some of the exotic candidates [e.g., $Z_{c}(4430)$] lie rather far from
relevant thresholds, and consequently have large decay widths. Such broad
resonances should properly be treated as poles in di-hadron scattering
amplitudes. However, the diabatic formalism has been developed, in the case of
mixing with conventional quarkonium, to include the calculation of strong
decay widths Bruschini and González (2021a) and di-hadron scattering
amplitudes Bruschini and González (2021b). These methods can be immediately
adapted to the case of mixing with diquark-antidiquark states, and will also
be incorporated into future work.

###### Acknowledgements.

This work was supported by the National Science Foundation (NSF) under Grants
No. PHY-1803912 and PHY-2110278.

## References

  * Choi  _et al._ (2003) S. Choi _et al._ (Belle Collaboration), Phys. Rev. Lett. 91, 262001 (2003), arXiv:hep-ex/0309032 .
  * Lebed  _et al._ (2017) R. Lebed, R. Mitchell, and E. Swanson, Prog. Part. Nucl. Phys. 93, 143 (2017), arXiv:1610.04528 [hep-ph] .
  * Chen  _et al._ (2016) H.-X. Chen, W. Chen, X. Liu, and S.-L. Zhu, Phys. Rep. 639, 1 (2016), arXiv:1601.02092 [hep-ph] .
  * Hosaka  _et al._ (2016) A. Hosaka, T. Iijima, K. Miyabayashi, Y. Sakai, and S. Yasui, Prog. Theor. Exp. Phys. 2016, 062C01 (2016), arXiv:1603.09229 [hep-ph] .
  * Esposito  _et al._ (2017) A. Esposito, A. Pilloni, and A. Polosa, Phys. Rep. 668, 1 (2017), arXiv:1611.07920 [hep-ph] .
  * Guo  _et al._ (2018) F.-K. Guo, C. Hanhart, U.-G. Meißner, Q. Wang, Q. Zhao, and B.-S. Zou, Rev. Mod. Phys. 90, 015004 (2018), arXiv:1705.00141 [hep-ph] .
  * Ali  _et al._ (2017) A. Ali, J. Lange, and S. Stone, Prog. Part. Nucl. Phys. 97, 123 (2017), arXiv:1706.00610 [hep-ph] .
  * Olsen  _et al._ (2018) S. Olsen, T. Skwarnicki, and D. Zieminska, Rev. Mod. Phys. 90, 015003 (2018), arXiv:1708.04012 [hep-ph] .
  * Karliner  _et al._ (2018) M. Karliner, J. Rosner, and T. Skwarnicki, Annu. Rev. Nucl. Part. Sci. 68, 17 (2018), arXiv:1711.10626 [hep-ph] .
  * Yuan (2018) C.-Z. Yuan, Int. J. Mod. Phys. A 33, 1830018 (2018), arXiv:1808.01570 [hep-ex] .
  * Liu  _et al._ (2019) Y.-R. Liu, H.-X. Chen, W. Chen, X. Liu, and S.-L. Zhu, Prog. Part. Nucl. Phys. 107, 237 (2019), arXiv:1903.11976 [hep-ph] .
  * Brambilla  _et al._ (2020) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C.-P. Shen, C. Thomas, A. Vairo, and C.-Z. Yuan, Phys. Rept. 873, 1 (2020), arXiv:1907.07583 [hep-ex] .
  * Chen  _et al._ (2022) H.-X. Chen, W. Chen, X. Liu, Y.-R. Liu, and S.-L. Zhu, (2022), arXiv:2204.02649 [hep-ph] .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Braaten and Kusunoki (2004) E. Braaten and M. Kusunoki, Phys. Rev. D 69, 074005 (2004), arXiv:hep-ph/0311147 .
  * Tiesinga  _et al._ (2021) E. Tiesinga, P. Mohr, D. Newell, and B. Taylor, Rev. Mod. Phys. 93, 025010 (2021).
  * Barnes  _et al._ (2005) T. Barnes, S. Godfrey, and E. Swanson, Phys. Rev. D 72, 054026 (2005), arXiv:hep-ph/0505002 [hep-ph] .
  * Suzuki (2005) M. Suzuki, Phys. Rev. D 72, 114013 (2005), arXiv:hep-ph/0508258 .
  * Brodsky  _et al._ (2014) S. Brodsky, D. Hwang, and R. Lebed, Phys. Rev. Lett. 113, 112001 (2014), arXiv:1406.7281 [hep-ph] .
  * Lebed (2017) R. Lebed, Phys. Rev. D 96, 116003 (2017), arXiv:1709.06097 [hep-ph] .
  * Giron  _et al._ (2019) J. Giron, R. Lebed, and C. Peterson, J. High Energy Phys. 05, 061 (2019), arXiv:1903.04551 [hep-ph] .
  * Giron  _et al._ (2020) J. Giron, R. Lebed, and C. Peterson, J. High Energy Phys. 01, 124 (2020), arXiv:1907.08546 [hep-ph] .
  * Giron and Lebed (2020a) J. Giron and R. Lebed, Phys. Rev. D 101, 074032 (2020a), arXiv:2003.02802 [hep-ph] .
  * Giron and Lebed (2020b) J. Giron and R. Lebed, Phys. Rev. D 102, 014036 (2020b), arXiv:2005.07100 [hep-ph] .
  * Giron  _et al._ (2021) J. Giron, R. Lebed, and S. Martinez, Phys. Rev. D 104, 054001 (2021), arXiv:2106.05883 [hep-ph] .
  * Giron and Lebed (2021) J. Giron and R. Lebed, Phys. Rev. D 104, 114028 (2021), arXiv:2110.05557 [hep-ph] .
  * Lebed (2015) R. Lebed, Phys. Lett. B 749, 454 (2015), arXiv:1507.05867 [hep-ph] .
  * Born and Oppenheimer (1927) M. Born and R. Oppenheimer, Ann. der Phys. 389, 457 (1927).
  * Baer (2006) M. Baer,  _Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections_ (Wiley, New Jersey, 2006).
  * Bruschini and González (2020) R. Bruschini and P. González, Phys. Rev. D 102, 074002 (2020), arXiv:2007.07693 [hep-ph] .
  * Bruschini and González (2021a) R. Bruschini and P. González, Phys. Rev. D 103, 074009 (2021a), arXiv:2101.04636 [hep-ph] .
  * Giron and Lebed (2020c) J. Giron and R. Lebed, Phys. Rev. D 102, 074003 (2020c), arXiv:2008.01631 [hep-ph] .
  * Bulava  _et al._ (2019) J. Bulava, B. Hörz, F. Knechtli, V. Koch, G. Moir, C. Morningstar, and M. Peardon, Phys. Lett. B 793, 493 (2019), arXiv:1902.04006 [hep-lat] .
  * Berwein  _et al._ (2015) M. Berwein, N. Brambilla, J. Tarrús Castellà, and A. Vairo, Phys. Rev. D 92, 114019 (2015), arXiv:1510.04299 [hep-ph] .
  * Prkacin  _et al._ (2006) Z. Prkacin, G. Bali, T. Dussel, T. Lippert, H. Neff, and K. Schilling, PoS LAT2005, 308 (2006), arXiv:hep-lat/0510051 .
  * Gens  _et al._ (2021) J. Gens, J. Giron, and R. Lebed, Phys. Rev. D 103, 094024 (2021), arXiv:2102.12591 [hep-ph] .
  * Swanson (2004) E. Swanson, Phys. Lett. B 598, 197 (2004), arXiv:hep-ph/0406080 .
  * Guo  _et al._ (2015) F.-K. Guo, C. Hanhart, Y. S. Kalashnikova, U.-G. Meißner, and A. Nefediev, Phys. Lett. B 742, 394 (2015), arXiv:1410.6712 [hep-ph] .
  * Cincioglu and Ozpineci (2019) E. Cincioglu and A. Ozpineci, Phys. Lett. B 797, 134856 (2019), arXiv:1901.03138 [hep-ph] .
  * Liu  _et al._ (2008) Y.-R. Liu, X. Liu, W.-Z. Deng, and S.-L. Zhu, Eur. Phys. J. C 56, 63 (2008), arXiv:0801.3540 [hep-ph] .
  * Bruschini and González (2021b) R. Bruschini and P. González, Phys. Rev. D 104, 074025 (2021b), arXiv:2107.05459 [hep-ph] .


# Solvable model of driven matter with pinning

Gourab Kumar Sar<EMAIL_ADDRESS>Physics and Applied Mathematics Unit,
Indian Statistical Institute, 203 B. T. Road, Kolkata 700108, India  Dibakar
Ghosh<EMAIL_ADDRESS>Physics and Applied Mathematics Unit, Indian
Statistical Institute, 203 B. T. Road, Kolkata 700108, India  Kevin O’Keeffe
<EMAIL_ADDRESS>Senseable City Lab, Massachusetts Institute of
Technology, Cambridge, MA, USA, 02139

###### Abstract

(Received XX MONTH XX; accepted XX MONTH XX; published XX MONTH XX)

We present a simple model of driven matter in a 1D medium with pinning
impurities, applicable to magnetic domains walls, confined colloids, and other
systems. We find rich dynamics, including hysteresis, reentrance,
quasiperiodicity, and two distinct routes to chaos. In contrast to other
minimal models of driven matter, the model is solvable: we derive the full
phase diagram for small $N$, and for large $N$, derive expressions for order
parameters and several bifurcation curves. The model is also realistic. Its
collective states match those seen in the experiments of magnetic domain
walls, and its force-velocity curve imitates those of superconductor vortices.

DOI: XXXXXXX

Driving matter through disordered environments has diverse applications in
science. Magnetic domain walls and other quasi-particles may be driven off
material defects and used as memory units in spintronics Allwood et al.
(2005); Luo et al. (2020); Wiesendanger (2016); Reichhardt et al. (2022).
Electromagnetic colloids may be forced to self-assemble into cargo carriers
for high precision medicine Reddy et al. (2012); Yang et al. (2017); Bricard
et al. (2015); Caleap and Drinkwater (2014). The colloids may also be used to
repair circuits Li et al. (2015), purify water Pinto et al. (2020), and
shatter blood clots Cheng et al. (2014); Manamanchaiyaporn et al. (2021).

All these applications rely on our ability to predict how a given matter
collective reacts to driving. To be concrete here, imagine a swarm of
particles being pushed around by an external field. We need to be able to
predict the swarm’s movements, and how those movements change as we change
parameters – to predict its collective dynamics and bifurcations. Predicting
these however is hard, because of swarms’ numerous degrees of freedom and
nonlinear particle interactions. It takes us into the world on nonequilibrium
statistical mechanics and many-body dynamical systems where standard tools and
techniques fail. Take magnetic domain walls. Each one obeys the integro-
differential Landau-Lifshitz-Gilbert equation, so $N$ obey a set of coupled
such equations whose solution is virtually impossible. (For $N=1,2$ walls
approximations such as the $(q,\phi)$ model Slonczewski (1972); Nasseri et al.
(2018) have been derived, but bifurcations and scaling beyond $N>2$ is
difficult Nasseri et al. (2018); Haltz et al. (2021)). And for magnetic
colloids, because of the coupling to the host fluid, there is an extra Navier-
Stokes type equation added to the mix. Vicsek-type models are sometimes used
as approximations here Liu et al. (2021); Zöttl and Stark (2023); Sándor et
al. (2017); Reichhardt and Reichhardt (2016), but are still largely
intractable; order parameters and bifurcation are often computed numerically
Levis et al. (2019); Liebchen and Levis (2017); Liu et al. (2021); Sándor et
al. (2017), leaving solvable models of driven matter scarce.

This Letter helps close this research gap by introducing a model of driven
matter which may be solved exactly. Our approach is to study a deliberately
simplified model which hopefully captures behavior with some universality, as
opposed to a detailed model specific to magnetic particles or colloids. The
model’s form is inspired from studies of coupled oscillators Strogatz et al.
(1989); O’Keeffe et al. (2022); Yoon et al. (2022); Strogatz et al. (1988);
Strogatz and Westervelt (1989) which allows us to leverage new tools from that
field to solve it.

Consider $N$ particles moving in a 1D periodic domain obeying

$\displaystyle\dot{x}_{i}$
$\displaystyle=E-b\sin(x_{i}-\alpha_{i})+\frac{J}{N}\sum_{j}\sin(x_{j}-x_{i})\cos(\theta_{j}-\theta_{i}),$
(1) $\displaystyle\dot{\theta}_{i}$
$\displaystyle=E-b\sin(\theta_{i}-\beta_{i})+\frac{K}{N}\sum_{j}\sin(\theta_{j}-\theta_{i})\cos(x_{j}-x_{i}).$
(2)

Here $x_{i},\theta_{i}$ are the $i$-th particle’s position and phase,
respectively. This phase could represent the orientation of a particle, the
(in-plane) orientation of an electric or magnetic dipole, or be associated
with an internal rhythm, like the chemical oscillation on the surface of an
autocatalytic particle Jin et al. (2023). Domain disorder is modeled by the
$b\sin(\cdot)$ terms, which pin $x_{i}$ and $\theta_{i}$ to sites
$\alpha_{i},\beta_{i}$ (we do not consider thermal disorder), and external
driving by the $E$ terms. The Kuramoto $\sin(\cdot)\cos(\cdot)$ terms capture
particle interactions. For the phases $\theta_{i}$, this creates
synchronization which depends on the distance, for the positions $x_{i}$,
aggregation which depends on the particles’ phases, natural choices of
interaction because they occur in diverse systems such as Janus particles Yan
et al. (2012) and Quincke rollers Zhang et al. (2020). We ignore excluded
volume interactions for simplicity.

Figure 1: States and bifurcations for $J=-K$. (a) Pinned ($E=0.4,K=0.5$), (b)
half-pinned ($E=0.4,K=1.1$), (c) sync ($E=0.4,K=2.1$), (d) periodic in $x$,
fixed in $\theta$ ($E=0.7,K=1.0$), (e) periodic in both $x$ and $\theta$
($E=0.7,K=1.8$), (f) chaotic ($E=1.1,K=1.0$). Lyapunov exponents highlighting
the bifurcations for two different values of $E$. (g) $E=0.7$, and (h)
$E=1.1$. The system is integrated with $(dt,T)=(0.01,100)$ using an RK4
method.

Small $N$ regime. We explore the little $N$ limit with a case study of $N=2$
particles. The physical system we have in mind here is two magnetic domain
walls moving on circular race track memory Zhang et al. (2018), similar to
recent experiments Hrabec et al. (2018), 111There, however, the spatial domain
was a straight line; here it is periodic. For simplicity, we study symmetric
pinning sites $(\alpha_{1},\alpha_{2})=(\beta_{1},\beta_{2})=(0,\pi)$ 222Which
may be appropriate for periodic substrates Reichhardt and Reichhardt (2016);
Hoffmann et al. (2000). Setting $b=1$ without loss of generality and moving to
coordinates $(x_{\pm},\theta_{\pm})=((x_{1}\pm
x_{2})/2,(\theta_{1}\pm\theta_{2}$)/2) yields

$\displaystyle\dot{x}_{+}=E-\cos x_{+}\sin x_{-},$ (3)
$\displaystyle\dot{x}_{-}=-\sin x_{+}\cos x_{-}-\frac{J}{2}\sin 2x_{-}\cos
2\theta_{-},$ (4)
$\displaystyle\dot{\theta}_{+}=E-\cos\theta_{+}\sin\theta_{-},$ (5)
$\displaystyle\dot{\theta}_{-}=-\sin\theta_{+}\cos x_{-}-\frac{K}{2}\sin
2\theta_{-}\cos 2x_{-}.$ (6)

The behavior of the model divides into two cases depending on the relative
sign of $J$ and $K$. We present the opposite sign case here (it contains the
most relevant physics) and the same-sign case in the supplementary material
(SM). We also set $J=-K$ for ease, since the magnitude of $J$ does not change
the overall phenomena (see SM).

First we develop some intuition for our system by visualizing its dynamics.
Imagine the particles as moving dots in the $(x,\theta)$ plane with periodic
boundary conditions 333Equivalently, the torus. Limit case dynamics are easy
to picture. When the driving dominates $E\gg b,K$, the particles will be swept
around the plane in uniform rotations. When the pinning is large $b\gg E,K$,
they will freeze into their pinning sites $x_{i}=\theta_{i}=\alpha_{i}$. And
when the coupling wins out $K\gg b,E$, they will unstick from $\alpha_{i}$ and
lock into synchrony 444since $J=-K$, it won’t be straightforward synchrony
$x_{1}=x_{2},\theta_{1}=\theta_{2}$, which we expect for $J=K$, but some type
of anti-sync. But what happens when the three effects have comparable strength
$b\approx K\approx E$? And how do the various states arise and disappear as
the parameters change?

To answer these questions, we ran numerical experiments. Figure 1 shows our
results, but before we talk through them, you should look at the bifurcation
diagram in Fig. 2 which encapsulates all the results at once. Having this
birds eye view of the systems dynamics in mind as you read will help you
follow the story. Supplementary Movie 1 also presents a live demo of our
experiments which is also useful to watch at this point.

Figure 2: Bifurcation diagram for $N=2$ particles. Black lines denote
theoretical predictions, colors signify the largest Lyapunov exponent. In the
periodic / chaotic parameter regions, bistability was sometimes observed.

To begin, we realized the pinned state – a natural ‘ground state’ to perturb
around – by turning off the driving $E=0$ and setting the phase coupling $K$
small. Fig. 1(a) shows time series of the coordinates $x_{\pm},\theta_{\pm}$
relaxing from random IC to the pinning sites, depicted as dotted lines (in the
$\pm$ coordinates, the sites are $\alpha_{\pm}=\pm\pi/2$). Then we gradually
increased $K$, expecting the particles to synchronize, in the sense of
minimizing their space differences $x_{-}$ and phase differences $\theta_{-}$.
(Note however that since driving is turned off, we don’t expect the particles
to oscillate here, we expect them to just shift to a new fixed point). We find
two transitions. First, a half pinned state emerges, where $\theta_{\pm}$ stay
pinned, but the phases $x_{\pm}$ are locked into a new sync fixed point (Fig.
1(b)). Notice the $\theta_{\pm}$ line settles onto the pinning sites (dotted
line), but $x_{\pm}$ does not. In the second transition, the remaining
coordinates depin and a full sync state is realized where all four coordinate
lie off the pinning sites (Fig. 1(c)).

Figure 3: Match to micromagnetic simulations of two domain walls. Top row
(a)-(f) high resolution micromagnetic simulations from Hrabec et al. (2018).
Middle row (g)-(l), $J=-K$ and bottom row (m)-(r), $J=K$ model Eqs. (1)-(2).
Numerical parameters: $(dt,T)=(0.25,200)$. Initial $(x_{i},\theta_{i})$ were
drawn uniformly at random from $[0,2\pi]$. Multistability was sometimes
observed. Reprinted Supplementary Figs. 7,8 with permission from [Aleš Hrabec
et al., Physical Review Letters, 120, 227204 (2018)]. Copyright (2023) by the
American Physical Society.

How much driving $E>0$ can the three states withstand before
$x_{\pm},\theta_{\pm}$ unlock and start moving? Look at the bifurcation
diagram in Fig. 2. For $E\lesssim 0.25$, the pinned $\rightarrow$ half-pinned
$\rightarrow$ sync transition persists, but for larger driving $E\approx 0.5$,
some periodic states (green region) arise between the half-pinned and sync
states. Here, either two of the four coordinates $(x_{\pm},\theta_{\pm})$
oscillate, the others remaining fixed (Fig. 1(d)), or all four coordinates
oscillate (Fig. 1(e)). For larger driving still $E\approx 0.8$, the pinned
state morphs directly into periodicity, which in turn undergoes an
intermittency transition to chaos (orangish region) depicted in Fig. 1(f).
Finally, for $E>1$, the static states vanish and the dynamic ones become
reentrant: the system flip flops between chaotic and periodic motion, then
settles into chaos (bistability between the two states was sometimes
observed). To confirm the motion was chaotic, we computed a heat map of
Lyapunov exponents $\lambda_{max}$ in the $(K,E)$ plane and saw
$\lambda_{max}>0$ where expected. We also compute power spectra which indicate
the chaos transition is of the intermittent type (Fig. S2). For convenience,
we also plot the Lyapunov exponents $\lambda(K)$ to show the bifurcation
sequence at $E=0.7$ and $E=1.1$ in Fig. 1(g)-(h).

Let’s take stock of our findings. We find three static states, a family of
periodic states, chaos, and various interstate transitions. To check these
states’ robustness, we run simulations for uneven coupling $J=-cK$ for $c\neq
1$ and asymmetric pinning $(\alpha_{1},\alpha_{2})=(0,a\pi)$ for $0\leq a\leq
1$ 555We can set $\alpha_{1}=0$ without loss of generality and find the same
physics (qualitatively identical bifurcation diagrams; see Figs. S5,S6,S7). A
surprise for asymmetric pinning (Fig. S6) is that the bifurcations of the
static states get richer; the half-pinned and sync states become reentrant.
Most of the collective states also appear for same-sign coupling $J=K$, albeit
with a different bifurcation structure (Fig. S1).

Our hope was that these states capture real world behaviour. They do. Figure 3
shows they mimic the behavior found in a recent study Hrabec et al. (2018) of
a pair magnetic domain walls. Briefly, their setup is this. Each wall is free
to move in the $x$ direction, and one wall is slightly higher than the other
in the $y$ direction, so that they do not collide. Then each wall $i$ may be
characterized by a single spatial degree of freedom $x_{i}\in\mathbb{R}^{1}$
and also a phase $\theta_{i}\in\mathbb{S}^{1}$ corresponding to the effective
magnetic dipole vector of the wall; thus the walls fall into our model class
666see Hrabec et al. (2018) for more details on the phase $\theta_{i}$. The
walls begin pinned at fixed points $(x_{i}^{*},\theta_{i}^{*})$. Then an
external magnetic field is turned on, which induces their positions $x_{i}$
and $\theta_{i}$ to unlock and interact. The top row of Fig. 3, taken from
Hrabec et al. (2018), shows the resulting dynamics for different parameters
with scatter plots of the dipole vectors for each wall $(m_{1},m_{2})$. These
$m_{i}$ relate to the phase via $m_{i}\propto\sin\theta_{i}$, as indicated by
the axes labels. Notice the walls settle into simple periodic behavior or more
complex dynamics represented by Lissajous curves and point clouds. The middle
column shows our model reproduces these states for the $J=-K$ case study we
have presented in the main text. The bottom row shows the same-sign $J=K$
case, whose analysis, recall, we present in the SM, can also reproduce the
states.

Beyond magnetic domain walls, our model can also reproduce realistic force-
velocity curves $\bar{v}(E)$ which are standard ways to analyze pinned-driven
systems (here $\bar{v}$ is the average system velocity and $E$ is the driving
force). Fig. S9 shows the $\bar{v}(E)$ of our model imitates that of pinned
superconductor vortices Reichhardt and Reichhardt (2016), for example.

Now we turn to analysis. We derive fixed point expressions for
$x_{\pm},\theta_{\pm}$ in the static states and derive their bifurcation
curves drawn as black lines in Fig. 2. We find the fixed points by setting the
RHS of Eqs. (3)-(5) to zero. Then we eliminate $(x_{-},\theta_{-})$ using Eqs.
(3), (5) and sub the result into Eqs. (4), (6):

$\displaystyle\sqrt{1-E^{2}\sec^{2}x_{+}}\Big{(}\frac{EK\left(2E^{2}\sec^{2}\theta_{+}-1\right)}{\cos
x_{+}}-\sin x_{+}\Big{)}=0,$ (7)
$\displaystyle\sqrt{1-E^{2}\sec^{2}\theta_{+}}\Big{(}\frac{EK\left(2E^{2}\sec^{2}x_{+}-1\right)}{\cos\theta_{+}}-\sin\theta_{+}\Big{)}=0.$
(8)

The systems has form $AB=CD=0$, implying four sets of fixed points. When
$(A,C)=(0,0)$ we get the pinned state. When the $(A,D)=(0,0)$ or $(B,C)=(0,0)$
we get the half pinned state (in the one $x_{i}$ stays pinned, $\theta_{i}$
syncs, in the other the reverse), and when $(B,D)=(0,0)$ we get the sync
state. The analysis of each state is complex, but the strategy is simple: find
the fixed points, and then their stability by linearization. So we present
just a sketch of the analysis for the pinned state here, and put the rest in
the SM.

The condition $(A,C)=(0,0)$ implies
$\sqrt{1-E^{2}\sec^{2}x_{+}}=\sqrt{1-E^{2}\sec^{2}\theta_{+}}=0$. Solving
these yields a family of 16 fixed points,
$(x_{+},x_{-},\theta_{+},\theta_{-})=(\pm\cos^{-1}(\pm
E),\pm\pi/2,\pm\cos^{-1}(\pm E),\pm\pi/2)$, four of which are stable,
reflecting a single macrostate. Linearization reveals a zero eigenvalue at
$E_{c}=\sqrt{1-K^{2}}$. Turning the same crank for the other states reveals
two-piece boundaries: for the half-pinned state, $E_{c}(K)=1/2K$ and
$E_{c}(K)=\sqrt{-4+5K^{2}-K^{4}}/3K$, intersecting at $K^{*}=\sqrt{5/2}$, and
for the sync state, $E_{c}=K/4$,
$E_{c}=1/2\sqrt{(K(\sqrt{K^{2}-4}+K)-1)/K^{2}}$ which meet at
$K^{*}=\sqrt{2(1+\sqrt{2})}$. The expressions for $x_{\pm,},\theta_{\pm}$ for
these states are too complex to display here. This completes our analysis of
the small $N$ regime.

Figure 4: Large $N$ limit for $J=K$. (a) Pinned state (blue, $K=1$) and anti-
pinned state (red, $K=-2$) plotted on same graph to save space for $E=0$ and
$N=200$ particles. (b) Sync state for $(K,E)=(3,0)$. (c) Order parameters for
$E=0$. (d) Bifurcation diagram in $(K,E)$ space. Black curves show theoretical
predictions. Colors denote Lyapunov exponents, computed using $N=10$
particles, which well approximated the $N\gg 1$ limit (simulations for larger
$N$ were prohibitive).

Large N regime. How do our results change as we scale up the population size
$N\gg 1$? A surprise is the $J=-K$ case we have thus far considered gets
simpler. The half-pinned and sync states disappear, leaving just the pinned
and unsteady states (Fig. S8). We used linearly spaced pinning
$\alpha_{i}=2\pi i/N$ here to facilitate analysis (randomly chosen
$\alpha_{i},\beta_{i}$ produce similar results Sar et al. (2023)).

Dynamics for same sign coupling $J=K$, in contrast, get richer. Fig. 4(a)
shows the pinned state (blue dots) persists, but now a new anti-pinned state
(red dots) arises where neighboring particles are shifted by an amount
$x_{i+1}-x_{i}=\theta_{i+1}-\theta_{i}=\Delta(K)$. Fig. 4(b) shows the large
$N$ analogue of the sync state with particles bunching into two sync clusters.
New periodic behavior is observed, along with quasiperiodicity, and the route
to chaos is now via period doubling. Supplementary Movies 2,3 depict the
evolution of all the states and Fig. S4 studies the transition to chaos.

We analyzed the $N\rightarrow\infty$ states in a previous work Sar et al.
(2023). The critical coupling for the anti-pinned state is
$E_{c}=\sqrt{1-K^{2}}$, found using a self-consistency analysis, and for the
pinned state is $E_{c}=\sqrt{1-K^{2}/4}$, found using a variational argument
777Note, in Sar et al. (2023) we called pinned, anti-pinned states the phase
wave, and split phase wave. These are tools from coupled oscillator theory we
advertised in the introduction.. These are the sides of the lop-sided bell in
the bifurcation diagram Fig. 4(d). We also derived expressions for the order
parameters

$\displaystyle
S_{\pm}e^{i\Phi_{\pm}}=\frac{1}{N}\sum_{j}e^{i(x_{j}\pm\theta_{j})}.$ (9)

In the split phase wave,

$\displaystyle S_{-}$
$\displaystyle=-\frac{1}{3}+\frac{K^{4}+(\Gamma_{1}+3\sqrt{3}\sqrt{\Gamma_{2}})^{2/3}}{3K^{2}(\Gamma_{1}+3\sqrt{3}\sqrt{\Gamma_{2}})^{1/3}},$
(10)

where

$\displaystyle\Gamma_{1}$ $\displaystyle:=-27\left(E^{2}-1\right)K^{4}-K^{6},$
(11) $\displaystyle\Gamma_{2}$
$\displaystyle:=\left(E^{2}-1\right)K^{8}\left(27E^{2}+2K^{2}-27\right).$ (12)

In the sync state, we derived a pair of self-consistency equations for $S_{+}$

$\displaystyle-$ $\displaystyle
2\sin\Big{(}\frac{\xi}{2}-\alpha\Big{)}+KS_{+}\sin(\Phi_{+}-\xi)=0,$ (13)
$\displaystyle
S_{+}e^{\tilde{i}\Phi_{+}}=\frac{1}{2\pi}\int_{0}^{2\pi}e^{\tilde{i}\xi(\alpha)}d\alpha,$
(14)

where $\xi_{i}:=x_{i}-\theta_{i}$ which we solved numerically 888We must solve
Eq. (13) for the fixed points $\xi^{*}(\alpha)$ then plug them into Eq. (14)
to find $S_{+}$. First we set $\Phi_{+}=0$ without loss of generality. Then by
applying various trig identities to Eq. (13) we arrive 4-th order polynomial
in $\cos\xi^{*}$ and plug the roots into Eq. (14), which when $\Phi_{+}=0$
reads $S_{+}=\int\cos(\xi^{*}(\alpha))d\alpha$. Then we computed the integral
over $\alpha$ numerically.. Figure 4(c) shows these $S_{\pm}(K)$ match
simulations and distinguish between the static states; in the pinned state
$(S_{+},S_{-})=(1,0)$ trivially by subbing $x_{i}=\theta_{i}$ into the
definition for $S_{\pm}$ Eq. (9) . There is also a small region of hysteresis
(not visible in the graph) between the pinned and sync states Sar et al.
(2023). This completes our analysis.

The reaction of particulate matter to external driving is crucial for
applications, yet difficult to understand theoretically. This Letter sheds
light on this class of dynamics with a toy model tractable in both the low and
large $N$ limit – moreover, the intermediary $N$ regime is surprisingly well
approximated by the $N\rightarrow\infty$ model; Figures S8, S10 show as few as
$N\approx 4$ particles give the same physics. The model also captures the
behavior of real world systems such as magnetic domain walls (Fig. 3),
Japanese tree frogs Aihara (2009) and Janus matchsticks Chaudhary et al.
(2014) (both realize the anti-pinned state Fig. 4(a)), and superconductor
vortices ($\bar{v}(E)$ curve in Fig. S9). We also suspect its chaos may be
connected to the active turbulence of biological microswimmers Das et al.
(2020), since the swimmers contain the same basics physics as the model:
driving and emergence in environments with pinning.

Given this balance between solvability and realism, we wonder if the model
‘could be the Kuramoto model’ for this class of driven matter, by which we
mean the simplest, representative model in a universality class. Future work
could explore this conjecture by deriving our model from a physically rigorous
model. Start with say the Landau-Lifshitz-Gilbert equations for magnetic
particles, exploit a small quantity like a weak coupling limit or a
separations of time scales using a perturbative method Kuramoto (1984), and
see if our model or something close to it pops out. This was the way the
Kuramoto model itself was derived – starting with a general reaction diffusion
equation and simplifying using phase reduction methods – and is the source so
to speak of its universality Kuramoto (1984).

Our model could guide experimental work on systems of magnetic domains walls
and other particles Hrabec et al. (2018); Ledesma-Motolinía et al. (2023).
Chaos has a niche application in such systems, it can be exploited for
hardware security Ghosh (2016)), but to our knowledge has not yet been
reported in multiparticle experiments 999Chaos has been observed in single
particle systems Chen et al. (2022); Shen et al. (2020). Our model predicts it
occurs for all $N>1$ and gives parameter regimes where it is likely to arise.
The model or a close variant may also be useful in studies of charge density
waves which couple to lattice vibrations Hedayat et al. (2019). Eq. (2) for
$\dot{\theta}_{i}$ has already been used to model the phase of the charge wave
Strogatz et al. (1988); Strogatz and Westervelt (1989); Strogatz et al.
(1989); the novelty would be Eq. (1) for $\dot{x_{i}}$ which could represent
the displacements of the lattice atoms from their equilibrium positions. Such
displacements form sinusoidal patterns Hedayat et al. (2019), encapsulated by
the $b\sin(x_{i}-\beta_{i})$ term, which competes with the tendency to
‘synchronize’ at equilibrium $x_{i}=x_{j}=0$, as per the $K\sin(x_{j}-x_{i})$
term.

Reproducibility. Code used for simulations and analytic calculations available
at github 101010https://github.com/Khev/swarmalators/tree/master/1D/on-
ring/random-pinning.

## References

  * Allwood et al. (2005) D. A. Allwood, G. Xiong, C. Faulkner, D. Atkinson, D. Petit, and R. Cowburn, Science 309, 1688 (2005).
  * Luo et al. (2020) Z. Luo, A. Hrabec, T. P. Dao, G. Sala, S. Finizio, J. Feng, S. Mayr, J. Raabe, P. Gambardella, and L. J. Heyderman, Nature 579, 214 (2020).
  * Wiesendanger (2016) R. Wiesendanger, Nature Reviews Materials 1, 1 (2016).
  * Reichhardt et al. (2022) C. Reichhardt, C. J. O. Reichhardt, and M. Milošević, Reviews of Modern Physics 94, 035005 (2022).
  * Reddy et al. (2012) L. H. Reddy, J. L. Arias, J. Nicolas, and P. Couvreur, Chemical Reviews 112, 5818 (2012).
  * Yang et al. (2017) T.-H. Yang, S. Zhou, K. D. Gilroy, L. Figueroa-Cosme, Y.-H. Lee, J.-M. Wu, and Y. Xia, Proceedings of the National Academy of Sciences 114, 13619 (2017).
  * Bricard et al. (2015) A. Bricard, J.-B. Caussin, D. Das, C. Savoie, V. Chikkadi, K. Shitara, O. Chepizhko, F. Peruani, D. Saintillan, and D. Bartolo, Nature Communications 6, 7470 (2015).
  * Caleap and Drinkwater (2014) M. Caleap and B. W. Drinkwater, Proceedings of the National Academy of Sciences 111, 6226 (2014).
  * Li et al. (2015) J. Li, O. E. Shklyaev, T. Li, W. Liu, H. Shum, I. Rozen, A. C. Balazs, and J. Wang, Nano Letters 15, 7077 (2015).
  * Pinto et al. (2020) M. Pinto, P. Ramalho, N. Moreira, A. Gonçalves, O. Nunes, M. Pereira, and O. Soares, Environmental Advances 2, 100010 (2020).
  * Cheng et al. (2014) R. Cheng, W. Huang, L. Huang, B. Yang, L. Mao, K. Jin, Q. ZhuGe, and Y. Zhao, ACS nano 8, 7746 (2014).
  * Manamanchaiyaporn et al. (2021) L. Manamanchaiyaporn, X. Tang, X. Yan, and Y. Zheng, IEEE Robotics and Automation Letters (2021).
  * Slonczewski (1972) J. Slonczewski, in _AIP Conference Proceedings_ (American Institute of Physics, 1972), vol. 5, pp. 170–174.
  * Nasseri et al. (2018) S. A. Nasseri, E. Martinez, and G. Durin, Journal of Magnetism and Magnetic Materials 468, 25 (2018).
  * Haltz et al. (2021) E. Haltz, S. Krishnia, L. Berges, A. Mougin, and J. Sampaio, Physical Review B 103, 014444 (2021).
  * Liu et al. (2021) Z. T. Liu, Y. Shi, Y. Zhao, H. Chaté, X.-q. Shi, and T. H. Zhang, Proceedings of the National Academy of Sciences 118, e2104724118 (2021).
  * Zöttl and Stark (2023) A. Zöttl and H. Stark, Annual Review of Condensed Matter Physics 14, 109 (2023).
  * Sándor et al. (2017) C. Sándor, A. Libal, C. Reichhardt, and C. O. Reichhardt, Physical Review E 95, 032606 (2017).
  * Reichhardt and Reichhardt (2016) C. Reichhardt and C. O. Reichhardt, Reports on Progress in Physics 80, 026501 (2016).
  * Levis et al. (2019) D. Levis, I. Pagonabarraga, and B. Liebchen, Physical Review Research 1, 023026 (2019).
  * Liebchen and Levis (2017) B. Liebchen and D. Levis, Physical Review Letters 119, 058002 (2017).
  * Strogatz et al. (1989) S. H. Strogatz, C. M. Marcus, R. M. Westervelt, and R. E. Mirollo, Physica D: Nonlinear Phenomena 36, 23 (1989).
  * O’Keeffe et al. (2022) K. O’Keeffe, S. Ceron, and K. Petersen, Phys. Rev. E 105, 014211 (2022).
  * Yoon et al. (2022) S. Yoon, K. P. O’Keeffe, J. F. F. Mendes, and A. V. Goltsev, Phys. Rev. Lett. 129, 208002 (2022).
  * Strogatz et al. (1988) S. Strogatz, C. Marcus, R. Westervelt, and R. Mirollo, Physical Review Letters 61, 2380 (1988).
  * Strogatz and Westervelt (1989) S. H. Strogatz and R. M. Westervelt, Physical Review B 40, 10501 (1989).
  * Jin et al. (2023) H. Jin, J. Cui, and W. Zhan, Langmuir (2023).
  * Yan et al. (2012) J. Yan, M. Bloom, S. C. Bae, E. Luijten, and S. Granick, Nature 491, 578 (2012).
  * Zhang et al. (2020) B. Zhang, A. Sokolov, and A. Snezhko, Nature communications 11, 4401 (2020).
  * Zhang et al. (2018) S. Zhang, W. Wang, D. Burn, H. Peng, H. Berger, A. Bauer, C. Pfleiderer, G. Van Der Laan, and T. Hesjedal, Nature Communications 9, 2115 (2018).
  * Hrabec et al. (2018) A. Hrabec, V. Křižáková, S. Pizzini, J. Sampaio, A. Thiaville, S. Rohart, and J. Vogel, Physical Review Letters 120, 227204 (2018).
  * Sar et al. (2023) G. K. Sar, D. Ghosh, and K. O’Keeffe, Physical Review E 107, 024215 (2023).
  * Aihara (2009) I. Aihara, Physical Review E 80, 011918 (2009).
  * Chaudhary et al. (2014) K. Chaudhary, J. J. Juárez, Q. Chen, S. Granick, and J. A. Lewis, Soft Matter 10, 1320 (2014).
  * Das et al. (2020) A. Das, J. Bhattacharjee, and T. Kirkpatrick, Physical Review E 101, 023103 (2020).
  * Kuramoto (1984) Y. Kuramoto, _Chemical Oscillations, Waves, and Turbulence_ (Springer, 1984).
  * Ledesma-Motolinía et al. (2023) M. Ledesma-Motolinía, J. Carrillo-Estrada, A. Escobar, F. Donado, and P. Castro-Villarreal, Physical Review E 107, 024902 (2023).
  * Ghosh (2016) S. Ghosh, Proceedings of the IEEE 104, 1864 (2016).
  * Hedayat et al. (2019) H. Hedayat, C. J. Sayers, D. Bugini, C. Dallera, D. Wolverson, T. Batten, S. Karbassi, S. Friedemann, G. Cerullo, J. van Wezel, et al., Physical Review Research 1, 023029 (2019).
  * Hoffmann et al. (2000) A. Hoffmann, P. Prieto, and I. K. Schuller, Physical Review B 61, 6958 (2000).
  * Chen et al. (2022) Y. Chen, X. Ge, W. Luo, S. Liang, X. Yang, L. You, and Y. Zhang, arXiv preprint arXiv:2207.02492 (2022).
  * Shen et al. (2020) L. Shen, J. Xia, X. Zhang, M. Ezawa, O. A. Tretiakov, X. Liu, G. Zhao, and Y. Zhou, Physical Review Letters 124, 037202 (2020).