# Fast Muon Tracking with Machine Learning Implemented in FPGA

Chang Sun  Takumi Nakajima  Yuki Mitsumori  Yasuyuki Horii
<EMAIL_ADDRESS>Makoto Tomoto  Nagoya University, Chikusa-ku,
Nagoya 464-8602, Japan  High Energy Accelerator Research Organisation (KEK),
Oho, Tsukuba 305-0801, Japan

###### Abstract

In this work, we present a new approach for fast tracking on multiwire
proportional chambers with neural networks. The tracking networks are
developed and adapted for the first-level trigger at hadron collider
experiments. We use Monte Carlo samples generated by Geant4 with a custom muon
chamber, which resembles part of the thin gap chambers from the ATLAS
experiment, for training and performance evaluations. The chamber has a total
of seven gas gaps, where the first and last gas gaps are displaced by $\sim
1.5$ m. Each gas gap has 50 channels with a size of 18–20 mm. Two neural
network models are developed and presented: a convolutional neural network and
a neural network optimized for the detector configuration of this study. In
the latter network, a convolution layer is provided for each of three groups
formed from 2–3 gas gaps of the chamber, and the outputs are fed into
multilayer perceptrons in sequence. Both networks are transformed into
hardware description language and implemented in Virtex UltraScale+ FPGA. The
angular resolution is 2 mrad, which is comparable to the maximum resolution of
the detector estimated by the minimum $\chi^{2}$ method. The latency achieved
by the implemented firmware is less than 100 ns, and the throughput rate is
160 MHz.

###### keywords:

Muon Tracking, Machine Learning, Field-Programmable Gate Array, Artificial
Neural Network

##  1 Introduction

Triggering muons is critically important in proton-proton collider
experiments. For instance, in 2012, a new particle with a mass of
approximately 125 GeV was observed at the ATLAS and CMS experiments [1, 2]
using a muon trigger [3, 4] at the CERN Large Hadron Collider (LHC).
Subsequent measurements indicate that the new particle is consistent with the
Higgs boson in the standard model [5, 6].

An upgrade to the LHC, which is known as the High-Luminosity LHC (HL-LHC), is
expected to become operational in 2027 [7]. The peak luminosity of the HL-LHC
would ultimately be increased to $7.5\times 10^{34}~{}{\rm cm}^{-2}{\rm
s}^{-1}$, which is 7.5 times the design luminosity of the LHC ($1\times
10^{34}~{}{\rm cm}^{-2}{\rm s}^{-1}$). The trigger and data-acquisition system
of the ATLAS experiment requires upgrades [8] to accommodate the increase in
luminosity. In the original first-level muon trigger, only selected and
processed hit data of the muon spectrometer are transferred from the frontend
to the backend electronics. However, the new system will send all hit data to
the backend electronics, exploiting advances in data-transfer technology, and
all hits will be processed using high-end FPGAs integrated into the backend
electronics.

The ability of backend FPGAs to retrieve full hit data provides a unique
opportunity to boost the first-level muon trigger’s performance. A throughput
of $\geq 40~{}{\rm MHz}$ and a latency within ${\cal O}(100~{}{\rm ns})$ are
required. One method proposed for the new system is a fast tracking algorithm
known as “pattern matching” [8], which compares the received hit data with
predefined hit lists with corresponding track information assigned. Another
possibility, which is the subject of this paper, is to use a neural network to
perform fast tracking, exploiting the augmented on-board computing power.
Here, the hit data are used as the neural network’s input, and the track
information is outputted.

Machine learning techniques are rapidly being developed for and applied in
high-energy physics experiments. A living review is provided in Ref. [9].
Refs. [10, 11, 12, 13] show how different architectures can be deployed on a
first-level trigger using the hls4ml library. These works have opened a new
era for using machine learning for trigger systems. Novel fast tracking on a
silicon detector has been and continues to be developed [8, 14, 15, 16];
however, the use of machine learning is challenging because of the large
number of channels. A muon chamber has fewer channels, and fast muon tracking
with machine learning has been attempted [17].

In this study, a novel approach of using neural networks for online muon
tracking is introduced. Neural network models are developed for the first-
level muon trigger at high-luminosity hadron colliders. Two neural network
models are implemented in Xilinx Virtex UltraScale+ FPGA and are tested on
Monte Carlo (MC) samples through Vivado post-implementation simulations [18].
The tracking performance, latency, and required FPGA resources are evaluated
with different hyper parameters. For the tracking performance, focus is placed
on the angular resolution, which is critical for the momentum determination
(and hence the trigger rate) at the first-level muon trigger [8]. The
dependence of the tracking performance on the detector noise level is also
studied.

The remainder of this paper is organized as follows. The model for the MC
sample as well as the tracking performance by a conventional minimum
$\chi^{2}$ method is described in Section 2. The workflow and software are
introduced in Section 3. The network design and evaluations for the two types
of networks are described in Sections 4 and 5. The conclusion is provided in
Section 6.

##  2 Simulation

###  2.1 Simulation samples

The networks presented in this work were trained and evaluated on events
simulated with Geant4 [19]. In the simulation, a gas chamber detector shape
similar to the thin gap chamber (TGC) of the ATLAS experiment [20] was formed
with an approximate pseudorapidity ($\eta$) coverage of $2.03\leq|\eta|\leq
2.26$. A full schematic of the detector model used is shown in Fig. 1.
Rectangular plates were used instead of the trapezoidal plates of the TGC
detector for simplicity. The angular resolution of the tracking shown in this
paper is independent from the simplifications of the detector geometry because
the wire orientations are identical. Three plates—one triplet and two doublet
plates—were used in the simulation. The triplet plate corresponds to the TGC’s
innermost plate and has dimensions of $900\text{ mm}\times 500\text{ mm}\times
70\text{ mm}$. The two doublet plates correspond to the TGC’s outer plates and
have dimensions of $990\text{ mm}\times 500\text{ mm}\times 44\text{ mm}$. The
triplet plate is named M1, and the doublet plates are named M2 and M3. All the
plates were positioned with their central points set to the TGC plates’
central positions. The $z$ positions111 The $z$ axis is taken as the proton
beam axis in the ATLAS experiment, and is taken accordingly in this study. of
the M1, M2, and M3 plates were identical to TGC’s plates (Fig. 2). The wire
pitch for all the plates was $1.8$ mm, and 10 (11) wires were grouped into
channels for triplet (doublet) plates. No magnetic field was introduced
because the TGC detector is located outside the toroidal magnet at the ATLAS
experiment.

Figure 1: Schematic of the detector model. The $z$ coordinates of the chambers
are identical to those of TGC detectors in the ATLAS experiment. The three
plates—M1, M2, and M3 in ascending $z$-position order—are aligned with their
normals parallel to the $z$ axis. The arrow represents an incoming muon, and
$\theta$ is defined as the polar angle. Figure 2: Cross-sectional view of the
triplet (left) and doublet (right) plates.

The muon source was placed in front of the M1 plate at $z=13.3$ m. Multiple
scattering and emissions before arriving at the muon chamber were ignored
because the intrinsic performance of the tracking with the muon chamber was
the topic of interest. At the ATLAS experiment, muons with transverse momentum
$p_{\rm T}\geq 20$ GeV are deflected, at most, $\pm 30$ mrad from the
trajectory of muons with infinite momentum when impacting the TGC detector.
The sign of the deflection angle depends on the electric charge of the muon.
In this study, muons’ direction was uniformly distributed in the range $\pm
30$ mrad from the angle of the straight track, which was independently
uniformly distributed on the polar angle $\theta$. The straight tracks were
aligned directly from the interaction point. A single muon was generated for
each event. Tracking of the minimum ionizing particle is critical for muon
triggerring at the ATLAS experiment, and all muons were generated with $p_{\rm
T}=20$ GeV in this study. All data associated with tracks passing through all
seven gas gaps were collected for both training and performance evaluation. A
channel generated a readout value of 1 if it was fired and 0 if not. The fired
channels are referred to as “hits” in this work.

With the aforementioned settings, $3\times 10^{6}$ events were generated for
this study. Twenty percent of the data was used as the test set, which was not
exposed to any models introduced in this study prior to the performance
analysis.

To achieve a more realistic background environment, we prepared extra
simulation samples with artificial noises. These samples are generated by
randomly setting the readout value to $1$ with probability $p$ for all
channels independently of the generated dataset. Hereafter, $p$ is represented
as the “noise level”.

The channel-wise correlation matrices of the channel outputs between the
plates are shown in Fig. 3. The elements of the correlation matrices, the
correlation coefficients $\rho_{i,j}$, are defined as $\rho_{i,j}\equiv
E[\\{x_{i}-E(x_{i})\\}\\{x_{j}-E(x_{j})\\}]/$
$\sqrt{E[\\{x_{i}-E(x_{i})\\}^{2}]E[\\{x_{j}-E(x_{j})\\}^{2}]}$, where $x$ is
$1$ if the corresponding channel has a hit and $0$ otherwise, $i$ and $j$ are
the channel identifiers for the given two plates, and $E$ represents the
expectation value. These correlation matrices were used for determining some
of the hyper parameters and for reducing the number of multiplication
operations required by the network.

Figure 3: Channel-wise correlation matrices for combinations of the M1, M2,
and M3 plates.

###  2.2 Reference tracking performance

For reference, the best possible resolution of this detector was estimated
using the minimum $\chi^{2}$ method. In each event, one straight line was
randomly placed in the detector’s volume according to the distribution of muon
trajectories, as described in Section 2.1, and all channels passed by the line
were marked as fired. Because there was no false hit, an ordinary $\chi^{2}$
fit was used to estimate the resolution. Specifically, when a gas gap had
exactly one hit, the coordinate of the center of the channel was used for the
fit with an error $\sigma$; if a gas gap had two adjacent hits, the coordinate
of the center of the two channels combined with an error $\sigma/\sqrt{2}$ was
used for the fit. This method is referred to as “ideal $\chi^{2}$” in this
work. The detector’s angular resolution determined by this method is 1.7 mrad.

The performance determined from a $\chi^{2}$-fitting method iterating through
all viable hit combinations was evaluated for comparison to the network in
this study on the simulated events. The exact algorithm used for this fit is
described as follows: 1) One hit was taken from each gas gap with non-zero
hits, and a simple $\chi^{2}$ fit was performed for all combinations. 2) The
fitted track with the lowest $\chi^{2}$ score was selected as the output.
Here, given the uniform channel sizes, the hits from the same plate were
regarded as having the same uncertainty. This method is referred to as
“iterative $\chi^{2}$” in this work. The angular resolution obtained in this
manner for the simulated events is 2.2 mrad. Although the iterative $\chi^{2}$
method is more resistant to noises, it uses a maximum of one hit per gas gap
and does not fully exploit the position information for a muon that provides
hits in two neighboring channels.

##  3 Workflow

All neural networks used in this study were first defined and trained using
Keras 2.4.3 [21] and QKeras 0.8.0 [13] with TensorFlow 2.2.0 [22] as the
backend. The networks were then translated into an intermediate C++ script
using hls4ml 0.5.0 [10] with the setup of ReuseFactor = 1, Strategy = Latency,
and IOType = io_parallel. The synthesis was performed with Vivado HLS 2020.1
[23]. The implementation was carried out with Vivado 2020.1 [18], targeting a
Xilinx Virtex UltraScale+ XCVU13P FPGA222XCVU13P-L2FHGA2104E with a clock
frequency of 160 MHz. The performance was evaluated with a post-implementation
functional simulation running on the Vivado software.

##  4 Muon tracking with convolutional neural networks

###  4.1 Network design

A convolutional neural network (CNN) is a class of deep neural network, widely
applied to the analyses of visual images. In this study, the detector outputs
were considered to be visual images and muon tracks were reconstructed by a
CNN. The input of the network was the hit information for all channels
$\langle 50\times 7\rangle$ and the output was the muon track angle $\theta$.

A large-scale network that maximizes the muon tracking performance was
designed; its schematic is shown in Fig. 4(a). This design was primarily
inspired by the VGG network [24]. Multiple convolution layers extract the
features of the muon track, and the affine (also known as fully connected)
layers located after the convolution layers extract the muon track angle from
the extracted features. This network is denoted as software (SW).

(a) Large-scale CNN model (SW) (b) Compact CNN model (BL, QF8, QF6)

Figure 4: Schematics of the large-scale and compact CNN models. Both are
composed of two convolution (“Conv1D”) layers and four affine (“Affine”)
layers. The kernel size for the Affine layers refers to $\langle$output
dimension, input dimension$\rangle$. The kernel size for Conv1D layers refers
to filter size and $\langle$input channel, filter count$\rangle$. The <*>
notation on arrows indicates that the array shape passed through the arrow.

A compact version of the network was designed by downscaling the SW network
with the objective of reducing FPGA resource utilization and achieving a
smaller latency value. This network is referred to as baseline (BL), and its
schematic is shown in Fig. 4(b).

The filter sizes of convolution layers were chosen to cover almost all hits
from a single muon track. The exact numbers were determined by the channel-
wise correlation matrix shown in Fig. 3.

Quantization was performed on the compact network to further reduce FPGA
resource utilization and latency. All parameters, including weights, biases,
and intermediate values, in the quantized networks were quantized to
ap_fixed<$\cdot,\cdot$> or ap_ufixed<$\cdot,\cdot$>. Overflow and round
behaviors were set to be default with wraparound (AP_WRAP) and truncation to
minus infinity (AP_TRN). The numbers of integer bits assigned to each layer’s
weights, biases, and intermediate values were selected to be the minimum
number such that no overflow would occur for the entire training set. By
contrast, the number of float bits assigned ($n$) was fixed for all parameters
except accumulators, which had $(n+2)$ float bits instead of $n$. As the only
exception to this rule, the output layer always had $9$ unsigned integer bits
and $7$ float bits. A network quantized with $n$ float bits assigned is
denoted as QF[n]. Here, $n=8$ and $6$ were used for the compact network, and
the resultant quantized compact networks are referred to as QF8 and QF6,
respectively. Because no significant performance difference was observed
between post-training quantization and quantization-aware training, only post-
training quantization was performed.

All the networks were fully pipelined and able to accept an event every
$1/(160~{}{\rm MHz})=6.25~{}{\rm ns}$. No initialization was required between
the processing of consecutive events, which is a key feature of the designed
network to be used at the first-level trigger at hadron-collider experiments.

###  4.2 Performance

The angular resolution was evaluated from the widths of the distributions of
the difference between the CNN output and the MC truth angle. Figure 5 shows
the distributions for the QF8 network. The distributions of the ideal and
iterative $\chi^{2}$ methods are also shown for reference. Table 1 shows the
values of the angular resolution for the networks examined in this study. The
SW network provides the best angular resolution of 1.7 mrad. The angular
resolution of the BL, QF8, and QF6 networks ranges from 2.0 to 3.4 mrad.

Figure 5: Distributions of the angle difference between the MC truth track
segment and the reconstructed track segment. The solid (black) histogram shows
the distribution for the track segments reconstructed by the QF8 network. The
dashed (blue) and dotted (red) histograms show the distributions for track
segments reconstructed by the ideal and iterative $\chi^{2}$ methods,
respectively. Table 1: The angular resolution, latency, and resource
utilization for CNN examined in this study. The numbers in parentheses show
the rates of the resource utilization per super logic region. The angular
resolutions for SW and BL are estimated using software, whereas those for QF8
and QF6 are estimated using a Vivado post-implementation simulation. The SW
and BL networks are provided without the intention to be run on the hardware;
thus, the values of the latency and resource utilization are omitted.

Model | Resolution [mrad] | Latency [ns] | DSP48 | LUT | FF | BRAM
---|---|---|---|---|---|---
SW | 1.7 | - | - | - | - | -
BL | 2.0 | - | - | - | - | -
QF8 | 2.2 | 81 | 2087 (68%) | 66,441 (15%) | 19,849 (2%) | 0 (0%)
QF6 | 3.4 | 81 | 607 (20%) | 61,749 (14%) | 11,702 (1%) | 0 (0%)

The latency and the resource utilization were obtained for the quantized
networks and are summarized in Table 1. Both the QF8 and QF6 networks were
implemented on Vivado without timing violation and achieved a latency less
than 100 ns.

##  5 Muon tracking with multistage neural networks

###  5.1 Network design

We designed another network to take maximum advantage of the strong
correlation of the hits between plates of the detector. A schematic of the
network is shown in Fig. 6. The hit information on the channels from plates
M1, M2, and M3 were fed into the network separately in the form of 1-bit
arrays of shapes $\langle 50\times 3\rangle$, $\langle 50\times 2\rangle$, and
$\langle 50\times 2\rangle$, respectively. The network consists of two
building blocks: a feature-extraction network that exploits the layered
structure of the detectors and a generic fully-connected net that outputs the
final values.

Figure 6: Schematic of the neural network optimized for the detector
configuration. The kernel size for the Affine layers refers to $\langle$output
dimension, input dimension$\rangle$. The kernel size for the Conv1D layers
refers to filter size and $\langle$input channel, filter count$\rangle$. The
<*> notation on arrows indicates that the array shape passed through the
arrow. Arrows without an array shape specified are passing arrays of shape 50.

The design of the feature extraction part was inspired by the track-following
scheme [25], where multiple hits are combined together by detector stations
and added sequentially in the inside-out direction to form the full track. The
convolution layers extract real vectors encoding hit-position information from
each plate, whereas the “masked linear” layers (i.e., highly sparse affine
layers) define search windows and project the track information from one plate
to the next as real vectors. The previous track information (i.e., the output
of the last masked linear layer) is merged with the current observation (i.e.,
the output of the convolution layer) via a simple addition followed by a
hyperbolic tangent activation function to allow some inefficiency for the
detector while still favoring the track with more hits.

Channel-wise correlation matrices between the plates (Fig. 3) were used to
reduce the number of multiplication operations required by the networks. We
defined the search windows encoded in the masked linear layers by requiring a
correlation coefficient greater than $0.01$. All weights outside the windows
were constrained to be zero. This practice was responsible for the highly
sparse nature of the masked linear layers, where the sparsity was $\sim 90\%$.
A graphical explanation of how the masked linear layers work is provided in
Fig. 7.

Figure 7: Demonstration of how the mask matrices work. The three black lines
on the plates represent the observation and state vectors’ position. The gray
shaded cones represent the search region for a single hit on M1 or M2 marked
red.

The network with full precision in the software simulation is referred to as
the BL network. Three quantized networks were built from the BL network
following the procedure specified in Section 4, except that, after
quantization was applied to the trained BL network, further quantization-aware
training was performed to mitigate the performance drop. These networks were
quantized with 7, 5, and 3 float bits and are referred to as QF7, QF5, and
QF3, respectively. All the networks were fully pipelined, like the networks
described in Section 4, and were able to accept an event every $6.25~{}{\rm
ns}$.

###  5.2 Performance

The angular resolution was evaluated from the widths of the distributions of
the difference between the network output and the MC truth angle. Figure 8
shows the distributions for the QF7 and QF5 networks. The distribution of the
ideal $\chi^{2}$ method is also shown for reference. Table 2 shows the values
of the angular resolution for the networks examined in this study. The angular
resolution of the QF7, QF5, and QF3 networks ranges from 2.0 to 2.8 mrad.

Figure 8: Distributions of the angle difference between the MC truth track
segment and the reconstructed track segment. The solid (black) and dotted
(red) histograms show the distributions for the track segments reconstructed
by the QF7 and QF5 networks, respectively. The dashed (blue) histogram shows
the distribution for track segments reconstructed by the ideal $\chi^{2}$
method. Table 2: The angular resolution, latency, and resource utilization for
the neural networks optimized for the detector configuration of this study.
The numbers in parentheses show the rates of resource utilization per super
logic region. The angular resolution for BL is estimated using software,
whereas those for QF7, QF5, and QF3 are estimated using a Vivado post-
implementation simulation. The BL network is provided without the intention to
be run on the hardware; thus, the values of the latency and resource
utilization are omitted.

Model | Resolution [mrad] | Latency [ns] | DSP48 | LUT | FF | BRAM
---|---|---|---|---|---|---
BL | $1.9$ | - | - | - | - | -
QF7 | $2.0$ | 69 | 1389 (45%) | 34,848 (8.1%) | 5433 (0.6%) | 75 (5.6%)
QF5 | $2.2$ | 69 | 88 (2.9%) | 40,039 (9.3%) | 3419 (0.4%) | 75 (5.6%)
QF3 | $2.8$ | 56 | 2 ($<0.1$%) | 21,682 (5.0%) | 2242 (0.3%) | 75 (5.6%)

For the TGC detector at the ATLAS experiment at the HL-LHC, an average angular
resolution of 4 mrad is assumed [8]. The angular resolution depends on the
detector region, and the value obtained for $2.13<|\eta|<2.16$, which covers
the $\eta$ range of this study, is 2.4 mrad [26]. Although the detector
configuration is not identical and thus the numbers cannot be compared
directly, an angular resolution slightly greater than 2 mrad achieved by the
QF7 and QF5 networks would be sufficient for the TGC detector.

The latency and the resource utilization were obtained for the quantized
networks; the results are summarized in Table 2. Each quantized network was
implemented on Vivado without timing violation and achieved a latency less
than 100 ns. The resources were allocated by Vivado HLS on the basis of the
number of bits involved in the calculations; some of multiplications performed
by DSP48 in the QF7 network are performed by LUTs in the QF5 and QF3 networks.

The QF7 network was trained and evaluated on the datasets with artificial
random noises. The model was retrained on the noise levels of $10^{-4}$ and
$10^{-3}$ and then implemented through Vivado. The angular resolution of the
model at each noise level is shown in Table 3. With a reasonable noise level
of $10^{-4}$, the network could still outperform the iterative $\chi^{2}$
method; it maintained comparable performance at a noise level of $10^{-3}$.

Table 3: Angular resolution in unit of milli-radians of the model trained on different datasets and implemented as QF7. The row headings are the noise level for the test sets, and the column headings are the noise level for the training sets. Training Test | 0 | $10^{-4}$ | $10^{-3}$
---|---|---|---
0 | 2.0 | 2.2 | 3.3
$10^{-4}$ | 2.0 | 2.1 | 2.6
$10^{-3}$ | 2.1 | 2.1 | 2.3

##  6 Conclusion

Neural networks for fast tracking with multiwire proportional chambers were
proposed and studied. A compact and quantized neural network was developed
with a conventional CNN structure, where all the detector hits were inputted
to a single convolution layer followed by multilayer perceptrons. A multistage
neural network was then provided, where the outputs of the convolution layers
for gas-gap groups were fed into multilayer perceptrons in sequence. The
multistage neural network demonstrated better resource utilization and shorter
latency for a given angular resolution of the tracking. The firmware was
implemented in Xilinx UltraScale+ FPGA without timing violation and provided
an angular resolution comparable with the best possible resolution of the
detector estimated by the $\chi^{2}$ method. A trigger turn-on curve similar
to or better than that in Fig. 8.14 in Ref. [8], which was obtained using the
$\chi^{2}$ method, is expected with the neural networks developed in this
study. The performance degraded moderately and was still acceptable with
random noise hits. The firmware was able to perform one reconstruction in
sub-$100$ ns with a fixed latency and have a high throughput of $160$ MHz.

This work demonstrates a technique of forming a neural network and reducing
the required resources for fast muon tracking using the hits of TGC-type
detectors. Although the present networks will need further modifications and
optimizations for a larger coverage of the detector and for a full trigger
chain, the techniques shown in this paper can be exploited to maintain muon
trigger thresholds in the HL-LHC era despite the greater luminosity.

## Acknowledgements

This work was supported by JSPS KAKENHI Grant Numbers 16H06493, 18K03675, and
21H05085.

## References

  * [1] ATLAS Collaboration, Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Physics Letters B 716 (1) (2012) 1–29.  doi:https://doi.org/10.1016/j.physletb.2012.08.020.
  * [2] CMS Collaboration, Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC, Physics Letters B 716 (1) (2012) 30–61.  doi:https://doi.org/10.1016/j.physletb.2012.08.021.
  * [3] ATLAS Collaboration, Performance of the ATLAS muon trigger in pp collisions at $\sqrt{s}=8$ TeV, The European Physical Journal C 75 (2015) 120.  doi:https://doi.org/10.1140/epjc/s10052-015-3325-9.
  * [4] CMS Collaboration, The CMS trigger system, Journal of Instrumentation 12 (01) (2017) P01020.  doi:https://doi.org/10.1088/1748-0221/12/01/p01020.
  * [5] ATLAS Collaboration, A detailed map of Higgs boson interactions by the ATLAS experiment ten years after the discovery, Nature 607 (2022) 52–59.  doi:https://doi.org/10.1038/s41586-022-04893-w.
  * [6] CMS Collaboration, A portrait of the Higgs boson by the CMS experiment ten years after the discovery, Nature 607 (2022) 60–68.  doi:https://doi.org/10.1038/s41586-022-04892-x.
  * [7] I. B. Alonso, O. Brüning, P. Fessia, M. Lamont, L. Rossi, L. Tavian, M. Zerlauth, High-Luminosity Large Hadron Collider (HL-LHC): Technical design report (2020).  doi:https://doi.org/10.23731/CYRM-2020-0010.
  * [8] ATLAS Collaboration, Technical Design Report for the Phase-II Upgrade of the ATLAS TDAQ System (2018).
URL https://cds.cern.ch/record/2285584

  * [9] A Living Review of Machine Learning for Particle Physics.
URL https://iml-wg.github.io/HEPML-LivingReview/

  * [10] J. Duarte, S. Han, P. Harris, S. Jindariani, E. Kreinar, B. Kreis, J. Ngadiuba, M. Pierini, R. Rivera, N. Tran, Z. Wu, Fast inference of deep neural networks in FPGAs for particle physics, Journal of Instrumentation 13 (07) (2018) P07027.  doi:https://doi.org/10.1088/1748-0221/13/07/p07027.
  * [11] S. Summers, G. D. Guglielmo, J. Duarte, P. Harris, D. Hoang, S. Jindariani, E. Kreinar, V. Loncar, J. Ngadiuba, M. Pierini, D. Rankin, N. Tran, Z. Wu, Fast inference of Boosted Decision Trees in FPGAs for particle physics, Journal of Instrumentation 15 (05) (2020) P05026.  doi:https://doi.org/10.1088/1748-0221/15/05/p05026.
  * [12] J. Ngadiuba, V. Loncar, M. Pierini, S. Summers, G. D. Guglielmo, J. Duarte, P. Harris, D. Rankin, S. Jindariani, M. Liu, K. Pedro, N. Tran, E. Kreinar, S. Sagear, Z. Wu, D. Hoang, Compressing deep neural networks on FPGAs to binary and ternary precision with hls4ml, Machine Learning: Science and Technology 2 (1) (2020) 015001.  doi:https://doi.org/10.1088/2632-2153/aba042.
  * [13] C. N. Coelho Jr, A. Kuusela, S. Li, H. Zhuang, J. Ngadiuba, T. K. Aarrestad, V. Loncar, M. Pierini, A. A. Pol, S. Summers, Automatic heterogeneous quantization of deep neural networks for low-latency inference on the edge for particle detectors, Nature Machine Intelligence 3 (2021) 675–686.  doi:https://doi.org/10.1038/s42256-021-00356-5.
  * [14] CMS Collaboration, The Phase-2 Upgrade of the CMS Level-1 Trigger (2020).
URL https://cds.cern.ch/record/2714892

  * [15] ATLAS Collaboration, The ATLAS Fast TracKer system, Journal of Instrumentation 16 (07) (2021) P07006.  doi:https://doi.org/10.1088/1748-0221/16/07/P07006.
  * [16] W. Ashmanskas, A. Bardi, M. Bari, S. Belforte, J. Berryhill, M. Bogdan, A. Cerri, A. Clark, G. Chlanchidze, R. Condorelli, R. Culbertson, M. Dell’Orso, S. Donati, H. Frisch, S. Galeotti, P. Giannetti, V. Glagolev, A. Leger, E. Meschi, F. Morsani, T. Nakaya, G. Punzi, L. Ristori, H. Sanders, A. Semenov, G. Signorelli, M. Shochet, T. Speer, F. Spinella, P. Wilson, X. Wu, A. Zanetti, Silicon vertex tracker: a fast precise tracking trigger for CDF, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 447 (1) (2000) 218–222.  doi:https://doi.org/10.1016/S0168-9002(00)00190-X.
URL https://www.sciencedirect.com/science/article/pii/S016890020000190X

  * [17] M. Migliorini, J. Pazzini, A. Triossi, M. Zanetti, A. Zucchetta, Muon trigger with fast Neural Networks on FPGA, a demonstrator (2021).  arXiv:2105.04428[hep-ex].
URL https://arxiv.org/abs/2105.04428

  * [18] Xilinx, Vivado Design Suite.
URL https://www.xilinx.com

  * [19] S. Agostinelli, et al., Geant4–a simulation toolkit, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 506 (3) (2003) 250–303.  doi:https://doi.org/10.1016/S0168-9002(03)01368-8.
  * [20] ATLAS Collaboration, The ATLAS Experiment at the CERN Large Hadron Collider, Journal of Instrumentation 3 (08) (2008) S08003.  doi:https://doi.org/10.1088/1748-0221/3/08/s08003.
  * [21] F. Chollet, et al., Keras (2015).
URL https://github.com/fchollet/keras

  * [22] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, software available from tensorflow.org (2015).
URL https://www.tensorflow.org/

  * [23] Xilinx, Vivado Design Suite User Guide: High-Level Synthesis (2021).
URL
https://www.xilinx.com/support/documentation/sw_manuals/xilinx2020_1/ug902-vivado-
high-level-synthesis.pdf

  * [24] K. Simonyan, A. Zisserman, Very deep convolutional networks for large-scale image recognition (2014).  arXiv:1409.1556[cs.CV].
URL https://arxiv.org/abs/1409.1556

  * [25] A. Hennequin, B. Couturier, V. Gligorov, S. Ponce, R. Quagliani, L. Lacassagne, A fast and efficient SIMD track reconstruction algorithm for the LHCb upgrade 1 VELO-PIX detector, Journal of Instrumentation 15 (06) (2020) P06018.  doi:https://doi.org/10.1088/1748-0221/15/06/p06018.
  * [26] ATLAS Collaboration, L0MuonTriggerPublicResults.
URL
https://twiki.cern.ch/twiki/bin/view/AtlasPublic/L0MuonTriggerPublicResults


# Reverse Engineered Diophantine Equations

Stevan Gajović

###### Abstract.

We answer a question of Samir Siksek, asked at the open problems session of
the conference “Rational Points 2022”, which, in a broader sense, can be
viewed as a reverse engineering of Diophantine equations. For any finite set
$S$ of perfect integer powers, using Mihăilescu’s theorem, we construct a
polynomial $f_{S}\in\mathbb{Z}[x]$ such that the set $f_{S}(\mathbb{Z})$
contains a perfect integer power if and only if it belongs to $S$. We first
discuss the easier case where we restrict to all powers with the same
exponent. In this case, the constructed polynomials are inspired by Runge’s
method and Fermat’s Last Theorem. Therefore we can construct a polynomial-
exponential Diophantine equation whose solutions are described in advance.

Key words: Diophantine equations, Fermat’s Last Theorem, Mihăilescu’s theorem,
Runge’s method, elliptic curves

##  1\. Introduction

Diophantine problems look innocent but often are tricky. They can be
formulated using very basic mathematics, but it turns out that solving many of
them took centuries of serious work by numerous mathematicians. Some of them
are still not solved; they are famous conjectures in mathematics and, more
precisely, in number theory.

The most famous Diophantine equations include some relations with perfect
integer powers, and we briefly survey them in §2.1. We will use some of these
equations in the main results of this paper in §3.2 and §4.

When we also include rational solutions, to simplify the notation, we will
talk about Diophantine problems. We also mention in §2.2 one instance of
famous Diophantine problems, classifying the set of rational points
$C(\mathbb{Q})$ on curves $C$ defined over $\mathbb{Q}$. We briefly use some
of the results and techniques mentioned there in §3.2.

In §2, we note that Diophantine problems are an active and very challenging
area of mathematics. On the other hand, there are equations that are trivial
to solve, such as linear Diophantine equations, or some equations that have
obvious ways to be solved, e.g., by considering them modulo $n$ for some
$n\in\mathbb{Z}_{>1}$, or, more generally, using local obstructions.

Therefore, it is an interesting problem to look into the middle case - can we
create similar Diophantine equations with a prescribed set of solutions so
that we can solve them more easily, but not that obviously? The main task of
this article is to take a different perspective - and construct Diophantine
equations of a specific shape with a described set of solutions. That was
precisely the question Samir Siksek asked during the conference “Rational
Points 2022”. Namely, if we consider a Diophantine equation of the type

$f(x)=y^{n},$

for a given $f\in\mathbb{Z}[x]$, with unknowns $x,y,n\in\mathbb{Z}$ and $n\geq
2$, can it happen that a triple $(x,y,n)$ is a solution if and only if $y^{n}$
belongs to a finite set of integer powers given in advance?

###  1.1. Main results

We now reformulate Siksek’s question and make it more precise. Let

$\mathcal{P}=\\{a^{m}\colon a\in\mathbb{Z},m\geq 2\\}\subseteq\mathbb{Z}$

denote the set of all perfect integer powers.

######  Question 1.1 (Siksek).

Let $S\subseteq\mathcal{P}$ be a finite set of perfect powers. Is there a
polynomial $f_{S}\in\mathbb{Z}[x]$ such that
$f_{S}(\mathbb{Z})\cap\mathcal{P}=S$?

Here we give an affirmative answer to Question 1.1 by constructing such a
polynomial in §4, see Theorem 4.1.

One could ask an easier question. Let $m\geq 2$ be a fixed integer. Denote by
$\mathcal{P}_{m}=\\{a^{m}\colon a\in\mathbb{Z}\\}$ the set of all perfect
$m$th powers of integers.

######  Question 1.2 (Siksek).

Let $S\subseteq\mathcal{P}_{m}$ be a finite set of $m$th powers. Is there a
polynomial $f_{S}\in\mathbb{Z}[x]$ such that
$f_{S}(\mathbb{Z})\cap\mathcal{P}_{m}=S$?

We recall that in Question 1.2, we fix $m\in\mathbb{Z}_{\geq 2}$. Curves of
the shape $C:y^{m}=f(x)$ are called superelliptic curves. We can also rephrase
Question 1.2 as a task to construct superelliptic curves whose integral points
have prescribed $y$-coordinates in advance.

We explain two different methods to construct a polynomial as asked for in
Question 1.2 (one method only works for $m\geq 3$) in §3.1 and §3.2, see
Theorem 3.1 and Corollary 3.3, respectively.

We briefly comment on the same question when integers are replaced by rational
numbers in section 5. We note that one of two approaches to solving Question
1.2 still works over rationals, as explained in §5.1. However, our approach
for Question 1.1 does not extend to rational numbers; but in the meantime,
Question 5.2 was solved by Santicola [20] by clever and precise refining of
the arguments presented here, see §5.2.

##  2\. Famous Diophantine problems

###  2.1. Famous Diophantine equations

Now we present several famous Diophantine equations or related conjectures. We
also mention the time needed to solve these equations. We recall that these
equations mentioned below ask for or are related to integer solutions.

  1. (1)

Fermat’s Last Theorem: $x^{n}+y^{n}=z^{n}$, with an integer $n\geq 3$. It took
more than 350 years until it was proven in series of papers by Wiles et al.
that this equation has only trivial solutions, i.e., such that $xyz=0$.

  2. (2)

Generalized Fermat’s Equation: $x^{k}+y^{l}=z^{m}$, where $k,l,m\geq 2$ are
integers, see [2] for more details. A special case of this equation is related
to the Beal conjecture which states that if

$k,l,m\geq 3,$

then for any solution $(x,y,z)\in\mathbb{Z}^{3}$, there is a prime number $p$
such that $p\mid x$, $p\mid y$, $p\mid z$, i.e., solutions cannot be triples
of coprime integers. Note that the condition that $k,l,m\geq 3$ is necessary,
as otherwise one can find coprime solutions, such as $1^{k}+2^{3}=(\pm
3)^{2}$, for any $k\geq 2$, or, for example, $2^{5}+7^{2}=3^{4}$. More
identities can be found in [2].

For fixed $k,l,m\geq 2$ such that

$1/k+1/l+1/m<1,$

Darmon and Granville [6] proved that there are only finitely many triples of
coprime integers $(x,y,z)$ such that $x^{k}+y^{l}=z^{m}$.

If we vary $k,l,m\geq 2$ such that $1/k+1/l+1/m<1$, there is a conjecture, the
Fermat–Catalan conjecture, stating that there are only finitely many sextuples

$(x,y,z,k,l,m)$

such that $x^{k}+y^{l}=z^{m}$ and $x,y,z$ are coprime.

  3. (3)

Former Catalan’s conjecture, now Mihăilescu’s theorem: The equation
$x^{a}-y^{b}=1$ with integers $a,b\geq 2$ and $x,y>0$ has only one solution
$3^{2}-2^{3}=1$. In other words, the only two positive integers which are
consecutive perfect powers of integers are 8 and 9. This statement was
conjectured by Catalan and proved by Mihăilescu [14] slightly more than 150
years later.

  4. (4)

(Generalized) Ramanujan-Nagell equation: The Ramanujan-Nagell equation is the
equation $x^{2}+7=2^{n}$, where $x$ and $n$ are integers. Ramanujan
conjectured that there are five values of $n$ for which the equation has a
solution, $n\in\\{3,4,5,7,15\\}$, and Nagell proved this 35 years later. It
was first published in Norwegian, and later in English [18]. The generalized
Ramanujan-Nagell equation is an equation of the shape $x^{2}+D=y^{n}$, where
$x,n,D\in\mathbb{Z}$ and $n\geq 3$. It is widely studied, for example, see
[5], or [11] for a recent survey on this equation. Note that some authors,
such as, in [16], consider the slightly different equation
$F(x)=p_{1}^{e_{1}}\cdots p_{s}^{e_{s}}$, for some $F\in\mathbb{Z}[x]$ and
$p_{1},...,p_{s}$ prescribed primes and $e_{i}\geq 0$ (unknown) exponents.

###  2.2. Rational points on curves

One of the most important examples of Diophantine problems is the following
trichotomy of rational points on curves defined over $\mathbb{Q}$ (and in
fact, over any number field $K/\mathbb{Q}$). Let $C/\mathbb{Q}$ be a nice
curve; here, the adjective nice is a well-known notation used for a smooth,
projective, and geometrically irreducible curve. Denote the genus of $C$ by
$g(C)$.

  1. (a)

(known for a long time) If $g(C)=0$, then either $C(\mathbb{Q})=\emptyset$
(e.g., consider $C:x^{2}+y^{2}=3$) or $C(\mathbb{Q})$ is infinite (e.g.,
consider $C:x^{2}+y^{2}=1$, this curve can be used to parametrize the
Pythagorean triples) and furthermore isomorphic to
$\mathbb{P}^{1}(\mathbb{Q})$.

  2. (b)

(Mordell, 1922 [15]) If $g(C)=1$ and $C(\mathbb{Q})\neq\emptyset$, then $C$ is
called an elliptic curve and there is a group law on $C(\mathbb{Q})$, which
makes $C(\mathbb{Q})$ a finitely generated group, i.e.,
$C(\mathbb{Q})\cong\mathbb{Z}^{r}\bigoplus T$, where $T$ is a finite (torsion)
subgroup, and $r$ is called the rank of $C$ over $\mathbb{Q}$.

There are only finitely many possibilities for the torsion subgroup $T$
(exactly 15), as these were classified by Mazur [12], [13]. To prove this
statement, Mazur determined rational points on certain types of curves, called
modular curves.

Computing the rank $r$ is still a difficult problem, and there are some ways
that might succeed in computing it, but there is still no guarantee that the
known algorithms can compute the rank of all elliptic curves, look, for
example, at Silverman’s book [23]. Also, these algorithms for computing ranks
are quite complicated, and in practice, it is challenging for humans to
perform them. Hence, these algorithms are implemented in computer algebra
systems, such as Magma [3].

  3. (c)

(Faltings, 1983 [8]) If $g(C)\geq 2$, then $C(\mathbb{Q})$ is finite. In [15],
Mordell conjectured this statement, so it took about 60 years until it was
proved. This statement was so difficult and significant that Faltings won a
Fields medal for this proof. By now, there are a few different proofs of this
statement.

However, there are no practical algorithms to compute $C(\mathbb{Q})$ for a
given curve $C$, and it is a very active area of research to find methods that
can compute precise sets of rational points on curves. As we have already seen
in the case of elliptic curves, there are significant problems that reduce to
computing rational points on curves (e.g., Mazur’s theorem), so it is indeed
very important nowadays to further develop existing methods for determining
rational points on curves.

Unlike elliptic curves, the set $C(\mathbb{Q})$ for $g(C)\geq 2$ has no
particular algebraic structure, so sometimes we want to study it by embedding
it into an object with more structure, called the Jacobian $J$ of the curve
$C$; this is an abelian variety. As for elliptic curves, Weil [26] in 1929
proved, now called Mordell-Weil theorem, that $J(\mathbb{Q})$ is a finitely
generated abelian group, i.e., $J(\mathbb{Q})\cong\mathbb{Z}^{r}\bigoplus T$,
where $T$ is a finite (torsion) subgroup and $r$ is called the rank of $J$. In
contrast to elliptic curves, much less is known about possibilities for $T$,
and the computation of $r$ is much more difficult. In some cases, there are
ways to do so, for example, for Jacobians of hyperelliptic curves, by Stoll
[25], which is implemented in Magma.

##  3\. The same exponent

Fix $m\geq 2$. Let $S=\\{a_{1}^{m},\ldots,a_{k}^{m}\\}$ be a set of $k$
distinct $m$th powers of integers $a_{1},\ldots,a_{k}$. We now construct a
polynomial $f_{S}\in\mathbb{Z}[x]$ such that
$f_{S}(\mathbb{Z})\cap\mathcal{P}_{m}=S$.

###  3.1. First approach

The solution is inspired by Runge’s method (invented in 1887 by Runge [19];
see [21, Chapter 5] for a nice survey). We first define an auxiliary
polynomial

$g(x)=(x-a_{1})\cdots(x-a_{k}).$

Consider a polynomial

$f_{S}(x)=(x(x^{2}+1)g(x))^{4m}+(x^{2m}-x^{2}+2)g(x)^{2m}+x^{m}.$

######  Theorem 3.1.

We have

  1. (i)

$S\subseteq f_{S}(\mathbb{Z})\cap\mathcal{P}_{m}$;

  2. (ii)

$f(x)\notin\mathcal{P}_{m}$ for any
$x\in\mathbb{Z}\setminus\\{0,a_{1},\ldots,a_{k}\\}$;

  3. (iii)

$f_{S}(\mathbb{Z})\cap\mathcal{P}_{m}=S$.

###### Proof.

  1. (i)

We evaluate (note that this statement is true regardless whether $a_{1}\cdots
a_{k}=0$ or not)

$f_{S}(0)=2(a_{1}\cdots
a_{k})^{2m},\;\;g(a_{i})=0,\;\;f_{S}(a_{i})=a_{i}^{m},\;\;\text{for $1\leq
i\leq k$}.$

  2. (ii)

Let $x\in\mathbb{Z}\setminus\\{0,a_{1},\ldots,a_{k}\\}$. We first note that
$|g(x)|=|(x-a_{1})\cdots(x-a_{k})|\geq 1$. Then

$(x^{2m}-x^{2}+2)g(x)^{2m}\geq x^{2m}-x^{2}+2\geq x^{2m}-|x|^{m}+2>|x|^{m}.$

It follows that $(x^{2m}-x^{2}+2)g(x)^{2m}+x^{m}>0$, so

$f_{S}(x)>(x(x^{2}+1)g(x))^{4m}=((x(x^{2}+1)g(x))^{4})^{m}.$

Now we prove that

(1) $f_{S}(x)<((x(x^{2}+1)g(x))^{4}+1)^{m},$

implying that $f(x)$ cannot be an $m$th power because it is nested between two
consecutive $m$th powers. To prove inequality (1), it suffices to show

(2) $m(x(x^{2}+1)g(x))^{4m-4}>(x^{2m}-x^{2}+2)g(x)^{2m}+x^{m}.$

The following two inequalities will help us to show (2). The first one is

(3) $(x(x^{2}+1)g(x))^{4m-4}>(x^{2m}-x^{2}+2)g(x)^{2m}$

which holds because $g(x)^{4m-4}\geq g(x)^{2m}$ and the inequality

$(x(x^{2}+1))^{4m-4}>x^{2m}-x^{2}+2$

is trivial (note that $x\neq 0$).

The second one is clear

(4) $(x(x^{2}+1)g(x))^{4m-4}\geq(x(x^{2}+1))^{4m-4}>|x|^{m}\geq x^{m}.$

Inequality (2) follows from $m\geq 2$ and inequalities (3) and (4).

  3. (iii)

Follows directly from (i) and (ii).

∎

###  3.2. Second approach

If $m\geq 3$, there is another way to construct the required polynomial
$f_{S}\in\mathbb{Z}[x]$. The term
$(x^{2m}-x^{2}+2)((x-a_{1})\cdots(x-a_{k}))^{2m}$ in the previous construction
was used to ensure that $f(0)$ in not an $m$th power if $a_{1}\cdots a_{k}\neq
0$. One could try a simpler construction

$g_{S}(x)=((x-a_{1})\cdots(x-a_{k}))^{m}+x^{m}.$

By Fermat’s Last Theorem, $g_{S}(x)$ is not an $m$th power unless
$x\in\\{0,a_{1},\ldots,a_{k}\\}$. However, we want to exclude the possibility
for $x=0$ if $a_{1}\cdots a_{k}\neq 0$. We instead use a more general Fermat’s
Equation, as suggested by Samir Siksek. Consider

$f_{S}(x)=3((x-a_{1})\cdots(x-a_{k}))^{m}+x^{m}.$

######  Lemma 3.2.

Let $m\geq 3$. If $(x_{1},x_{2},x_{3})\in\mathbb{Z}^{3}$ satisfy
$3x_{1}^{m}+x_{2}^{m}=x_{3}^{m}$, then $x_{1}=0$.

###### Proof.

We distinguish three cases according to whether $m$ has a prime divisor $p\geq
5$, $m$ is a power of 3, or $m\geq 4$ is a power of 2.

  1. (i)

Let $p\geq 5$ be a prime number dividing $m$. By a generalized approach to
Fermat’s Last Theorem by Kraus [10], the equation
$3x_{1}^{p}+x_{2}^{p}=x_{3}^{p}$ has only trivial solutions
$x_{1}x_{2}x_{3}=0$, which implies that $x_{1}=0$. See also notes by Siksek
[24] for a nice explanation of the modular method. In this concrete case, see
[24, Theorems 1, 15] for an argument that a non-trivial solution to
$3x_{1}^{p}+x_{2}^{p}=x_{3}^{p}$ corresponds to a certain newform of weight 2
and level 6, which does not exist. Hence, if $3x_{1}^{m}+x_{2}^{m}=x_{3}^{m}$,
then $x_{1}=0$.

  2. (ii)

If $m$ is a power of 3, it suffices to consider $m=3$. The cubic curve

$X\colon 3x_{1}^{3}+x_{2}^{3}=x_{3}^{3}$

has genus one and $X(\mathbb{Q})\neq\emptyset$, hence it is isomorphic over
$\mathbb{Q}$ to an elliptic curve. As pointed out in §2.2 (b), using Magma
[3], we prove that the rank of $X$ is zero, and that $X(\mathbb{Q})$ consists
only of the point at infinity, which corresponds precisely to $x_{1}=0$.

  3. (iii)

If $m\geq 4$ is a power of 2, it suffices to consider $m=4$. Integral
solutions to $3x_{1}^{4}+x_{2}^{4}=x_{3}^{4}$ correspond to the rational
points on a curve

$X^{\prime}\colon 3x^{4}+y^{4}=1.$

One can consider a curve

$X/\mathbb{Q}\colon y^{2}=1-3x^{4}.$

As we stated in §2.2 (c), we can embed $X$ into its Jacobian $J$, and, using
Magma, we prove $J(\mathbb{Q})\simeq\mathbb{Z}/2\mathbb{Z}$. Hence, we
conclude $\\#X(\mathbb{Q})\leq\\#J(\mathbb{Q})=2$, so $X(\mathbb{Q})=\\{(0,\pm
1)\\}$, implying $X^{\prime}(\mathbb{Q})=\\{(0,\pm 1)\\}$, and hence, that
$x_{1}=0$.

∎

######  Corollary 3.3.

We have $f_{S}(\mathbb{Z})\cap\mathcal{P}_{m}=S$.

###### Proof.

It follows by Lemma 3.2 that $f_{S}(x)\in\mathcal{P}_{m}$ only when $x\in S$,
and we compute directly that $f_{S}(a_{i})=a_{i}^{m}$ for every $1\leq i\leq
k$. ∎

We note that this approach does not work directly for $m=2$ as the equation

$qx_{1}^{2}+x_{2}^{2}=x_{3}^{2}$

has solutions with $x_{1}\neq 0$ for any $q\in\mathbb{Z}$. We may assume that
$q>0$. If $q$ is not a square of an integer, then, consider $x_{2}=1$ and note
that Pell’s equation $x_{3}^{2}-qx_{1}^{2}=1$ has infinitely many solutions.
If $q=r^{2}$ is a square of an integer $r$, then $(rx_{1},x_{2},x_{3})$ is a
Pythagorean triple, and we can find infinitely many of them (for example, let
$x_{1}=2s$, $x_{2}=s^{2}-r^{2}$, $x_{3}=s^{2}+r^{2}$, for $s\in\mathbb{Z}$).

##  4\. General case

In this section, using Mihăilescu’s theorem, we construct a polynomial with
the desired property in the general case. Let
$S=\\{b_{1},\ldots,b_{k}\\}\subseteq\mathcal{P}$ be a finite set of perfect
integer powers. We construct a polynomial $f_{S}\in\mathbb{Z}[x]$ such that
$f_{S}$ is the identity on $S$, and $f_{S}(x)\notin\mathcal{P}$, for all
$x\in\mathbb{Z}\setminus S$. Define an auxiliary polynomial

$g(x)=((x-b_{1})\cdots(x-b_{k}))^{4}+1,$

and let

$f_{S}(x)=g(x)((x-1)g(x)+1).$

We now prove that $f_{S}$ satisfies the property asked in Question 1.1, hence
giving a positive answer to Question 1.1.

######  Theorem 4.1.

  1. (i)

Let $x\in\mathbb{Z}$. Then $f_{S}(x)\in\mathcal{P}$ if and only if $x\in S$.

  2. (ii)

We have $f_{S}(\mathbb{Z})\cap\mathcal{P}=S$.

###### Proof.

  1. (i)

Let $x\in\mathbb{Z}$ be such that $f_{S}(x)=y^{n}$, for some integers $y,n$
with $n\geq 2$. Since $g(x)$ and $(x-1)g(x)+1$ are coprime integers, we
conclude that there is $z\in\mathbb{Z}$ such that $g(x)=\pm z^{n}$, and since
$g(x)>0$, we may assume that $g(x)=z^{n}$. Denote

$c\colonequals(x-b_{1})\cdots(x-b_{k}).$

Then we have

$z^{n}-c^{4}=1.$

Since the exponent of $c$ is greater than 3, by Mihăilescu’s theorem, the only
possibility is that $z=1$ and $c=0$, implying that $x=b_{i}$ for some $1\leq
i\leq k$. If $x\in S$, we evaluate $f_{S}(x)=x\in\mathcal{P}$.

  2. (ii)

Follows from (i) and the fact that $f_{S}$ is the identity on $S$.

∎

##  5\. Generalizations.

It is interesting to see whether our approaches work when we ask the same
questions with integers replaced by rational numbers. Denote by
$\mathcal{Q}=\\{a^{m}\colon a\in\mathbb{Q},m\geq 2\\}$ the set of all rational
powers, and, for a fixed integer $m\geq 2$, let
$\mathcal{Q}_{m}=\\{a^{m}\colon a\in\mathbb{Q}\\}$ be the set of all $m$th
powers of rational numbers. We have naturally two questions.

###  5.1. Rational numbers: The same exponent

######  Question 5.1.

Let $S\subseteq\mathcal{Q}_{m}$ be a finite subset of $m$th rational powers.
Is there a polynomial $f_{S}\in\mathbb{Q}[x]$ such that
$f_{S}(\mathbb{Q})\cap\mathcal{Q}_{m}=S$?

We first note that for $m\geq 3$ we can give a positive answer to Question
5.1. It is clear that the first approach cannot be used because there is no
version of Theorem 3.1 that covers rational numbers. In this approach, we use
the property that distinct integers differ by at least 1, which is not true
for rational numbers; their absolute difference can be arbitrarily small.
However, Lemma 3.2 remains true for rational numbers, i.e., if we replace the
condition $(x_{1},x_{2},x_{3})\in\mathbb{Z}^{3}$ by
$(x_{1},x_{2},x_{3})\in\mathbb{Q}^{3}$. This follows because the equation
$3x_{1}^{m}+x_{2}^{m}=x_{3}^{m}$ is homogeneous, so any rational solution
easily leads to an integer solution. Hence, if
$S=\\{a_{1}^{m},\ldots,a_{k}^{m}\\}$ is a set of $k$ distinct $m$th powers of
rational numbers $a_{1},\ldots,a_{k}$, where $m\geq 3$, then again, for the
polynomial

$f_{S}(x)=3((x-a_{1})\cdots(x-a_{k}))^{m}+x^{m},$

we have that $f_{S}(\mathbb{Q})\cap\mathcal{Q}_{m}=S$.

###  5.2. Rational numbers: General case

On the other hand, the approach from section 4 does not work directly over
rational numbers. We cannot use the coprimality argument in the factorization
so easily; we would need to take care of possible denominators. Hence, we can
formulate the question:

######  Question 5.2.

Let $S\subseteq\mathcal{Q}$ be a finite set of perfect rational powers. Is
there a polynomial $f_{S}\in\mathbb{Q}[x]$ such that
$f_{S}(\mathbb{Q})\cap\mathcal{Q}=S$?

This question was answered affirmatively by Santicola [20], who noted that it
is sufficient to use special cases of Mihăilescu’s theorem proven by Lebesgue
[9] to answer Question 1.2. Furthermore, Santicola’s construction uses the
results of [1, 4, 7]. To handle possible denominators and to prove the key
lemma [20, Lemma 6], Santicola used the result of Pethő [17], and
independently of Shorey and Stewart [22].

###  5.3. Challenge

We see that this article has already inspired further research. It is also
interesting to consider this question over other rings and fields. We
challenge the interested reader to try to answer the analogous questions over
$\mathbb{Z}[i]$ or $\mathbb{Q}[i]$, or, in general, over $\mathcal{O}_{K}$ or
$K$, where $K/\mathbb{Q}$ is any number field.

### Acknowledgment.

The author would like to thank Samir Siksek for raising the questions,
discussing the topic and giving useful comments and suggestions, Michael Stoll
for organizing the conference “Rational Points 2022” and for the financial
support for the conference, Ludwig Fürst and Himanshu Shukla for conversations
about this topic and writing a report on the problem session, and Pieter
Moree, Oana Padurariu, and Alain Togbe for many valuable comments that
improved the article. The author was supported by a guest postdoc fellowship
at the Max Planck Institute for Mathematics in Bonn, by Czech Science
Foundation GAČR, grant 21-00420M, and by Charles University Research Centre
program UNCE/SCI/022 during various stages of this project.

## References

  * [1] Bennett, M. and Ellenberg, J. S. and Ng, N. (2010). The Diophantine Equation $A^{4}+2^{\delta}B^{2}=C^{n}$, International Journal of Number Theory. 06: 311–338.
  * [2] Bennett, M. and Mihăilescu, P. and Siksek, S. (2016). The generalized Fermat equation. Open problems in mathematics: Springer, [Cham]. 173–205.
  * [3] Bosma, W. and Cannon, J. and Playoust, C. (1997). The Magma algebra system. I. The user language. Journal of Symbolic Computation. 24(3–4): 235–265.
  * [4] Bruin, N. (1999). The Diophantine Equations $x^{2}±\pm y^{4}=\pm z^{6}$ and $x^{2}+y^{8}=z^{3}$. Compositio Mathematica. 118(3): 305–321.
  * [5] Bugeaud, Y. and Mignotte, M. and Siksek, S. (2006). Classical and modular approaches to exponential Diophantine equations. II. The Lebesgue-Nagell equation. Compositio Mathematica. 142(1): 31–62.
  * [6] Darmon, H. and Granville, A. (1995). On the equations $z^{m}=F(x,y)$ and $Ax^{p}+By^{q}=Cz^{r}$. The Bulletin of the London Mathematical Society. 27(6): 513–543.
  * [7] Ellenberg, J. S. (2004). Galois Representations Attached to $\mathbb{Q}$-Curves and the Generalized Fermat Equation $A^{4}+B^{2}=C^{p}$. American Journal of Mathematics. 126(4): 763–787.
  * [8] Faltings, G. (1983). Endlichkeitssätze für abelsche Varietäten über Zahlkörpern. Inventiones Mathematicae. 73(3): 349–366.
  * [9] Lebesgue, V. A. (1850). Sur l’impossibilité, en nombres entiers, de l’équation $x^{m}=y^{2}+I$. Nouvelles annales de mathématiques : journal des candidats aux écoles polytechnique et normale 9: 178–181.
  * [10] Kraus, A. (1997). Majorations effectives pour l’équation de Fermat généralisée. Journal Canadien de Mathématiques. 49(6): 1139–1161.
  * [11] Le, M. and Soydan, G. (2020). A brief survey on the generalized Lebesgue-Ramanujan-Nagell equation. Surveys in Mathematics and its Applications. 15: 473–523.
  * [12] Mazur, B. (1977). Modular curves and the Eisenstein ideal, Publications Mathématiques I.H.E.S. 47: 33–186.
  * [13] Mazur, B. (1977). Rational points on modular curves. In: Serre, JP., Zagier, D.B. (eds) Modular Functions of one Variable V. Lecture Notes in Mathematics, vol 601. Springer, Berlin, Heidelberg. 107–148.
  * [14] Mihăilescu, P. (2004). Primary cyclotomic units and a proof of Catalan’s conjecture. Journal für die Reine und Angewandte Mathematik. [Crelle’s Journal]. 572: 167–195.
  * [15] Mordell, L, J. (1922), On the rational solutions of the indeterminate equations of the third and fourth degrees. Proc. Camb. Phil. Soc. 21: 179–192.
  * [16] Moree, P. and Stewart, C. L. (1990) Some Ramanujan-Nagell equations with many solutions, Indag. Math. (N. S.) (4) 1: 465–472.
  * [17] Pethő, A. (1982). Perfect powers in second order linear recurrences. Journal of Number Theory. 15(1): 5–13.
  * [18] Nagell, T. (1961). The diophantine equation $x^{2}+7=2^{n}$. Arkiv för Matematik. Journal Canadien de Mathématiques. 4: 185–187.
  * [19] Runge, C. (1887). Ueber ganzzahlige Lösungen von Gleichungen zwischen zwei Veränderlichen. Journal für die Reine und Angewandte Mathematik. [Crelle’s Journal]. 100: 425–435.
  * [20] Santicola, K. (2022). Reverse engineered Diophantine equations over $\mathbb{Q}$. Journal de Théorie des Nombres de Bordeaux, to appear.
  * [21] Schoof, R. (2008). Catalan’s Conjecture. Springer-Verlag London, Ltd., London.
  * [22] Shorey, T. N. and Stewart, C. L. (1983). On the Diophantine equation $ax^{2t}+bx^{t}y+cy^{2}=d$ and pure powers in recurrence sequences. Mathematica Scandinavica 52(1): 24–36.
  * [23] Silverman, J. H. (2009), The arithmetic of elliptic curves, 2nd edition. Springer: Graduate Texts in Mathematics. Volume 106.
  * [24] Siksek, S. (2016). The Modular Approach to Diophantine Equations, http://homepages.warwick.ac.uk/~maseap/sarajevo/notes.pdf. (accessed 18 May 2022)
  * [25] Stoll, M. (2001). Implementing 2-Descent for Jacobians of Hyperelliptic Curves. Acta Arithmetica. 98(3): 245–277.
  * [26] Weil, A. (1929). L’arithmétique sur les courbes algébriques, Acta Mathematica. 52(1): 281–315.


[a,b,c]S. Prelovsek

# Charmonium-like resonances in coupled $D\bar{D}$ \- $D_{s}\bar{D}_{s}$
scattering

S. Collins  D. Mohler  M. Padmanath  S. Piemonte

###### Abstract

Charmonium-like resonances and bound states with isospin zero and
$J^{PC}=0^{++},~{}1^{--},~{}2^{++},~{}3^{--}$ are extracted on the lattice.
Coupled $D\bar{D}$ and $D_{s}\bar{D}_{s}$ scattering suggests three
charmonium-like states with $J^{PC}=0^{++}$ in addition to $\chi_{c0}(1P)$: a
so far unobserved $D\bar{D}$ bound state just below threshold, a conventional
resonance likely related to $\chi_{c0}(3860)/\chi_{c0}(2P)$ and a narrow
resonance just below the $D_{s}\bar{D}_{s}$ threshold with a large coupling to
$D_{s}\bar{D}_{s}$ likely related to $X(3915)/\chi_{c0}(3930)$. One-channel
$D\bar{D}$ scattering renders resonances and bound states with
$J^{PC}=1^{--},~{}2^{++},~{}3^{--}$ related to the observed conventional
charmonia. Lattice QCD ensembles from the CLS consortium with $m_{\pi}\simeq
280$ MeV are utilized.

##  1 Introduction

We study charmonium-like resonances with isospin zero by simulating the
scattering on the lattice. Studies of charmonium-like states are of particular
interest as most of the exotic hadrons discovered contain a $\bar{c}c$ pair.
Furthermore, all these exotic hadrons are resonances and can decay strongly.
We consider the resonances and bound states with
$J^{PC}=0^{++},~{}1^{--},~{}2^{++},~{}3^{--}$ indicated in Fig. 1, where the
relevant thresholds are also shown. The scalars are extracted from the coupled
channel scattering $D\bar{D}$ \- $D_{s}\bar{D}_{s}$ since two resonances
reside near the $D_{s}\bar{D}_{s}$ threshold. The negative parity states of
interest lie significantly below the $D_{s}\bar{D}_{s}$ threshold and are
explored with one-channel $D\bar{D}$ scattering. The spin-two resonance is
also extracted in the approximation of one-channel $D\bar{D}$ scattering.
Details are provided in publications [1] and [2] which consider the negative
and positive parity sectors, respectively.

Figure 1:  Charmonium-like states from PDG [3] are shown in black. The
ellipses indicate states considered in the present lattice study:
$J^{PC}=0^{++},~{}2^{++}$ and $J^{PC}=1^{--},~{}3^{--}$ were extracted in [2]
and [1] respectively. The state indicated in blue is the lattice prediction of
yet-unobserved unconventional state that couples strongly to $D\bar{D}$ and
lies near this threshold. Other states have been experimentally discovered:
all of them seem to be conventional charmonia, except for
$X(3915)/\chi_{c0}(3930)$ which is unconventional state that couples strongly
to $D_{s}\bar{D}_{s}$ according to this lattice study.

##  2 Extracting scattering amplitudes from lattice

The infinite-volume scattering matrices $t_{ij}(E_{cm})$ are determined from
the finite-volume eigen-energies via Lüscher’s formalism. The eigen-energies
are obtained from the correlation matrices based on a number of $\bar{c}c$,
$D\bar{D}$, $D_{s}\bar{D}_{s}$ and $J/\psi\omega$ interpolating fields with
appropriate quantum numbers. The channel $J/\psi\omega$ is treated as
decoupled, while the channel $\eta_{c}\eta$ is omitted from the simulation. We
consider the systems with total momenta
$|\vec{P}|=0,~{}2\pi/L,~{}\sqrt{2}(2\pi/L)$ on lattices with $L=24a,~{}32a$ in
order to constrain the scattering matrix at more values of the center-of-
momentum energy $E_{cm}$. The $c\bar{c}$ annihilation is omitted like in most
of the charmonium studies on the lattice, while other Wick contractions as
evaluated using the distillation method. The simulation is performed on the
$N_{f}=2+1$ ensembles U101 and H105 generated by the CLS consortium [4] with
$a=0.08636(98)(40)$ fm and unphysical $m_{u/d}>m_{u/d}^{exp}$ and
$m_{s}<m_{s}^{exp}$, such that $2m_{u/d}+m_{s}$ is close to the physical
value. The value of the charm quark mass (related to $\kappa_{c}=0.12315$) is
slightly higher than physical and the resulting masses of the relevant stable
hadrons are collected in Table 1. In particular, the $D\bar{D}$ and
$D_{s}\bar{D}_{s}$ thresholds are closer to each other than in experiment as
shown in Fig. 3.

The energy dependence of scattering matrices $t_{ij}(E_{cm})$ is parametrized
and the parameters are fitted as to best satisfy the Lüscher quantization
condition for all lattice eigen-energies simultaneously. The TwoHadronsInBox
package [5] is utilised for this. The typical parametrizations used for the
real part of the inverse scattering amplitude are of Breit-Wigner type
$\tilde{K}^{-1}/E_{cm}=p^{2l+1}\cot\delta/E_{cm}=a+bE_{cm}^{2}$ and
$\tilde{K}_{ij}^{-1}/E_{cm}=a_{ij}+b_{ij}E_{cm}^{2}$ for one-channel and two-
channel scattering, respectively, with more details given in [2, 1]. The study
makes several simplifying assumptions (detailed in Section 5 of [2]) necessary
for a first investigation of this coupled-channel system.

$m_{\pi}$ [MeV] |  $m_{K}$ [MeV] |  $m_{D}$ [MeV] |  $m_{D^{*}}$ [MeV] |  $m_{D_{s}}$ [MeV] |  $M_{av}$ [MeV]
---|---|---|---|---|---
280(3) | 467(2) | 1927(2) | 2050(2) | 1981(1) | 3103(3)
Table 1: Hadron masses for the gauge configurations used in this study, where
$M_{av}=(m_{\eta_{c}}+3m_{J/\psi})/4$.  Figure 2:  Pole singularities of the
scattering amplitude/matrix in the complex energy plane that are extracted
from this lattice study. These singularities are associated with the
charmonium-like states and various symbol shapes denote the corresponding
$J^{PC}$. The $3^{--}$ state is not presented as we did not extract its width.

(a) Left pane: The masses $m$ obtained from our lattice study with unphysical
quark masses $m_{u/d}>m_{u/d}^{exp}$, $m_{s}<m_{s}^{exp}$ and $m_{c}\gtrsim
m_{c}^{exp}$. The green lines denote positions of the $D\bar{D}$ and
$D_{s}\bar{D}_{s}$ thresholds on our lattice with $m_{D}\simeq 1927~{}$MeV and
$m_{Ds}\simeq 1981~{}$MeV.

(b) Left pane: The same masses $m$ as above, but shifted to
$m-E^{ref}+E^{ref}_{exp}$ in order to account for the dominant effect of
unphysical quark masses in the simulation. The reference energy is
$E^{ref}=2m_{D}$ ($2m_{D_{s}}$) for the state closest to the $D\bar{D}$
($D_{s}\bar{D}_{s}$) threshold, while
$E^{ref}=M_{av}=\tfrac{1}{4}(3m_{J/\psi}+m_{\eta_{c}})$ for the remaining four
states. The green lines denote experimental thresholds.

Figure 3:  Masses of charmonium-like states with isospin zero from a lattice
simulation [2, 1] (left) compared to experiment (right). Left: The magenta
symbols correspond to hadrons extracted via the scattering analysis on the
lattice: diamonds represent resonances and triangles represent bound states.
The blue crosses are extracted directly from the lattice energies. Right:
experimental spectrum from the PDG [3], where $\chi_{c0}(3930)$ [6] and
$X(3915)$ are now identified as the same state.

##  3 Summary of the resulting hadrons and their relation to experiment

Let us summarize the properties of the charmonium-like hadrons found in this
simulation. All results were obtained at unphysical quark masses and are not
extrapolated to the continuum limit. The resonances and bound states are
related to the poles of the scattering matrix. The locations of the poles in
the complex energy plane are given in Fig. 2. The masses of the charmonium-
like states correspond to $\mathrm{Re}(E_{cm}^{pole})$ and are compared to the
experimental spectrum in Fig. 3. The resonance decay widths $\Gamma$ are given
by $2\mathrm{Im}(E_{cm}^{pole})$. We don’t compare them directly to experiment
since they depend on the phase space and therefore on the position of the
threshold, which is different in the simulation. Figure 4 shows the coupling
$g$ that parametrizes the full width of a resonance and compares it to
experiment

$\Gamma\equiv g^{2}p_{D}^{2l+1}/m^{2}~{}\quad\mathrm{with}\quad
l=0,1,2,3\quad\mathrm{for}\ J^{PC}=0^{++},1^{--},2^{++},3^{--}~{}.$ (1)

Due to the unphysical quark masses in the simulation, the difference
$m-E^{ref}$ is compared to experiment below and $m-E^{ref}+E^{ref}_{exp}$ is
shown in Fig. 3b. The reference energy $E^{ref}$ is the spin-averaged
charmonium mass $M_{av}=\tfrac{1}{4}(3m_{J/\psi}+m_{\eta_{c}})$ for
conventional candidates and a nearby threshold for states dominated by
$D_{(s)}\bar{D}_{(s)}$.

Below we review the spectroscopic properties of the charmonium-like states,
which are summarized also in Figs. 2, 3 and 4:

Figure 4:  The coupling $g$ that parametrizes the decay width for various
resonances $R\to D\bar{D}$ via $\Gamma\equiv g^{2}p_{D}^{2l+1}/m^{2}$ from
this lattice simulation and from experiment. The quantum numbers of $R$ and
partial waves $l$ in the decay $R\to D\bar{D}$ are provided at the bottom.

  * •

Ground states $\bm{J/\psi,~{}\psi(2S),~{}\chi_{c0}(1P),~{}\bm{\chi_{c2}(1P)}}$

These states lie significantly below the $D\bar{D}$ threshold and their masses
are extracted from the energies $m\\!=\\!E(P\\!=\\!0)$, resulting in a
reasonable agreement with experimental masses.

  * •

$\bm{0^{++}}$ state dominated by $\mathbf{D\bar{D}}$ slightly below threshold

The $D\bar{D}$ scattering near threshold in Fig. 5 indicates the presence of a
shallow $D\bar{D}$ bound state with the binding energy
$m-2~{}m_{D}=-4.0^{~{}+3.7}_{~{}-5.0}$ MeV. Our results suggest this state is
likely not a conventional charmonium $\bar{c}c$, but owes its existence to a
large interaction in the $D\bar{D}$ channel near threshold. Such a state has
not been claimed by experiments (yet). It would feature as a sharp increase of
the $D\bar{D}$ rate just above the threshold if it has a binding energy of a
few MeV. Various strategies for its experimental search were proposed by E.
Oset et al. and are listed in Section 8.1 of [2].

Figure 5:  $D\bar{D}$ scattering in partial wave $l\\!=\\!0$, where the
threshold position is shown in green. Left: The bound state occurs at the
energy where the red and orange curves intersect. Right: The shallow bound
state is responsible for a sharp rise of the $D\bar{D}$ rate
$N_{D\bar{D}}\propto p\sigma\propto\rho|t|^{2}$ ($\rho=2p/E_{cm}$ and $t$ is
the scattering amplitude).

The existence of a shallow $D\bar{D}$ bound state dubbed $X(3720)$ was already
suggested by an effective phenomenological model based on the exchanges of
light pseudoscalar and vector mesons [7]111This state with $m\simeq
3.718~{}$GeV is listed in Table 4 of Ref. [7].. In a molecular picture, a
$0^{++}$ state is expected as a partner of $X(3872)$ via heavy-quark symmetry
[8, 9, 17]. A similar virtual bound state with a binding energy of 20 MeV
follows from the data of the only previous lattice simulation of $D\bar{D}$
scattering [10]222The presence of this state was not mentioned in Ref. [10],
as such virtual bound states were not searched for..

  * •

$\bm{2^{++}}$ resonance and relation to $\bm{\chi_{c2}(3930)}$

We find a resonance with $J^{PC}\\!=\\!2^{++}$ $D\bar{D}$ scattering with
$l=2$ that is most likely related to the conventional $\chi_{c2}(3930)$ aka
$\chi_{c2}(2P)$ discovered by Belle [11]

$\displaystyle\mathrm{lat}:\ \ $ $\displaystyle
m-M_{av}=904^{~{}+14}_{~{}-22}~{}\mathrm{MeV}~{},\quad g=4.5^{+0.7}_{-1.5}\
\mathrm{GeV^{-1}}~{}$ (2) $\displaystyle\mathrm{exp}\ \chi_{c2}(3930):\ \ $
$\displaystyle m-M_{av}=854\pm 1~{}\mathrm{MeV}~{},\quad g=2.65\pm
0.12~{}\mathrm{GeV^{-1}}~{}.$

The masses are reasonably close, while the coupling from lattice QCD is larger
than that measured experiment, but not inconsistent given the statistical and
systematic uncertainties.

  * •

Broad $\bm{0^{++}}$ resonance and relation to $\bm{\chi_{c0}(3860)}$

The coupled $D\bar{D}-D_{s}\bar{D}_{s}$ scattering shows a broad peak in the
$D\bar{D}$ channel featured in the left pane in Fig. 6). It is related to a
resonance that couples mostly to $D\bar{D}$. We compare it to
$\chi_{c0}(3860)$ discovered by Belle [12] that is a candidate for the
conventional $\chi_{c0}(2P)$

$\displaystyle\mathrm{lat}:\ \ $ $\displaystyle
m-M_{av}=880^{+28}_{-20}~{}\mathrm{MeV}\ ,\quad g=1.35~{}^{+0.04}_{-0.08}\
\mathrm{GeV}~{}$ (3) $\displaystyle\mathrm{exp}\ \chi_{c0}(3860):\ $
$\displaystyle m-M_{av}=793~{}^{+48}_{-35}~{}\mathrm{MeV}\ ,\quad
g=2.5^{+1.2}_{-0.9}~{}\mathrm{GeV}~{}.$

The mass and coupling are reasonably consistent with experiment, in
particular, when considering the experimental errors and the systematic
uncertainties in the lattice results.

Figure 6:  Coupled channel $D\bar{D}$, $D_{s}\bar{D}_{s}$ scattering in
partial wave $l\\!=\\!0$. The quantity
$\sqrt{\rho_{i}\rho_{j}}|t_{ij}|^{2}\propto p|\sigma|$ related to the number
of events in experiment is shown ($\rho=2p/E_{cm}$ and $t$ is the scattering
matrix).

  * •

Narrow $\bm{0^{++}}$ state dominated by $\mathrm{D_{s}\bar{D}_{s}}$ slightly
below threshold and relation to $\bm{X(3915)}/\bm{\chi_{c0}(3930)}$

The coupled $D\bar{D}-D_{s}\bar{D}_{s}$ scattering gives also an indication
for a narrow $0^{++}$ resonance slightly below the $D_{s}\bar{D}_{s}$
threshold, which shows up as the sharp increase at the threshold in Fig. 6. It
has a large coupling to $D_{s}\bar{D}_{s}$ and a very small coupling to
$D\bar{D}$, as shown in Fig. 6 of [2]. The latter is responsible for its small
decay rate to $D\bar{D}$ and the small total width. We compare the resulting
resonance parameters with the experimental state $X(3915)/\chi_{c0}(3930)$
that shares similar features to the state we find

$\displaystyle\mathrm{lat}:\ \ $ $\displaystyle
m-2m_{D_{s}}=-0.2^{~{}+0.16}_{~{}-4.9}~{}\mathrm{MeV}\ ,\quad\
g=0.10^{~{}+0.21}_{~{}-0.03}~{}\mathrm{GeV}~{}$ (4)
$\displaystyle\mathrm{exp}\ X(3915)/\chi_{c0}(3930):\ $ $\displaystyle\
m-2m_{D_{s}}=-14.3\pm 1.8~{}\mathrm{MeV}~{},\enskip g=0.69\pm
0.07~{}\mathrm{GeV}~{}.$

The $\chi_{c0}(3930)$ with $J^{PC}=0^{++}$ was recently discovered in
$D\bar{D}$ decay by LHCb [6], while the $X(3915)$ with $J^{PC}\\!=\\!0^{++}$
or $2^{++}$ was found in $J/\psi\omega$ decay by Belle [13] and BaBar [14, 15,
16]: both states are now listed as the same state in the PDG [3]. It lies high
above the $D\bar{D}$ threshold, so one would naturally expect a broad width if
it was a conventional charmonium, given the large phase space available to
$D\bar{D}$ decay. Its narrow experimental width $\Gamma=18.8\pm 3.5~{}$MeV can
only be explained if its decay to $D\bar{D}$ is suppressed by some mechanism.

Our lattice results suggest that the $X(3915)/\chi_{c0}(3930)$ charmonium-like
state is likely not a conventional charmonium $\bar{c}c$, but owes its
existence to a large interaction in the $D_{s}\bar{D}_{s}$ channel near
threshold. It lies just below this threshold on the lattice and in experiment.
Our results indicate that it has a very small coupling to $D\bar{D}$, which
explains why its width is small and its decay to $D\bar{D}$ is suppressed in
experiment.

A state dominated by $\bar{c}c\bar{s}s$ in this energy region was found also
in phenomenological studies of molecular [17] and diquark-antidiquark [18, 19,
20] states.

  * •

$\bm{1^{--}}$ state related to $\bm{\psi(3770)}$

The $\psi(3770)$ related to conventional $\bar{c}c$ state $1^{3}D_{1}$ appears
in experiment just above the $D\bar{D}$ threshold. In our lattice simulation
with $m_{u/d}>m_{u/d}^{exp}$ and $m_{c}\gtrsim m_{c}^{exp}$, we find it as a
bound state in the $D\bar{D}$ scattering amplitude for partial wave $l=1$
slightly below threshold. The resulting resonance parameters show nice
agreement with experiment333The coupling $g$ is defined in (1), while the
value of $\sqrt{6\pi}g$ is more commonly listed for vector resonances. The
coupling for the bound state was extracted using
$p^{3}\cot\delta_{1}/E_{cm}=(m^{2}-E_{cm}^{2})/g^{2}$. :

$\displaystyle\mathrm{lat}:\ \ $ $\displaystyle m-M_{av}=707\pm
7~{}\mathrm{MeV}\ ,\quad\qquad\sqrt{6\pi}g=18.9^{~{}+0.8}_{~{}-0.7}~{}$ (5)
$\displaystyle\psi(3770)\ \ \mathrm{exp}:\ $ $\displaystyle\
m-M_{av}=704.25\pm 0.35~{}\mathrm{MeV}~{},\enskip\sqrt{6\pi}g=18.7\pm 0.9~{}.$

We note that $\psi(3770)$ is found as a resonance for our lighter charm quark
mass in [1] and the corresponding coupling $g$ has a very similar value.

  * •

$\bm{3^{--}}$ resonance related to $\bm{\psi_{3}(3842)}$

The $D\bar{D}$ scattering with partial wave $l=3$ shows a resonance pole that
corresponds to the charmonium with $J^{PC}=3^{--}$:

$\displaystyle\mathrm{lat}:\ \ $ $\displaystyle
m-M_{av}=754~{}{{}^{+4}_{-7}}~{}\mathrm{MeV}$ (6)
$\displaystyle\psi_{3}(3842)\ \ \mathrm{exp}:\ $ $\displaystyle\
m-M_{av}=773.9\pm 0.2~{},\enskip\Gamma=2.8\pm 0.6~{}.$

We are not able to determine its width since there is no $D\bar{D}$ lattice
energy level within the energy region of this narrow resonance. The candidate
for conventional charmonium with $J^{PC}=3^{--}$ was discovered recently by
LHCb [21] and the mass is close to our lattice value.

##  4 Conclusions

Charmonium-like states with isospin zero are studied by simulating one-channel
and coupled-channel scattering on the lattice. We find all conventional
charmonium resonances and bound states with
$J^{PC}=0^{++},~{}1^{--},~{}2^{++},~{}3^{--}$ up to the $D_{s}\bar{D}_{s}$
threshold. In addition, the results suggest the existence of unconventional
states just below $D_{s}\bar{D}_{s}$ and $D\bar{D}$ thresholds, where the
first is likely related to $X(3915)/\chi_{c0}(3930)$, while the second state
has not been discovered yet. The study makes several simplifying assumptions
necessary for a first investigation of this coupled-channel system.

Acknowledgments

We thank G. Bali, V. Baru, T. Gershon, F.-K. Guo, B. Hörz, D. Johnson, C. B.
Lang, R. Molina, J. Nieves, E. Oset, S. Paul, A. Schäfer and J. Simeth for
useful discussions. The simulations were performed on the Regensburg iDataCool
and Athene2 clusters, and the SFB/TRR 55 QPACE 2 and QPACE 3 machines. The
work is supported by DFG SFB/TRR-55, the European Union’s Horizon 2020
Research and Innovation programme under the Marie Sklodowska-Curie grant
813942 (ITN EuroPLEx), the STRONG-2020 project grant 824093, EU grant MSCA-IF-
EF-ST-744659 (XQCDBaryons) and Slovenian Research Agency grants P1-0035 and
J1-8137.

## References

  * [1] S. Piemonte, S. Collins, D. Mohler, M. Padmanath and S. Prelovsek,  Phys. Rev. D 100, 074505 (2019), [arXiv:1905.03506].
  * [2] S. Prelovsek, S. Collins, D. Mohler, M. Padmanath and S. Piemonte,  JHEP 06, 035 (2021), [arXiv:2011.02542].
  * [3] Particle Data Group, P. Zyla et al.,  PTEP 2020, 083C01 (2020).
  * [4] M. Bruno et al.,  JHEP 02, 043 (2015), [arXiv:1411.3982].
  * [5] C. Morningstar et al.,  Nucl. Phys. B 924, 477 (2017), [arXiv:1707.05817].
  * [6] LHCb, R. Aaij et al.,  2009.00026.
  * [7] D. Gamermann, E. Oset, D. Strottman and M. Vicente Vacas,  Phys. Rev. D 76, 074016 (2007), [arXiv:hep-ph/0612179].
  * [8] C. Hidalgo-Duque, J. Nieves, A. Ozpineci and V. Zamiralov,  Phys. Lett. B 727, 432 (2013), [arXiv:1305.4487].
  * [9] V. Baru et al.,  Phys. Lett. B 763, 20 (2016), [arXiv:1605.09649].
  * [10] C. Lang, L. Leskovec, D. Mohler and S. Prelovsek,  JHEP 09, 089 (2015), [arXiv:1503.05363].
  * [11] Belle, S. Uehara et al.,  Phys. Rev. Lett. 96, 082003 (2006), [arXiv:hep-ex/0512035].
  * [12] Belle, K. Chilikin et al.,  Phys. Rev. D 95, 112003 (2017), [arXiv:1704.01872].
  * [13] Belle, S. Uehara et al.,  Phys. Rev. Lett. 104, 092001 (2010), [arXiv:0912.4451].
  * [14] BaBar, B. Aubert et al.,  Phys. Rev. Lett. 101, 082001 (2008), [arXiv:0711.2047].
  * [15] BaBar, P. del Amo Sanchez et al.,  Phys. Rev. D 82, 011101 (2010), [arXiv:1005.5190].
  * [16] BaBar, J. Lees et al.,  Phys. Rev. D 86, 072002 (2012), [arXiv:1207.2651].
  * [17] X.-K. Dong, F.-K. Guo and B.-S. Zou,  Progr. Phys. 41, 65 (2021), [arXiv:2101.01021].
  * [18] R. F. Lebed and A. D. Polosa,  Phys. Rev. D 93, 094024 (2016), [arXiv:1602.08421].
  * [19] J. F. Giron and R. F. Lebed,  Phys. Rev. D 102, 014036 (2020), [arXiv:2005.07100].
  * [20] W. Chen, H.-X. Chen, X. Liu, T. Steele and S.-L. Zhu,  Phys. Rev. D 96, 114017 (2017), [arXiv:1706.09731].
  * [21] LHCb, R. Aaij et al.,  JHEP 07, 035 (2019), [arXiv:1903.12240].