“Rosenbluth” separation of the $J/\psi$ near-threshold photoproduction - an
access to the gluon Gravitational Form Factors at high $t$

###### Abstract

We perform analysis of the near-threshold $J/\psi$ photoproduction data off
the proton based on the high skewness expansion ( $\xi\to 1$) of the cross-
section in Ref. Guo  _et al._ (2024). This expansion allows to separate the
corresponding gluonic form factors, in much the same way as this is done in
the electro-magnetic case with the Rosenbluth separation. We examine the
independence of the extracted form factors with the photon beam energy. In
Ref. Guo  _et al._ (2024) these form factors, in the leading-moment
approximation, are related to the proton’s gluon Gravitational Form Factors
(gGFF). Using this approximation, we compare the extracted data points with
lattice calculations in the region of momentum transfer squared of $1<|t|<2$
GeV2, where the data and lattice results overlap. Such analysis demonstrates
the possibility of extracting some combinations of the gGFFs from the data at
high $t$, complementary to the lattice calculations available in the low $t$
region. However, higher statistics is needed to more accurately check the
predicted $\xi$-scaling behavior of the data and compare with the lattice
results, thus testing the theoretical assumptions used in Ref. Guo  _et al._
(2024).

##  I Introduction

The charmonium photoproduction near threshold provides a unique way to probe
the gluonic structure of the nucleon. Due to the high mass of the charm quark
and proximity to the threshold, the reaction is expected to be dominated by
two-gluon exchange Brodsky and Miller (1997); Sun  _et al._ (2021). In
addition, assuming factorization where the hard scale is defined by the high
charmonium mass, the process involves the gluonic Generalized Parton
Distribution (GPD) of the nucleon Guo  _et al._ (2021, 2023). As the two
gluons can mimic graviton-like exchange, the gluonic GPD in turn can be
related to the gluon Gravitational Form Factors (gGFFs) Guo  _et al._ (2021);
Hatta and Strikman (2021). To be able to access the gluonic content of the
nucleon, the above chain of assumptions requires experimental testing.

In the recent years, thanks to the 12 GeV Jefferson Lab electron accelerator,
several experiments measured the $J/\psi$ photoproduction near threshold with
much higher precision than the previous experiments from the 1970s. The GlueX
experiment in Hall D performed measurements with nearly full acceptance
spectrometer in a wide range of the photon beam energy starting from the
threshold of $8.2$ GeV to the maximum available energy of about $11.8$ GeV Ali
et al. (2019) (GlueX collaboration); Adhikari et al. (2023) (GlueX
collaboration). In addition to the total cross-section, the differential cross
sections were extracted in three slices of the photon beam energy, providing
full coverage of the near-threshold kinematic region. The $J/\psi$-007
experiment in Hall C accumulated similar statistics in a narrower kinematic
region with several settings of their two low-acceptance spectrometers that
however allowed very high luminosity Duran et al. (2023) (J/$\psi$-007
collaboration). The differential cross sections were measured in 10 narrow
energy bins with very high statistics in the forward region.

In Ref. Duran et al. (2023) (J/$\psi$-007 collaboration) the gGFFs have been
extracted by fitting the $J/\psi$-007 data with the theoretical predictions
from the holographic Mamo and Zahed (2021) and GPD Guo  _et al._ (2021)
approaches. Note that the fits used the full $J/\psi$-007 data set dominated
by the cross sections at low $|t|$. In the GPD analysis of Ref. Guo  _et al._
(2023), the data from both GlueX and $J/\psi$-007 experiments were used to
extract the gGFFs. The GPD approach is based on expansion in the skewness
parameter111The skewness is one of the three GPD variables (in addition to $t$
and the parton longitudinal momentum fraction $x$) defined as the longitudinal
momentum transfer fraction - see Eq.(7) in Ref. Guo  _et al._ (2021) for
$\xi\to 1$ (see also Hatta and Strikman (2021)) therefore its validity is
justified for high $|t|$ values (due to kinematic correlation between $t$ and
$\xi$). A common feature of the above two analyses in Refs. Duran et al.
(2023) (J/$\psi$-007 collaboration); Guo  _et al._ (2023) is the modeling of
the gGFFs with dipole/tripole functions, as proposed in Ref. Frankfurt and
Strikman (2002) in analogy with the electro-magnetic form factors. The choice
of this functional form comes also from the parametrization of the lattice
results in Refs. Pefkou  _et al._ (2022); Hackett  _et al._ (2023) so that
some parameters in these parametrizations can be fixed from the lattice
calculations when fitting the data. While the GPD approach works at high
$|t|$, the lattice calculations of the gGFFs have been done so far in the
$0<|t|<2$ GeV2 region. Therefore, using the lattice parametrization of the
gGFFs and assuming the same functional form to extend to high $|t|$ may result
in significant bias.

The analysis in this work is based on one of the important results of Ref. Guo
_et al._ (2024), that the matrix element squared of the charmonium
photoproduction for high $\xi$ values can be presented in this general form:

$\displaystyle|G(\xi,t)|^{2}$
$\displaystyle=\xi^{-4}\left[G_{0}(t)+\xi^{2}G_{2}(t)+\xi^{4}G_{4}(t)\right],$
(1)

where $G_{0}$, $G_{2}$, and $G_{4}$ are functions only of $t$, they are
related to the gluonic Compton-like Form Factors (gCFFs). We aim to extract
the individual terms in Eq.(1) from the experimental data with some minimal
additional assumptions. In contrast to the studies in Refs. Duran et al.
(2023) (J/$\psi$-007 collaboration); Guo  _et al._ (2023), the results of the
procedure in this work will not be model functions fitted to the data, but
actual data points measuring the $t$-dependence of the above form factors.

We point out the analogy with the Rosenbluth separation of the electric,
$G_{E}$, and magnetic, $G_{M}$, form factors. Based on the Rosenbluth formula
Rosenbluth (1950) for the electron-proton elastic cross section,
$d\sigma/d\Omega$, the method relies on the linear dependence of the reduced
cross section, $\sigma_{R}$, on the kinematic variable $\epsilon$:

$\sigma_{R}=\frac{d\sigma}{d\Omega}/\left(\frac{d\sigma}{d\Omega}\right)_{M}\frac{\epsilon(1+\tau)}{\tau}=\frac{\epsilon}{\tau}G_{E}^{2}(t)+G_{M}^{2}(t),$
(2)

where $\left(\frac{d\sigma}{d\Omega}\right)_{M}$ is the Mott cross section,
$\tau=-t/4m^{2}$, where $m$ is the proton mass. By measuring $\sigma_{R}$ at
the same momentum transfer $t$ as function of $\epsilon/\tau$, one can
estimate the slope and the intercept, being respectively the electric and
magnetic form factors squared at this value of $t$.

The analogy of Eq.(1) with Eq.(2) is clear: instead of $\epsilon$, we can use
the variations with $\xi$ at a fixed $t$ to extract the form factors in
Eq.(1). In this paper we will apply a similar “Rosenbluth” technique, however
with some noticeable differences. Due to the limited experimental data so far,
it is practically impossible to analyze separately data sets in narrow ranges
of $t$ versus some function of $\xi$. At the same time global analysis of the
data is complicated due to kinematic constraints that impose correlations
between the ranges of $t$ and $\xi$.

In our analysis we will use the results from the two Jefferson Lab
experiments, GlueX and $J/\psi$-007, selecting data points with $\xi>0.4$.
Such limit may not be high enough for the assumptions in Ref. Guo  _et al._
(2024) to be valid, and the selection is driven by the limited amount of data
at high $\xi$. In addition we will select measurements with photon energies
above $9.3$ GeV, to stay away from the open-charm thresholds of
$\Lambda_{c}\bar{D}$ and $\Lambda_{c}\bar{D}^{*}$. As discussed recently in
Refs.Du  _et al._ (2020); Adhikari et al. (2023) (GlueX collaboration); Winney
_et al._ (2023), there is possible evidence for contributions from reactions
with open-charm exchange.

In section II we will formulate the “Rosenbluth” technique. It will be used in
section III to extract the $G_{0}(t)$ and $G_{2}(t)$ form factors and check
the energy independence of the results. In section IV we will work within the
leading-moment approximation to demonstrate the possibility of extracting some
combinations of the gGFFs and compare them to the lattice results. Some
additional remarks and summary will be given in section V.

##  II The “Rosenbluth” technique

In Ref. Guo  _et al._ (2024) the differential cross section of the $J/\psi$
photoproduction near threshold is presented as a function of the skewness
parameter $\xi$ and the momentum transfer $t$:

$\frac{d\sigma}{dt}(\xi,t)=\frac{\alpha
e_{q}^{2}}{4(W^{2}-m^{2})^{2}}\frac{(16\pi\alpha_{S})^{2}}{3M^{3}}|\Psi_{NR}(0)|^{2}\cdot|G(\xi,t)|^{2}=F\cdot|G(\xi,t)|^{2},$
(3)

where the factor $F$ includes a kinematic factor and the $J/\psi$ non-
relativistic wave function at origin, $\Psi_{NR}(0)$, and depends only on the
photon beam energy, $E$. The matrix element squared, $|G(\xi,t)|^{2}$, given
in Eq.(1), has three terms related to the gluonic Compton Form Factors (gCFF),
${\cal{A}}_{g}$, ${\cal{B}}_{g}$, and ${\cal{C}}_{g}$, as:

$\displaystyle G_{0}(t)$
$\displaystyle=\left({\cal{A}}_{g}(t)\right)^{2}-\frac{t}{4m^{2}}\left({\cal{B}}_{g}(t)\right)^{2}$
(4) $\displaystyle G_{2}(t)$
$\displaystyle=2{\cal{A}}_{g}(t){\cal{C}}_{g}(t)+2\frac{t}{4m^{2}}{\cal{B}}_{g}(t){\cal{C}}_{g}(t)-\left({\cal{A}}_{g}(t)+{\cal{B}}_{g}(t)\right)^{2}$
(5) $\displaystyle G_{4}(t)$
$\displaystyle=\left(1-\frac{t}{4m^{2}}\right)\left({\cal{C}}_{g}(t)\right)^{2}.$
(6)

From Eqs.(3) and (1) we construct the reduced cross section:

$\displaystyle\sigma_{R0}(E,t)=\frac{d\sigma}{dt}(E,t)\frac{\xi^{2}(E,t)}{F(E)}=\xi^{-2}(E,t)G_{0}(t)+G_{2}(t)+\xi^{2}(E,t)G_{4}(t),$
(7)

where we have indicated the energy and $t$ dependence of the corresponding
variables222Note that all reduced cross sections defined in this work are
dimensionless quantities. We have verified, initially based on the lattice
results from Refs. Pefkou  _et al._ (2022); Hackett  _et al._ (2023) discussed
later, that the last term in Eq.(1), and respectively Eq.(7), is much smaller
in absolute value than each of the other two terms. Partially this is because
of the $\xi^{4}$/$\xi^{2}$ suppression factor (as $\xi<1$) with respect to the
first/second term. This assumption will be checked later based on the results
of this analysis. Note that this does not necessarily mean that the last term
in Eqs.(1) and (7) can be neglected. It will be shown that the other two terms
have different signs and it is possible they largely cancel each other.

The second term in Eq.(7) depends only on $t$ and we can eliminate it by
taking the difference of the cross sections at two different photon energies,
$E_{i}$ and $E_{j}$:

$\displaystyle\sigma_{R0}(E_{i},t)-\sigma_{R0}(E_{j},t)$
$\displaystyle=\left[\xi^{-2}(E_{i},t)-\xi^{-2}(E_{j},t)\right]G_{0}(t)$ (8)
$\displaystyle+\left[\xi^{2}(E_{i},t)-\xi^{2}(E_{j},t)\right]G_{4}(t).$

Now the last term can be neglected with respect to the first one, and we have
the $G_{0}$ function extracted from experimental data:

$\displaystyle
G_{0}^{exp.}(t)=\left[\sigma_{R0}(E_{i},t)-\sigma_{R0}(E_{j},t)\right]/\left[\xi^{-2}(E_{i},t)-\xi^{-2}(E_{j},t)\right]$
(9)

We will choose a data set at a given energy $E_{i}$ as a reference. Thus, the
above equation will give us the slope of the reduced cross section as function
of $\xi^{-2}$ for each additional data point with transfer momentum $t$ and
energy $E_{j}$ with respect to the reference energy, for which the slope is
$G_{0}(t)$. With this method the uncertainties of the reference data set can
be considered to be systematic errors common for all data points.

In much the same way we can extract $G_{2}(t)$ in Eq.(1) as:

$\displaystyle
G_{2}^{exp.}(t)=\left[\sigma_{R2}(E_{i},t)-\sigma_{R2}(E_{j},t)\right]/\left[\xi^{2}(E_{i},t)-\xi^{2}(E_{j},t)\right],$
(10)

where we have defined another reduced cross section as:

$\displaystyle\sigma_{R2}(E,t)=\frac{d\sigma}{dt}(E,t)\frac{\xi^{4}(E,t)}{F(E)}$
(11)

##  III Energy independence of the form factors

Figure 1:  Data points of the GlueX Adhikari et al. (2023) (GlueX
collaboration) (red solid) and $J/\psi$-007 Duran et al. (2023) (J/$\psi$-007
collaboration) (black open) experiments on the $E$-vs-$t$ plane. The vertical
error bars represent the relative errors of the differential cross sections
measured at these points in arbitrary units (not related to the $y$ axis). The
$\xi$ values on this plot are given on the $z$-axis (color coded). Some of the
$J/\psi$-007 points are slightly shifted up in energy for visibility. The
colored (shaded in print) area corresponds to the $\xi>0.4$ and $E>9.3$ GeV
region that is used in this analysis.

The practical implementation of the above technique is limited by the quality
and range of the available data. The data sets from the two Jefferson Lab
experiments are illustrated in Fig.1. As discussed above for the main part of
our analysis we will use the data with $\xi>0.4$ and $E>9.3$ GeV. This region
on the $E$-vs-$t$ plane in Fig.1 is shown with color. We will choose the
highest energy data set from the GlueX experiment, at an average energy
$E_{i}=10.82$ GeV, as a reference and subtract the other data at different
energies $E_{j}$. Such a choice gives maximum leverage in $\xi$ when
calculating the slopes from Eqs.(9) and (10), at the same time the highest
energy data set is more precise and has the widest range in $t$.

Figure 2:  Global four-parameter fit of the reduced cross sections,
$\xi^{2}/F\cdot d\sigma/dt$, as function of $|t|$ (in GeV2), using the model
of Ref.Guo  _et al._ (2024), from GlueX (red) and $J/\psi$-007 (black)
experiments for different photon energies as indicated. For the further
analysis only the fitted function at $10.82$ GeV is used.

In principle, we can fit the reference reduced cross section with some
function without any model, so that we can use it in Eqs.(9) and (10) to
subtract the reduced cross sections at the $t$ values of the data points at
the other energies. Still, as seen in Fig.1, due to the $\xi>0.4$ constraint,
the $1-2$ GeV2 region in $|t|$ is not covered by the reference cross section,
therefore it requires extrapolation to lower $|t|$. For that we will use the
results from Ref. Guo  _et al._ (2024) in the leading-moment approximation to
fit all the data simultaneously. An alternative model-independent approach
will be discussed later. In the leading-moment approximation the gCFFs
${\cal{A}}_{g}$, ${\cal{B}}_{g}$, and ${\cal{C}}_{g}$ in Eqs.(4-6) are related
to the gluonic Gravitational Form Factors (gGFF) $A_{g}$, $B_{g}$, and $C_{g}$
as:

$\displaystyle{\cal{A}}_{g}(t)$ $\displaystyle=2A_{1}^{conf}A_{g}(t)$ (12)
$\displaystyle{\cal{B}}_{g}(t)$ $\displaystyle=2A_{1}^{conf}B_{g}(t)$ (13)
$\displaystyle{\cal{C}}_{g}(t)$ $\displaystyle=8A_{1}^{conf}C_{g}(t),$ (14)

using the conformal expansion of the gCFFs where $A_{1}^{conf}=5/4$ Guo  _et
al._ (2024). For this fit we use a tripole parametrization for all three
gGFFs, $F_{g}(0)/(1-t/m_{g}^{2})^{3}$, each with two parameters, the form
factor at $t=0$, $F_{g}(0)$, and the mass scale, $m_{g}$. We also fix:

$\displaystyle A_{g}(0)$ $\displaystyle=0.414$ (15) $\displaystyle C_{g}(0)$
$\displaystyle=-0.642$ (16)

from the CT18 global analysis of the parton distribution functions Hou  _et
al._ (2021) for $A_{g}(0)$, and lattice results in Ref. Hackett  _et al._
(2023) for $C_{g}(0)$, keeping the two parameters for $B_{g}$ free. In
Refs.Duran et al. (2023) (J/$\psi$-007 collaboration); Guo  _et al._ (2023)
when extracting the gGFFs, $B_{g}$ is neglected based on the lattice results
for $|t|<2$ GeV2. Here we will keep it as it improves the fit quality at high
$|t|$. We demonstrate this procedure when fitting the reduced cross section,
$\sigma_{R0}$, as function of $t$ for all energies using Eq. (7), where the
$G$ functions are taken from Eqs.(4-6,12-14). The results of this global fit
are shown in Fig.2 exhibiting a good fit quality. In the further analysis for
the extraction of $G_{0}^{exp.}(t)$, in Eq.(9) we will use only the fitted
function for the reference energy $E_{i}=10.82$ GeV (Fig.2 bottom right). The
global fit of the other reduced cross section, $\sigma_{R2}$, gives very
similar results, from which fit only the fitted function for the reference
energy will be used in Eq.(10) to extract $G_{2}^{exp.}(t)$.

Figure 3:  Different stages of applying the “Rosenbluth” extraction of the
$G_{0}(t)$ function using the GlueX data: the differential cross sections at
three energy slices normalized by $F(E)$ \- see Eq.(3) (top), the reduced
cross sections - Eq.(7) (middle), and the extracted $G_{0}(t)$ function with
Eq.(9) using the highest energy as a reference (bottom). The fitted function
of the reduced cross section at the highest energy from the global fit (Fig.2
bottom right), is shown in the middle panel (blue curve), as well. The points
are connected with lines to guide the eye only.

(a)

(b)

Figure 4:  $G_{0}(t)$ (top) and $G_{2}(t)$ (bottom) functions as extracted
from the GlueX (solid) and $J/\psi$-007 (open symbols) data for $\xi>0.4$
using Eqs.(9,10) for different photon beam energies (different colors
connected with lines). The data are fitted with dipole functions squared,
$G(0)/(1-t/m_{G}^{2})^{4}$, shown with one-sigma error band.

To illustrate clearly how the “Rosenbluth” technique works, in Fig. 3 we show
the different stages of the extraction of the $G_{0}(t)$ function from the
GlueX data. Here, just for this demonstration, we have included the
differential cross section at the lowest beam energy. Starting from the
original cross sections, going to the $\xi$-scaled cross sections, and finally
calculating the slopes from Eq.(9), it is remarkable to see how the results at
the different energies converge into a function, $G_{0}^{exp.}(t)$, that is
largely energy independent.

The results for $G_{0}^{exp.}$ and $G_{2}^{exp.}$ using all the data are shown
in Figs. 4(a) and 4(b). We see that the extracted functions do not depend on
energy within the experimental errors. Numerically, this is confirmed by the
good quality of the fit of all the data using a dipole function. The
uncertainties of the fitted function at the reference energy (as in Fig. 2
bottom right panel) are scaled by the factors in Eqs.(9,10) and shown in
Figs.4(a) and 4(b) as global systematic errors.

While we applied the GPD model of Ref. Guo  _et al._ (2024) in the leading-
moment approximation to fit all the data, in the above procedure we have only
used the fitted function for the reference energy mainly to extrapolate the
reduced cross section from $|t|\sim 2$ GeV2, where the lowest-$|t|$ data point
is, down to $|t|\sim 1$ GeV2 \- see the curve in the middle panel in Fig.3. We
found that similar results can be obtained when fitting only the reference
reduced cross section (though with bigger fit uncertainties) with some model-
independent function if its value is fixed at a point $|t|\leq 1$ GeV2 by some
model constraints like Eqs.(15,16), which are important for the extrapolation
to the lower $|t|$ region. Nevertheless, when enough statistics is available,
it would be possible to chose a free parameter that constrains the fit of the
reference cross section at some $t$ value. We determine this parameter by
minimizing the $\chi^{2}$ of the fits of $G_{0}(t)$ and $G_{2}(t)$, i.e.
requiring $G_{0}(t)$ and $G_{2}(t)$ to be energy independent.

##  IV Relation to the gluon Gravitational Form Factors

In this section we will use the leading-moment approximation of Ref. Guo  _et
al._ (2024) that relates the $G$ form factors in Eq.(1) to the gluon
Gravitational Form Factors, as given by Eqs.(4-6) and Eqs.(12-14). In Fig.5(a)
and Fig.5(b), we compare our results from the previous section to the
calculations of $G_{0}$ and $G_{2}$ in the leading-moment approximation using
the lattice results for the gGFFs from Ref. Hackett  _et al._ (2023). We see a
good agreement between the lattice and experimental results in the common
$1-2$ GeV2 region, especially for $G_{0}$. However, this is partially a result
of the constraints from Eqs.(15,16), as discussed above. To demonstrate the
deviation from the $\xi$-scaling, in addition to the data with $\xi>0.4$ used
in the analysis, we show also the results for $\xi<0.4$. Notably, the
$\xi<0.4$ data points for both $G_{0}$ and $G_{2}$ exhibit larger spread and
deviate significantly from the lattice calculations.

(a)

(b)

Figure 5:  $G_{0}(t)$ (top) and $G_{2}(t)$ (bottom) functions as extracted
from the GlueX and $J/\psi$-007 data for $\xi>0.4$ and $\xi<0.4$ using
Eqs.(9,10). The data for $\xi>0.4$ GeV are fitted with dipole function
squared, $G(0)/(1-t/m_{G}^{2})^{4}$, shown with one-sigma error band. The
results are compared with the calculations of $G_{0}$ and $G_{2}$ using the
lattice results for the gGFFs in Ref. Hackett  _et al._ (2023) in leading-
moment approximation of Ref. Guo  _et al._ (2024).

(a)

(b)

Figure 6:  As in Figs.4,5 but for the $G_{0}+G_{2}$ sum, using Eq.(18). The
data are fitted with a dipole function squared,
$(G_{0}(0)+G_{2}(0))/(1-t/m_{G02}^{2})^{4}$.

The $A_{g}$ form factor squared is proportional to $G_{0}$ due to the small
$B_{g}^{2}$ contribution (see Eqs.(4,12)) and can be extracted from the
experimental data for $|t|>1$ GeV2. We can use the above results for $G_{0}$
and $G_{2}$ to estimate the $C_{g}$ form factor from Eqs.(4,5,12,13). However,
as the calculations of $G_{0}$ and $G_{2}$ use the same data, the best way to
do this is to extract directly the sum $G_{0}+G_{2}$ in which $A_{g}^{2}$
cancels out. For that we rearrange Eq.(7) as:

$\displaystyle\frac{d\sigma}{dt}(E,t)\frac{1}{F(E)}\frac{\xi^{4}(E,t)}{1-\xi^{2}(E,t)}$
$\displaystyle=G_{0}(t)+\frac{\xi^{2}(E,t)}{1-\xi^{2}(E,t)}(G_{0}(t)+G_{2}(t))$
$\displaystyle+\frac{\xi^{4}(E,t)}{1-\xi^{2}(E,t)}G_{4}(t).$

In much the same way as for $G_{0}^{exp.}(t)$ and $G_{2}^{exp.}(t)$, we can
extract $G_{0}+G_{2}$ as:

$\displaystyle(G_{0}(t)+G_{2}(t))^{exp.}=\left[\sigma_{R02}(E_{i},t)-\sigma_{R02}(E_{j},t)\right]/\left[\frac{\xi^{2}(E_{i},t)}{1-\xi^{2}(E_{i},t)}-\frac{\xi^{2}(E_{j},t)}{1-\xi^{2}(E_{j},t)}\right],$
(18)

where we have defined another reduced cross section as:

$\displaystyle\sigma_{R02}(E,t)=\frac{d\sigma}{dt}(E,t)\frac{\xi^{4}(E,t)}{F(E)(1-\xi^{2}(E,t))}.$
(19)

The results for the sum $(G_{0}+G_{2})^{exp.}$ are shown in Fig.6, in the same
way as this was done in Figs.4,5 for $G_{0}^{exp.}$ and $G_{2}^{exp.}$.

Figure 7:  The $A_{g}$ form factor squared (top), and the product of the two
form factors with negative sign $-A_{g}C_{g}$ (bottom), extracted from the
experimental data within the leading-moment approximation of Ref. Guo  _et
al._ (2024), fitted with dipole functions (one-sigma error band shown) and
compared to the lattice results of Ref. Hackett  _et al._ (2023).

In the leading-moment approximation, when neglecting the $B_{g}$ form factor,
we have:

$\displaystyle G_{0}(t)=(2A_{1}^{conf})^{2}A_{g}^{2}(t)$ (20) $\displaystyle
G_{0}(t)+G_{2}(t)=2(2A_{1}^{conf}A_{g}(t))(8A_{1}^{conf}C_{g}(t)).$ (21)

Thus, from the data we can extract directly $A_{g}^{2}(t)$ and the product of
the two form factors, $A_{g}(t)C_{g}(t)$, being proportional to $G_{0}$ and
$G_{0}+G_{2}$, respectively. The results are shown in Fig.7.

Note that the mass scale parameter for $G_{0}$ is bigger than the one for
$G_{0}+G_{2}$, $m_{G0}>m_{G02}$. This means that $C_{g}$ vanishes faster with
$|t|$ than $A_{g}$. Thus, in addition to the $\xi$ suppression at low $|t|$,
this explains qualitatively the smaller contribution of $G_{4}\sim C_{g}^{2}$
in Eq.(1) at high $|t|$ values. Quantitatively, we have verified based on the
results in this paper, that for $|t|>1$ GeV2, $G_{4}$ is much smaller in
absolute value than $\xi^{-4}G_{0}$ and $\xi^{-2}G_{2}$. Thus, in the leading-
moment approximation, this justifies neglecting the $G_{4}$ term when
extracting the form factors.

##  V Additional remarks and summary

An interesting result can be deduced from the general $\xi$-scaling
expression, Eq.(1), though not directly related to the extraction of the form
factors. First, we note that $G_{0}$ must be positive as a sum of squares of
gCFFs, according to Eq.(4) with $t$ being negative. Therefore, from Eq.(9),
for $E_{i}>E_{j}$ at a fixed $t$, as $\xi(E_{i},t)<\xi(E_{j},t)$ and therefore
$\xi^{-2}(E_{i},t)-\xi^{-2}(E_{j},t)>0$, it follows that:

$\displaystyle
d\sigma/dt(E_{i},t)\frac{\xi^{2}(E_{i},t)}{F(E_{i})}>d\sigma/dt(E_{j},t)\frac{\xi^{2}(E_{j},t)}{F(E_{j})},\;E_{i}>E_{j}.$
(22)

In particular, as $\xi^{2}(E_{i},t)/F(E_{i})<\xi^{2}(E_{j},t)/F(E_{j})$, it
means that:

$\displaystyle d\sigma/dt(E_{i},t)>d\sigma/dt(E_{j},t),\;E_{i}>E_{j}.$ (23)

In other words, the differential cross section at a fixed $t$ must increase
with the energy, something that is not trivial. Based on Eq.(22) and using the
same arguments in reverse we can claim from Eq.(10) that $G_{2}(t)<0$, as
$\xi^{2}(E_{i},t)-\xi^{2}(E_{j},t)<0$. This statement however is not strict
mathematically as, compared to Eq.(22), we need to take into account the
additional factor of $\xi^{2}$ in the reduced cross section, Eq.(11), though
this effect turned out to be relatively small. In the leading-moment
approximation, $G_{2}(t)<0$ would mean that $A_{g}$ and $C_{g}$ have opposite
signs.

To summarize, our “Rosenbluth” procedure relies on Eq.(1) that comes in Ref.
Guo  _et al._ (2024) from some general theoretical assumptions about the
$J/\psi$ photoproduction mechanisms, such as gluon exchange and factorization.
The observation of the energy independence of the extracted $G_{0}$ and
$G_{2}$ functions means that the experimental data, within their
uncertainties, are consistent with the assumptions in Ref. Guo  _et al._
(2024).

We conclude, based on the analysis of the current experimental data (Fig.5),
that the $\xi$-scaling behavior is valid for $\xi>0.4$. As discussed in Ref.
Guo  _et al._ (2024), such a low limit is not justified theoretically and may
result in a significant bias. However, this value should be considered as a
lower limit of the validity of the $\xi$-scaling, as more precise data could
move this value up.

We have compared the extracted $G_{0}$ and $G_{2}$ functions with the
calculations in the leading-moment approximation when using lattice gGFFs.
Note that our procedure relies partially on the lattice constraints of the
form factors when extrapolating the reference cross section to the lower $|t|$
region. Therefore, the data points in Figs.5,7 cannot be considered
independent when compared to the lattice results. In addition, the $\xi>0.4$
limit discussed above was also partially based on comparison with the lattice
calculations. A minimization procedure was proposed that would make it
possible to avoid the use of the lattice constraints, which however, requires
higher statistics. More data are needed to study the region of the validity of
the $\xi$-scaling. If with higher statistics the agreement between the data
and lattice calculations in the overlapping $t$-region is confirmed, this
could justify the extraction of the gGFFs at high $t$ directly from the data,
complementary to the lattice results available at low $t$.

It is critical to further investigate the near-threshold region of $E<9.3$
GeV. For that we would need more statistics and further theoretical and
phenomenological studies to understand the possible contribution of the open
charm channels. If gluon exchange still dominates at least in some kinematic
region, it will be very important to have this region included in the
presented analysis. It is believed that close to threshold, being close to
both $|t|_{min}$ and $|t|_{max}$ kinematic limits, the factorization would be
better justified, at the same time the higher $\xi$ values in this region
would better satisfy the requirements for the high-$\xi$ expansion in Ref.Guo
_et al._ (2024).

The results in this paper should not be considered as a proof of the validity
of the assumptions in Ref. Guo  _et al._ (2024), rather just as an indication
within the present data uncertainties, of the validity of the $\xi$-scaling as
defined by Eq.(1). Nevertheless, the presented results are encouraging and
should stimulate the continuation of the experimental studies and further
theoretical analysis.

##  VI ACKNOWLEDGMENTS

We would like to thank Yuxun Guo for the fruitful discussions and explanations
of the theory on which the results of this work are based and Simon Taylor for
valuable comments.

Notice: Authored by Jefferson Science Associates, LLC under U.S. DOE Contract
No. DE-AC05-06OR23177. The U.S. Government retains a non-exclusive, paid-up,
irre- vocable, world-wide license to publish or reproduce this manuscript for
U.S. Government purposes.

## References

  * Guo  _et al._ (2024) Y. Guo, X. Ji, and F. Yuan, Phys. Rev. D 109, 014014 (2024).
  * Brodsky and Miller (1997) S. J. Brodsky and G. A. Miller, Physics Letters B 412, 125 (1997).
  * Sun  _et al._ (2021) P. Sun, X.-B. Tong, and F. Yuan, Physics Letters B 822, 136655 (2021).
  * Guo  _et al._ (2021) Y. Guo, X. Ji, and Y. Liu, Physical Review D 103 (2021), 10.1103/physrevd.103.096010.
  * Guo  _et al._ (2023) Y. Guo, X. Ji, Y. Liu, and J. Yang, Physical Review D 108 (2023), 10.1103/physrevd.108.034003.
  * Hatta and Strikman (2021) Y. Hatta and M. Strikman, Phys. Lett. B 817, 136295 (2021).
  * Ali et al. (2019) (GlueX collaboration) A. Ali et al. (GlueX collaboration), Phys. Rev. Lett. 123, 072001 (2019).
  * Adhikari et al. (2023) (GlueX collaboration) S. Adhikari et al. (GlueX collaboration), Physical Review C 108 (2023), 10.1103/physrevc.108.025201.
  * Duran et al. (2023) (J/$\psi$-007 collaboration) B. Duran et al. (J/$\psi$-007 collaboration), Nature 615, 813 (2023).
  * Mamo and Zahed (2021) K. A. Mamo and I. Zahed, Phys. Rev. D 103 (2021), 10.1103/physrevd.103.094010.
  * Frankfurt and Strikman (2002) L. Frankfurt and M. Strikman, Phys. Rev. D 66, 031502 (2002).
  * Pefkou  _et al._ (2022) D. A. Pefkou, D. C. Hackett, and P. E. Shanahan, Physical Review D 105 (2022), 10.1103/physrevd.105.054509.
  * Hackett  _et al._ (2023) D. C. Hackett, D. A. Pefkou, and P. E. Shanahan, “Gravitational form factors of the proton from lattice qcd,” (2023), arXiv:2310.08484 [hep-lat] .
  * Rosenbluth (1950) M. N. Rosenbluth, Phys. Rev. 79, 615 (1950).
  * Du  _et al._ (2020) M.-L. Du, V. Baru, F.-K. Guo, C. Hanhart, U.-G. Meißner, A. Nefediev, and I. Strakovsky, Eur.Phys.J.C 80 (2020), 10.1140/epjc/s10052-020-08620-5.
  * Winney  _et al._ (2023) D. Winney, C. Fernández-Ramírez, A. Pilloni, A. N. Hiller Blin, M. Albaladejo, L. Bibrzycki, N. Hammoud, J. Liao, V. Mathieu, G. Montaña, R. J. Perry, V. Shastry, W. A. Smith, and A. P. Szczepaniak (Joint Physics Analysis Center), Phys. Rev. D 108, 054018 (2023).
  * Hou  _et al._ (2021) T.-J. Hou, J. Gao, T. J. Hobbs, K. Xie, S. Dulat, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, I. Sitiwaldi, D. Stump, and C.-P. Yuan, Phys. Rev. D 103, 014013 (2021).


∎

11institutetext:  Karl Jansen 22institutetext: Julia Volmer - Speakers,
33institutetext: DESY Zeuthen, Platanenallee 6, 15738 Zeuthen, Germany
33email<EMAIL_ADDRESS>33email<EMAIL_ADDRESS>44institutetext:
Tobias Hartung 55institutetext: Department of Mathematics, Kings College
London, Strand, London WC2R 2LS, United Kingdom
55email<EMAIL_ADDRESS>66institutetext: Hernan Leövey
77institutetext: Structured Energy Management, Axpo Trading, Parkstrasse 23,
5400 Baden, Germany
77email<EMAIL_ADDRESS>
# Avoiding the sign-problem in lattice field theory

Tobias Hartung  Karl Jansen  Hernan Leövey  Julia Volmer

###### Abstract

In lattice field theory, the interactions of elementary particles can be
computed via high-dimensional integrals. Markov-chain Monte Carlo (MCMC)
methods based on importance sampling are normally efficient to solve most of
these integrals. But these methods give large errors for oscillatory
integrands, exhibiting the so-called sign-problem. We developed new quadrature
rules using the symmetry of the considered systems to avoid the sign-problem
in physical one-dimensional models for the resulting high-dimensional
integrals. This article gives a short introduction to integrals used in
lattice QCD where the interactions of gluon and quark elementary particles are
investigated, explains the alternative integration methods we developed and
shows results of applying them to models with one physical dimension. The new
quadrature rules avoid the sign-problem and can therefore be used to perform
simulations at until now not reachable regions in parameter space, where the
MCMC errors are too big for affordable sample sizes. However, it is still a
challenge to develop these techniques further for applications with physical
higher-dimensional systems.

##  1 Introduction

Monte Carlo (MC) methods are in general very efficient to solve high-
dimensional integrals. They use the law of large numbers to approximate an
integral with quadrature rules that use random sampling points. But MC methods
are highly _inefficient_ for oscillatory integrand functions, e.g. the
function shown in Figure 1a. An exact integration of oscillatory functions
would, of course, result in the cancellation of large negative and positive
contributions to the integral - in the example in Figure 1a this would give an
integral of zero. However, the random choice of sampling points in MC methods,
shown as black points in Figure 1a, does lead only to approximate cancellation
when the number of points is relatively big and hence, it is very difficult to
obtain accurate results with affordable sample sizes. This non-perfect
cancellation of negative and positive parts in the integration method, usually
resulting in large quadrature rule errors, is called the _sign-problem_. The
sign-problem is for example the reason why important physical interactions in
the early universe cannot be simulated which could explain why there is more
matter than anti-matter in our universe today. To acquaint better knowledge of
these fundamental phenomena, it is essential to develop alternative quadrature
rules to MC that avoid the sign-problem.

In physical applications, the function to-be-integrated describes some
characteristic in a given physical model. We investigated methods that use
some symmetry of the physical model to result in the exact cancellation of
positive and negative parts in the quadrature rule. If the model behind the
function in Figure 1b has a reflection symmetry, few MC sampling points - in
black - can be chosen and together with their reflected - white - points they
form a set of sampling points that results in an exact quadrature rule. In
this specific example even one MC point with its reflection point would give
an exact result, for more complicated functions more sampling points are
needed.

(a) MC sampling points

(b) Symmetric sampling points

Figure 1: MC integration of an oscillatory function results in large errors,
known as the sign-problem. This problem is due to the non-cancellation of
positive and negative contributions to the quadrature rule (a). Choosing
sampling points by using the symmetry of the underlying model results in an
exact quadrature rule (b).

This article first gives a short introduction to the high-dimensional
integrals that have to be solved in particle physics, more precisely in
lattice QCD. Readers that are mostly interested in the integration methods can
easily skip this part. The main part of this article presents the methods we
developed and tested to avoid the sign-problem for high-dimensional
integration in physical one-dimensional systems.

We found that symmetrically chosen quadrature rules can avoid the sign-problem
and can efficiently be applied also to high-dimensional integrals. These rules
can help to perform simulations in important, not-yet reachable regions in
parameter space, at least in physical one-dimensional systems so far. To apply
them to higher physical dimensions, in particular to physical four-dimensional
systems in high energy physics as lattice QCD, they clearly need to be
developed further.

##  2 Integration in lattice QCD

In theoretical physics the interaction between elemetary particles such as the
electron, is described by quantum field theories (QFT), see e.g. Peskin:1995ev
. The mathematical formalism in QFT defines particles as classical fields that
are functions in three space dimensions and one time dimension, $P(x,y,z,t)$.
Operators, $O[P]$, are functionals of these fields and describe the
interactions between them. An expectation value $A$ of this interaction or
operator $O[P]$, also called amplitude, is computed via the path integral,

$\displaystyle A=\frac{\int O[P]B[P]\,\text{d}P}{\int B[P]\,\text{d}P}.$ (1)

$\int\\!\text{d}P\,$ is the infinite-dimensional integration over all possible
states of the field $P$ in time and space. The path integral becomes a well
defined expression, if a Euclidean metric is used and the fields are defined
on a finite dimensional, discrete lattice111For an alternative definition
using the $\zeta$-regularization see Hartung:2018usn ; Hartung:2019fcf .. In
(1), $B[P]$ is called the Boltzmann-weight and provides a probability which
weights the particle (field) interactions. The denominator in (1) insures the
proper normalization of $A$. The expectation value $A$ is interesting because
physical observables can be derived from it and their numerical values can be
compared with experimental results or can give new results that are not yet
possible to reach with experiments.

In lattice field theory, space-time and the involved functionals $O[P]$ and
$B[P]$ are discretized in Euclidean space, such that (1) can be computed
numerically. Often, the Boltzmann-weight is a highly peaked function
suggesting that this computation can be done using importance sampling
techniques. In most computations, this importance sampling is done by a Markov
chain MC (MCMC) algorithm using a Markov chain that leaves the distribution
density $\frac{B[P]}{\int B[P]\text{d}P}$ invariant. To compute $A$
numerically, four-dimensional space-time is discretized on a four-dimensional
lattice with four directions $\mu\in\\{1,2,3,4\\}$ , lattice sites
$\bm{n}\in\Lambda=\\{(n_{1},n_{2},n_{3},n_{4})|n_{1},n_{2},n_{3},n_{4}\in\\{1,...,d\\}\\}$
and discretized fields $P$. This results in an $4d$-dimensional integration
over the Haar measure of the compact group $\mathcal{SU}(3)$. For real
applications, $d$ can be very large, reaching orders of magnitude of several
thousands nowadays. Thus, we are left with an extremely high dimensional
integration problem. Moreover, for some physically very important questions
MCMC methods cannot be applied succesfully. This concerns, for example, the
very early universe or the matter anti-matter asymmetry which leads to our
sheer existence. Thus, a number of interesting questions remain completely
unanswered and it is exactly here where new high dimensional integration
methods could be extremely helpful

Still, the MCM methods have led to very successful computations already. By
performing numerical computations on massively parallel super computers a very
impressive result of such a lattice MCMC can be obtained: namely, the mass
spectrum of the lightest composite particles made out of quarks and gluons
that agrees completely with the experimental values, see Figure 2. To get
similar precise results for other, more error-prone observables, research is
going on to develop new methods to make this high-dimensional integration
faster and the results more precise.

Figure 2: The via lattice QCD computed masses of different composite particles
(dots with vertical error bars) agree with the experimentally measured values
(horizontal lines with error boxes) Durr:2008zz . The masses of $\pi$, $K$ and
$\sum$ (dots without error bars) were input values to the computation.

A more detailed introduction to lattice QCD is for example given in the
textbooks Gattringer:2010zz ; DeGrand:2006zz ; Rothe:1992nt .

##  3 Quadrature rules for one-dimensional lattices

In lattice QCD, the amplitude of interactions between quarks and gluons in
physical four-dimensional space-time can be computed via a high-dimensional
integral using the Haar-measure over the compact group $\mathcal{SU}(3)$, see
section 2. This integral is typically solved numerically using MCMC methods.
If the integrand is an oscillatory function, this method results in the sign-
problem that gives large errors and avoids physical insights in important
processes. We developed alternative methods that avoid the sign-problem and at
the same time are efficient for high-dimensional integration over compact
groups. Due to various complications with physical four-dimensional lattice
QCD, we developed and tested the methods for physical one-dimensional models
that involve low-dimensional and high-dimensional integration over compact
groups. As suggested in section 1, we developed quadrature rules using the
symmetry of the models.

This section is structured from low-dimensional to high-dimensional
integration: First, it introduces symmetric quadrature rules for one-
dimensional integration over compact groups to avoid the sign-problem here.
Then, it presents the recursive numerical integration (RNI), a method to
reduce high-dimensional integrals to nested one-dimensional integrals.
Finally, it shows how to combine both methods to avoid the sign-problem for
high-dimensional integration over compact groups. For all three presented
methods, the section shows results of applying them to simple physical, one-
dimensional models. More detailed explanations of the methods and applications
can be found in Volmer2018New .

###  3.1 Avoiding the sign-problem in physical one-dimensional systems

The sign-problem can already arise in a one-dimensional integration, solving

$\displaystyle I(f)=\int_{G}f(U)\,\text{d}U$ (2)

with MC methods over the Haar-measure of
$G\in\\{\mathcal{U}(N),\mathcal{SU}(N)\\}$. Finding an alternative suitable
quadrature rule $Q(f)$ ad-hoc to approximate this integral is not
straightforward. The articles Bloch:2013ara ; Bloch:2013qva  suggest that
using symmetrically distributed sampling points can be beneficial for avoiding
the sign-problem, possibly resulting in an exact cancellation of positive and
negative contributions to the integral, as stated in section 1. The article of
Genz genz  gives efficient quadrature rules for integrations over spheres,
choosing the sampling points symmetrically on the spheres. We searched for
measure preserving homeomorphisms to apply the symmetric quadrature rules on
spheres to the integration over compact groups. This section describes the two
steps to create the symmetric quadrature rules $Q(f)$ for (2):

  1. Sym 1.

Rewrite the integral $I(f)$ over the compact group $G$ into an integral over
spheres. We restricted ourselves to
$G\in\\{\mathcal{U}(1),\mathcal{U}(2),\mathcal{U}(3),\mathcal{SU}(2),\mathcal{SU}(3)\\}$.

  2. Sym 2.

Approximate each integral over one spheres by a symmetric quadrature rule as
proposed in Genz genz , and combine them to a product rule $Q(f)$.

Finally, this section shows results of applying $Q(f)$ to the one-dimensional
QCD model with a sign-problem. A more detailed explanation of the method can
be found in Ammon:2016jap ; Ammon:2016ztz .

By finding measure preserving homeomorphisms between the compact groups and
products of spheres we created polynomially exact quadrature rules for compact
groups. The application of these rules to the one-dimensional QCD model gave
results on machine precision where the standard MC method shows a sign-
problem. Therefore the symmetric quadrature rules avoid the sign-problem and
give rise to solve integrals in beforehand non-reachable parameter regions.

#### Sym 1. Rewriting the integral

The symmetric quadrature rules of Genz genz  are designed for the integration
over $k$-dimensional spheres $S^{k}$. To use them for the integration over the
compact groups $\mathcal{U}(N)$ and $\mathcal{SU}(N)$ with $N\in\\{2,3\\}$ in
(2), the compact groups have to be associated with spheres. The facts that
$\mathcal{U}(N)$ is isomorphic to the semidirect product of $\mathcal{SU}(N)$
acting on $\mathcal{U}(1)$ $\big{(}U(N)\cong SU(N)\rtimes U(1)\big{)}$, that
$\mathcal{U}(1)$ is isomorphic to $S^{1}$ $\big{(}\mathcal{U}(1)\cong
S^{1}\big{)}$ and that $\mathcal{SU}(N)$ is a principal $\mathcal{SU}(N-1)$
bundle over $S^{2N-1}$ result in

$\displaystyle\mathcal{SU}(N)$ $\displaystyle\simeq S^{3}\times
S^{5}\times...\times S^{2N-1},$ (3) $\displaystyle\mathcal{U}(N)$
$\displaystyle\simeq S^{1}\times S^{3}\times...\times S^{2N-1}.$ (4)

Then, the integral over the Haar-measure of $G$ in (2) can be rewritten as the
integral over products of spheres,

$\displaystyle\int_{G}\\!\text{d}U\,f(U)=$
$\displaystyle\int_{S^{2N-1}}\Bigg{(}\int_{S^{2N-3}}\bigg{(}\cdots\int_{S^{n+2}}\Big{(}\int_{S^{n}}$
$\displaystyle\qquad\qquad
f\big{(}\Phi(\bm{x}_{S^{2N-1}},\bm{x}_{S^{2N-3}},\ldots,\bm{x}_{S^{n+2}},\bm{x}_{S^{n}})\big{)}$
(5)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\text{d}\bm{x}_{S^{n}}\Big{)}\text{d}\bm{x}_{S^{n+2}}\cdots\Big{)}\text{d}\bm{x}_{S^{2N-3}}\Big{)}\text{d}\bm{x}_{S^{2N-1}},$

with $n=1$ for $\mathcal{U}(N)$ and $n=3$ for $\mathcal{SU}(N)$ Ammon:2016zur
. Here, $\bm{x}_{S^{k}}$ is an element on the $k$-sphere and $\Phi:\
\bigtimes_{j}S^{2j-1}\rightarrow G$ with
$G\in\\{\mathcal{U}(N),\mathcal{SU}(N)\\}$ is a measure preserving
homeomorphism. We found the homeomorphisms $\Phi_{G}\equiv\Phi$ for the
compact groups
$G\in\\{\mathcal{U}(1),\mathcal{U}(2),\mathcal{U}(3),\mathcal{SU}(2),\mathcal{SU}(3)\\}$:

  * •

For $\mathcal{SU}(2)$, $\Phi$ is an isomorphism, given by

$\displaystyle\Phi_{\mathcal{SU}(2)}:S^{3}$
$\displaystyle\rightarrow\mathcal{SU}(2),$ $\displaystyle\bm{x}$
$\displaystyle\mapsto\begin{pmatrix}x_{1}+ix_{2}&-(x_{3}+ix_{4})^{*}\\\
x_{3}+ix_{4}&\phantom{-}(x_{1}+ix_{2})^{*}\end{pmatrix}.$ (6)

  * •

For $\mathcal{SU}(3)$, spherical coordinates of $S^{5}$ are needed,

$\displaystyle\Psi:[0,2\pi)^{3}\times[0,\frac{\pi}{2})$
$\displaystyle\rightarrow S^{5},$
$\displaystyle(\alpha_{1},\alpha_{2},\alpha_{3},\phi_{1},\phi_{2})$
$\displaystyle\mapsto\begin{pmatrix}\cos\alpha_{1}\sin\phi_{1}\\\
\sin\alpha_{1}\sin\phi_{1}\\\ \sin\alpha_{2}\cos\phi_{2}\sin\phi_{2}\\\
\cos\alpha_{2}\cos\phi_{2}\sin\phi_{2}\\\
\sin\alpha_{3}\cos\phi_{1}\cos\phi_{2}\\\
\cos\alpha_{3}\cos\phi_{1}\cos\phi_{2}\end{pmatrix}.$ (7)

Then, $\Phi$ is given by

$\displaystyle\Phi_{\mathcal{SU}(3)}:S^{5}_{1}\times S^{3}$
$\displaystyle\rightarrow\mathcal{SU}(3),$ $\displaystyle(\bm{x},\bm{y})$
$\displaystyle\mapsto A(\Psi^{-1}(\bm{x}))\cdot B(\bm{y}),$ (8)

with the matrices

$\displaystyle A(\Psi^{-1}(\bm{x}))$
$\displaystyle=\scalebox{0.7}{\mbox{$\displaystyle\begin{pmatrix}\operatorname{\mathrm{e}}^{i\alpha_{1}}\cos\phi_{1}&0&\operatorname{\mathrm{e}}^{i\alpha_{1}}\sin\phi_{1}\\\
-\operatorname{\mathrm{e}}^{i\alpha_{2}}\sin\phi_{1}\sin\phi_{2}&\operatorname{\mathrm{e}}^{-i(\alpha_{1}+\alpha_{3})}\cos\phi_{2}&\operatorname{\mathrm{e}}^{i\alpha_{2}}\cos\phi_{1}\sin\phi_{2}\\\
-\operatorname{\mathrm{e}}^{i\alpha_{3}}\sin\phi_{1}\cos\phi_{2}&-\operatorname{\mathrm{e}}^{-i(\alpha_{1}+\alpha_{2})}\sin\phi_{2}&\operatorname{\mathrm{e}}^{i\alpha_{3}}\cos\phi_{1}\cos\phi_{2}\end{pmatrix}$}},$
(9) $\displaystyle B(\bm{y})$
$\displaystyle=\begin{pmatrix}x_{1}+ix_{2}&-(x_{3}+ix_{4})^{*}&0\\\
x_{3}+ix_{4}&\phantom{-}(x_{1}+ix_{2})^{*}&0\\\ 0&0&1\end{pmatrix}.$ (10)

$\Psi^{-1}(\bm{x})$ is the inverse transformation of (7) from Euclidean to
spherical coordinates. $S^{5}_{1}$ denotes $S^{5}$ without its poles,
$\phi_{1}=0$ or $\phi_{2}=0$, because at these points the inverse
transformation is not unique. The therefore excluded set is a null set, thus
$\Phi_{\mathcal{SU}(3)}$ can still be used in (3.1).

  * •

For $\mathcal{U}(1)$, $\Phi$ is an isomorphism,

$\displaystyle\Phi_{\mathcal{U}(1)}:S^{1}$
$\displaystyle\rightarrow\mathcal{U}(1),$ $\displaystyle\alpha$
$\displaystyle\mapsto\operatorname{\mathrm{e}}^{i\alpha},$ (11)

with $\alpha\in[0,2\pi)$.

  * •

For $\mathcal{U}(2)$, $\Phi$ is an isormophism,

$\displaystyle\Phi_{\mathcal{U}(2)}:S^{3}\times S^{1}$
$\displaystyle\rightarrow\mathcal{U}(2),$ $\displaystyle(\bm{x},\alpha)$
$\displaystyle\mapsto\Phi_{\mathcal{SU}(2)}(\bm{x})\cdot\text{diag}{(}\operatorname{\mathrm{e}}^{i\alpha},1).$
(12)

  * •

For $\mathcal{U}(3)$, $\Phi$ is given by

$\displaystyle\Phi_{\mathcal{SU}(3)}:S^{5}_{1}\times S^{3}\times S^{1}$
$\displaystyle\rightarrow\mathcal{U}(3),$
$\displaystyle(\bm{x},\bm{y},\alpha)$
$\displaystyle\mapsto\Phi_{\mathcal{SU}(3)}(\bm{x},\bm{y})\cdot\text{diag}{(}\operatorname{\mathrm{e}}^{i\alpha},1,1).$
(13)

#### Sym 2. Quadrature rule for spheres

With the measure preserving homeomorphism $\Phi$ in section 3.1, the integral
(2) can be written as an integral over a product of spheres as in (3.1). To
approximate the full integral numerically, one can use a product quadrature
rule with quadratures $Q_{S^{k}}(g)$ that are specifically designed for
integrations over spheres. The full integral can be computed efficiently if
the number of involved spheres is small. As pointed out in the last
subsection, in practice we are interested to build product rules for at most
$S^{5}_{1}\times S^{3}\times S^{1}$. The quadratures over each sphere can be
built in many ways. Since we are aiming for resulting quadratures that exhibit
some symmetry characteristics to hopefully overcome the sign-problem, it seems
that quadrature rules given in genz  exhibit all requiered properties, i.e.
high accuracy due to polynomial exactness over spheres, numerical stability of
the resulting weights, and beeing fully symmetric. The quadratures over each
sphere take the form

$\displaystyle
Q_{S^{k}}(g)=\sum_{\gamma=1}^{N_{\text{sym}}}w_{\gamma}\,\,g(\bm{t}_{\gamma}).$
(14)

The sampling points $\bm{t}\in S^{k}$ are chosen symmetrically on the
$k$-sphere and are weighted via $w\in{\mathbb{R}}$. The specific definitions
of $\bm{t}$, $w$ and $N_{\text{sym}}$ for different $k$ are given in genz .
(Note that in this reference, the notation $U_{k}$ is equivalent to the here
used $S^{k-1}$.) It is possible to randomize these quadrature rules, such that
an error estimate for each quadrature rule can be computed via independent
replication genz .

The final quadrature rule $Q(f)$ of the full integral in (3.1) is a
combination of different single-sphere quadrature rules given in (14). Due to
the symmetric choice of the sampling points on spheres, the rule $Q(f)$ is in
the following called _symmetrized quadrature rule_. A more detailed
description of $Q_{S^{k}}(g)$ and $Q(f)$ is given in Volmer2018New , section
6.1.

#### Application to one-dimensional QCD

We applied these constructed quadrature rules to physical one-dimensional QCD
problems Bilic:1988rw , which is a simplified model of strong interactions in
elementary particle physics. This model is a good test model because it can be
solved analytically, giving a well defined measure for the uncertainties
computed by different numerical integration methods. This model has one
integration variable $U\in G$ and three real input parameters: a mass $m$, a
chemical potential $\mu$ and a length scale $d$. A small mass ($m\ll d\mu$)
introduces a sign-problem which makes it very hard for standard methods as MC
to compute amplitudes as in (1) numerically.

We computed the chiral condensate in this model, given by

$\displaystyle\chi=\frac{\int_{G}\partial_{m}B[U]\,\text{d}U}{\int_{G}B[U]\,\text{d}U},$
(15)

with the Boltzmann-weight

$\displaystyle
B[U]=\det\left(c_{1}(m)+c_{2}(d,\mu)U^{\dagger}+c_{3}(d,\mu)U\right),$ (16)

expressed via the parameters

$\displaystyle c_{1}(m)$ $\displaystyle=\prod_{j=1}^{L}\tilde{m}_{j},$
$\displaystyle\tilde{m}_{1}$ $\displaystyle=m,$ $\displaystyle\tilde{m}_{j}$
$\displaystyle=m+\frac{1}{4\tilde{m}_{j-1}}\quad\forall
j\in\\{2,3,...,d-1\\},$ $\displaystyle\tilde{m}_{d}$
$\displaystyle=m+\frac{1}{4\tilde{m}_{d-1}}+\sum_{j=1}^{d-1}\frac{(-1)^{j+1}2^{-2j}}{\tilde{m}_{j}\prod_{k=1}^{j-1}\tilde{m}_{k}^{2}},$
(17) $\displaystyle c_{2}(d,\mu)$
$\displaystyle=2^{-d}\operatorname{\mathrm{e}}^{-d\mu},$ (18) $\displaystyle
c_{3}(d,\mu)$ $\displaystyle=(-1)^{d}2^{-d}\operatorname{\mathrm{e}}^{d\mu}.$
(19)

For brevity, the dependencies of these parameters are in the following only
written when needed.

In all numerical calculations, we first computed both numerator and
denominator of (15) separately and then divided them. We computed the
numerator by symbolically differentiating $B[U]$ and computing the integral
over the result numerically.

We compared the results for $\chi$ using the symmetrized quadrature rules that
are described in 3.1, with a standard integration method, ordinary MC
sampling. The latter quadrature rule is given by

$\displaystyle
Q(f)=\frac{1}{N_{\text{MC}}}\sum_{\gamma=1}^{N_{\text{MC}}}\,\,f(V_{\gamma}),$
(20)

where the $V$ are matrices that are chosen randomly from a uniform
distribution. We chose $N_{\text{MC}}$ to be as large as the number of used
symmetric sampling points.

Because the analytic results of $\chi$ can be calculated straightforwardly, we
computed the error estimates of the numerical solutions - MC and symmetrized
quadrature rules - directly via the relative deviation from the analytic
value,

$\displaystyle\Delta\chi=\frac{|\chi_{\text{numerical}}-\chi_{\text{analytic}}|}{|\chi_{\text{analytic}}|}$
(21)

and derived the standard deviation of this error by repeatedly using on the
one hand the MC quadrature rules with different random matrices $V$’s and on
the other hand the randomized symmetrized quadrature rules as indicated in
section 3.1.

The results for $\Delta\chi$ of both MC and symmetrized quadrature rule can be
roughly split into a small $m$ ($m<10^{-1}$), a large $m$ ($m>10^{0.5}$) and a
transition region, shown in Figure 3 for constant $\mu=1$ and $d=8$, extended
1024-bit machine precision and different compact groups. For both quadrature
rules we used the sampling sizes $N\equiv N_{\text{sym}}=N_{\text{MC}}=8$ for
$\mathcal{SU}(2)$, $N=96$ for $\mathcal{SU}(3)$, $N=4$ for $\mathcal{U}(1)$,
$N=32$ for $\mathcal{U}(2)$ and $N=384$ for $\mathcal{U}(3)$.

Figure 3: The sign-problem arises for MC results with small $m$ constants,
giving errors of the order of one. On the contrast, the symmetrized quadrature
rules avoid the sign-problem in this region completely, giving errors
approximately at machine precision for all shown groups.

First, we describe the MC results: In the small $m$ region, $\Delta\chi$ for
all groups are large - equal or larger than one. It can be shown that in this
region the numerator of $\chi$ is such small that the MC evaluation cannot
resolve these values for affordable sample sizes, resulting in large errors
Volmer2018New . This is the manifestation of the sign-problem, making it
almost impossible to compute reasonable values of $\chi$ with MC in the small
$m$ region. On the other side, for large $m$ all groups have a smaller MC
error estimate than in the small $m$ region. Here the numerator of $\chi$
tends to be larger and especially the denominator becomes very large, both
resulting in a slightly better error estimate for the MC results.

Opposed to MC results, the symmetrized quadrature rules give error estimates
approximately at machine precision up to very small $m$ values, see Figure 3.
These numerical results show that the symmetrized quadrature rules give
significant results in the sign-problem region in practice, where MC
simulations have error estimates of order one.

###  3.2 Reducing high-dimensional integrals to nested one-dimensional
integrals

The previous section shows efficient quadrature rules for physical one-
dimensional integration to avoid the sign-problem. Most physical models have
more than one integration variable. In general, it is not straightforward to
find an efficient quadrature rule, and usually restricted Monte Carlo methods
are applied to high-dimensional integrals. As a first alternative, we
investigated the recursive numerical integration(RNI) method. This method
reduces the $d$-dimensional integral

$\displaystyle I(f)=\int_{D^{d}}f[\varphi]\,\text{d}\varphi$ (22)

with $\text{d}\varphi=\prod_{i=1}^{d}\text{d}\varphi_{i}$ and $D=[0,2\pi)$
into many recursive one-dimensional integrals, and can be applied for several
physical models of interest.

This is done by utilizing the typical structure of the integrand $f[\varphi]$.
This section focuses on the RNI method and how to find an efficient quadrature
rule for a high-dimensional integral. It does not discuss the sign-problem
which is investigated further in section 3.3. More specifically, this section
describes the two steps to create an efficient quadrature rules $Q(f)$ for the
integral $I(f)$ in (22):

  1. RNI 1.

Use the structure of the integrand of the high-dimensional integral to rewrite
it into recursive one-dimensional integrals.

  2. RNI 2.

Choose an efficient quadrature rule to compute each one-dimensional integral
numerically. Recursively doing this results in the full quadrature rule
$Q(f)$.

Finally, this section shows results of applying the method to a physical model
called the topological osciallator. A more detailed explanation of the method
and the results can be found in Ammon:2015mra ; Ammon:2016dfi .

#### RNI 1. Using the structure of the integrand

Many models in one physical dimensional have integrands with the structure

$\displaystyle f[\varphi]=\prod_{i=1}^{d}f_{i}(\varphi_{i+1},\varphi_{i}),$
(23)

with periodic boundary conditions $\varphi_{d+1}=\varphi_{1}$. These models
have only next-neighbor couplings.

The integral of (23) can be rewritten using recursive integration as described
in Genz86 ; Hayter06 : Because of next-neighbor couplings, each variable
$\varphi_{i}$ appears only twice in $f[\varphi]$, in $f_{i}$ and $f_{i-1}$,
and therefore the integral can be written as $d$ nested one-variable integrals
$I_{i}$,

$\displaystyle
I(f)=\int_{D}...\int_{D}\prod_{i=1}^{d}f_{i}(\varphi_{i},\varphi_{i+1})\,\text{d}\varphi_{d}\cdots\text{d}\varphi_{1}$
(24)
$\displaystyle\leavevmode\resizebox{422.77661pt}{}{$\displaystyle{=\underbrace{\int_{D}\Bigg{(}...\underbrace{\bigg{(}\int_{D}f_{d-2}(\varphi_{d-2},\varphi_{d-1})\cdot\underbrace{\bigg{(}\int_{D}f_{d-1}(\varphi_{d-1},\varphi_{d})\cdot
f_{d}(\varphi_{d},\varphi_{d+1})\,\text{d}\varphi_{d}\bigg{)}}_{I_{d}}\,\text{d}\varphi_{d-1}\bigg{)}}_{I_{d-1}}\cdots\Bigg{)}\,\text{d}\varphi_{1}}_{I_{1}}}$}.$

This full integral can be computed recursively: $I_{d}$ integrates out
$\varphi_{d}$ first, then $I_{d-1}$ integrates out $\varphi_{d-1}$ and so on
until finally $I_{1}=I(f)$ integrates out $\varphi_{1}$.

To avoid under- and overflow of the single quadrature rule results, we
actually used quadrature rules to approximate $I_{i}^{*}=\frac{1}{c_{i}}I_{i}$
with $c_{i}>0$ chosen adaptively. Then, the final integral is computed via
$I=\left(\prod_{i=1}^{d}c_{i}\right)I^{*}$. For brevity, the method is
described in the following without this trick.

Each integral is approximated by using an $N_{\text{quad}}$-point quadrature
rule. The first integrand in (24) (last from the right) depends on three
variables $\varphi_{d-1}$, $\varphi_{d}$ and $\varphi_{d+1}$. The variable
$\varphi_{d}$ is integrated out, therefore the quadrature rule
$Q_{d}(f_{d-1}\cdot f_{d})\equiv Q_{d}$ of $I_{d}$ depends on two variables,

$\displaystyle
Q_{d}(\varphi_{d-1},\varphi_{d+1})=\sum_{\gamma=1}^{N_{\text{quad}}}w_{\gamma}\,\,f_{d-1}(\varphi_{d-1},t_{\gamma})\,\,f_{d}(t_{\gamma},\varphi_{d+1}),$
(25)

with sampling points $t$ and weights $w$. The next integral $I_{d-1}$ is
approximated by the quadrature rule

$\displaystyle
Q_{d-1}(\varphi_{d-2},\varphi_{d+1})=\sum_{\gamma=1}^{N_{\text{quad}}}w_{\gamma}\,\,f_{d-2}(\varphi_{d-2},t_{\gamma})\,\,Q_{d}(t_{\gamma},\varphi_{d+1}),$
(26)

and includes the quadrature rule $Q_{d}$ given in (25). The quadrature rules
$Q_{d-2}$, …, $Q_{1}$ are created analogically to (26). Using the same
sampling points $w_{\gamma}$ and weights $t_{\gamma}$,
$\gamma\in\\{1,...,N_{\text{quad}}\\}$ in all quadrature rules $Q_{i}$ results
in the full quadrature rule for (24),

$\displaystyle
Q=Q_{1}=\sum_{\gamma=1}^{N_{\text{quad}}}w_{\gamma}Q_{2}(t_{\gamma},t_{\gamma})=\text{tr}{\left[{{\prod_{i=1}^{d}\left(M_{i}\cdot\text{diag}{(w_{1},\dots,w_{N_{\text{quad}}})}\right)}}\right]},$
(27)

with $M_{i}$ beeing an $N_{\text{quad}}\times N_{\text{quad}}$ matrix with
entries $(M_{i})_{\alpha\beta}=f_{i}(t_{\alpha},t_{\beta})$.

#### RNI 2. Choosing an efficient quadrature rule

We used the Gaussian-Legendre $N_{\text{quad}}$-point quadrature rule, see
stoer2013introduction  to define the sampling points $t$ and weights $w$. For
this rule, the error scales asymptotically (for large $N_{\text{quad}}$) as
$\sigma\sim\mathcal{O}\left(\frac{1}{(2N_{\text{quad}})!}\right)$ (For
Legendre polynomials the correct asymptotic error scaling is
$\frac{(N_{\text{quad}}!)^{4}}{((2N_{\text{quad}})!)^{3}}$
kahaner1989numerical  which is slightly improved over
$\frac{1}{(2N_{\text{quad}})!}$.). The Stirling formula
($N_{\text{quad}}!\approx\sqrt{2\pi
N_{\text{quad}}}\left(\frac{N_{\text{quad}}}{e}\right)^{N_{\text{quad}}}$
asymptotically) approximates the factorial to give

$\displaystyle\sigma\sim\mathcal{O}\left(\exp(-2N_{\text{quad}}\ln
N_{\text{quad}})\frac{1}{\sqrt{N_{\text{quad}}}}\right)$ (28)

asymptotically. This is a huge improvement over the MC error scaling
$1/\sqrt{N_{\text{MC}}}$.

#### Application to the topological oscillator

We applied the RNI method to the topological oscillator Bietenholz:2010xg ,
also called quantum rotor, which is a simple, physically one-dimensional model
that has non-trivial characteristics which are also present in more complex
models. It has $d$ variables $\phi_{i}\in[0,2\pi)$, a length scale $T$ and a
coupling constant $c$. We investigated the topological charge susceptibility
of this model,

$\displaystyle\chi_{\text{top}}$ $\displaystyle=\frac{\int
O[\varphi]B[\varphi]\,\text{d}\varphi}{\int B[\varphi]\,\text{d}\varphi},$
(29)

with Boltzmann-weight

$\displaystyle B[\varphi]$
$\displaystyle=\exp{\left(-c\sum_{i=1}^{d}\left(1-\cos(\varphi_{i+1}-\varphi_{i})\right)\right)},$
(30)

and a squared topological charge

$\displaystyle O[\varphi]$
$\displaystyle=\frac{1}{T}\left(\frac{1}{2\pi}\sum_{i=1}^{d}(\varphi_{i+1}-\varphi_{i})\mod
2\pi\right)^{2}.$ (31)

With RNI, we computed both numerator and denominator of $\chi_{\text{top}}$
separately, both differing in the factorization (23) of their integrands.
Straightforwardly, the denominator integrand consists out of local exponential
factors. The numerator consists out of summands with varying factorization
schemes, each of these summands is computed separately with RNI and they are
presented in more detail in Volmer2018New , section 5.2. We estimated the
error of $\chi_{\text{top}}$ by choosing a large number of samples
$N_{\text{quad}}^{g}$ in (25), (26) and similar ones for which we assumed that
$\chi_{\text{top}}(N_{\text{quad}}^{g})$ has converged to the actual value and
computed the difference of $\chi_{\text{top}}(N_{\text{quad}})$ for
$N_{\text{quad}}<N_{\text{quad}}^{g}$ to this value,

$\displaystyle\Delta\chi_{\text{top}}(N_{\text{quad}})=|\chi_{\text{top}}(N_{\text{quad}})-\chi_{\text{top}}(N_{\text{quad}}^{g})|.$
(32)

We tested beforehand that this truncation error behaves exponentially for
large $N_{\text{quad}}$ in practice, as expected from (28), Ammon:2015mra .

We compared the results of the RNI method with results using the Cluster
algorithm Wolff:1988uh , which we found is an optimal MCMC method for the
application to the topological oscillator Ammon:2015mra . Due to the
exponential error scaling of the Gauss-Legendre rule, the new method advances
MCMC for large enough $N_{\text{quad}}$. We found that the RNI method is also
advantageous for lower $N_{\text{quad}}$-values: our simulations showed that
the RNI method needs orders of magnitude less runtime than the Cluster
algorithm to result in a specified error estimate on an observable, compare
Figure 4 for $c=2.5$, $T=20$, $d=200$.

Figure 4: The runtime to arrive at a given error estimate is orders of
magnitudes smaller when using the RNI method with Gauss-Legendre points than
using the Cluster MCMC algorithm.

The Cluster algorithm measurements resulted in an error estimate that
decreases proportional to $t^{-1/2}$ for runtime $t$, consistent with the
typical MC error scaling probTheory . We used between $10^{2}$ and $10^{6}$
sampling points here. The RNI method, using between 10 and 300 sampling points
with $N_{\text{quad}}^{g}=400$, resulted in orders of magnitude smaller
errors. The exponential error scaling in (28) is not visible here, the
asymptotic regime of the method is not yet reached with the used numbers of
sampling points.

All in all, the RNI method results in orders of magnitude smaller errors than
the Cluster algorithm for a fixed runtime or equivalently, the RNI method
needs orders of magnitude less runtime than the Cluster algorithm to arrive at
a fixed error estimate, even for a number of sampling points where the RNI
error does not yet scale exponentially.

###  3.3 Avoiding the sign-problem in high-dimensional integrals

Section 3.1 shows that the sign-problem can be avoided for one-dimensional
integrals using symmetric quadrature rules. But what about the sign-problem
for high-dimensional integrals? A quadrature rule for high-dimensional
integrals over compact groups,

$\displaystyle I(f)=\int_{G^{d}}f[U]\,\text{d}U,$ (33)

with $\text{d}U=\prod_{i=1}^{d}\text{d}U_{i}$ is needed that also avoids the
sign-problem. We combined both already presented methods, the symmetric
quadrature rules in section 3.1 and the RNI in section 3.2 to find an
efficient quadrature rule $Q(f)$ for $I(f)$ in (33). An alternative attempt to
generalize the symmetrized quadrature rules to high-dimensional integrals is
discussed in 10.1007/978-3-319-75996-8_15 .

#### Combining recursive numerical integration and symmetric quadrature rules

RNI can be used to transform the high-dimensional integral $I(f)$ in (33) into
one-dimensional integrals. These one-dimensional integrals can be approximated
recursively, using the symmetric quadrature rules. In the following, these
steps are described in more detail:

  1. RNI 1.

Find the structure, i.e. all $f_{i}$, of the integrand

$\displaystyle f[U]=\prod_{i=1}^{d}f_{i}(U_{i+1},U_{i}),$ (34)

to be able to write the full integral as nested one-dimensional integrals,
similar to (24).

  2. RNI 2.

Apply symmetric quadrature rules to each one-dimensional integration over
$U_{i}$. Here is an example how to do this for the innermost integral $I_{d}$,
integrating over $U_{d}$:

    1. Sym 1.

Rewrite the integral over $U_{d}$ into an integral over the products of
spheres as done in (3.1).

    2. Sym 2.

Approximate each iterated integral $I_{d}(g)$ by a product rule of quadratures
over spheres parametrising the group $U_{d}$ to be integrated. Note that the
group $U_{d}$ is parametrised at most as the product of $S^{1}$, $S^{3}$ and
$S^{5}$.

#### Application to topological oscillator with sign-problem

We applied this combined method again to the topological oscillator discussed
in 3.2. This time we transformed the variables $\varphi_{i}$ to new variables
$U_{j}=\operatorname{\mathrm{e}}^{i\varphi_{j}}\in\mathcal{U}(1)$.
Additionally, we added a sign-problem to the model by using an additional
factor $\prod_{j=1}^{d}U_{j}^{-\theta}$ in the Boltzmann-weight,

$\displaystyle B[U]$
$\displaystyle=\exp{\left(-c\sum_{i=1}^{d}\Re(1-U_{i+1}U_{i}^{*})\right)}\cdot\prod_{j=1}^{d}U_{j}^{-\theta},$
(35)

with a new parameter $\theta\in{\mathbb{R}}$. If this parameter is larger than
zero, the sign-problem arises and is most severe for $\theta=\pi$.

In this model we computed the plaquette,

$\displaystyle plaquette$ $\displaystyle=\frac{\int O[U]B[U]\,\text{d}U}{\int
B[U]\,\text{d}U},$ (36)

with

$\displaystyle
O[U]=\frac{1}{d}\Re\left(\sum_{i=1}^{d}U_{i+1}U_{i}^{*}\right).$ (37)

For the combined method, we computed both numerator and denominator of (36)
separately and divided the values. We used a truncation error, similar to the
one given in (32). We compared the method with a standard MC method as used in
3.1. The MC error is computed via the standard deviation.

For $\theta=\pi$ we found that the combined method avoids the sign-problem
that is visible with the MC computation, compare Figure 5. It gives orders of
magnitude smaller errors that shrink the more symmetrization points are used.

Figure 5: The combined method avoids the sign-problem that exists when using
the MC method.

Therefore the combination of RNI and symmetric quadrature rules is suitable to
avoid the sign-problem for high-dimensional integration.

##  4 Conclusion

In this contribution we have demonstrated that through symmetric quadrature
rules exact symmetrization and recursive numerical integration techniques
problems in high energy physics can be solved which constitute a major, if not
unsurmountable obstacle for standard Markov chain Monte Carlo methods. The
examples we have considered here invole only a time lattice and are hence
0+1-dimensional in space-time, where as real physical problem include spatials
dimensions of up to 3. We are presently investigating whether the methods we
have presented here can be extended to higher, i.e. including also spacial,
dimensions. While for the recursive numerical integration technique we have
first results which are promising, for the full symmetrization method we were
so far not successful.

Also combining the symmetrized quadrature rules with MC methods did not lead
to a practically feasible method in higher dimensions. However, we are
following a path to combine Quasi Monte Carlo, recursive numerical integration
and a full symmetrization to overcome this problem and hope to report about
these attempts in the future.

## References

  * (1) Ammon, A., Genz, A., Hartung, T., Jansen, K., Leövey, H., Volmer, J.: Applying recursive numerical integration techniques for solving high dimensional integrals.  PoS LATTICE2016, 335 (2016)
  * (2) Ammon, A., Genz, A., Hartung, T., Jansen, K., Leövey, H., Volmer, J.: On the efficient numerical solution of lattice systems with low-order couplings.  Comput. Phys. Commun. 198, 71–81 (2016).  DOI 10.1016/j.cpc.2015.09.004
  * (3) Ammon, A., Hartung, T., Jansen, K., Leövey, H., Volmer, J.: New polynomially exact integration rules on $U(N)$ and $SU(N)$.  PoS LATTICE2016, 334 (2016)
  * (4) Ammon, A., Hartung, T., Jansen, K., Leövey, H., Volmer, J.: New polynomially exact integration rules on U(N) and SU(N) (2016).  URL https://inspirehep.net/record/1490030/files/arXiv:1610.01931.pdf
  * (5) Ammon, A., Hartung, T., Jansen, K., Leövey, H., Volmer, J.: Overcoming the sign problem in one-dimensional QCD by new integration rules with polynomial exactness.  Phys. Rev. D94(11), 114,508 (2016).  DOI 10.1103/PhysRevD.94.114508
  * (6) Bietenholz, W., Gerber, U., Pepe, M., Wiese, U.J.: Topological Lattice Actions.  JHEP 1012, 020 (2010).  DOI 10.1007/JHEP12(2010)020
  * (7) Bilic, N., Demeterfi, K.: One-dimensional QCD With Finite Chemical Potential.  Phys. Lett. B212, 83–87 (1988).  DOI 10.1016/0370-2693(88)91240-3
  * (8) Bloch, J., Bruckmann, F., Wettig, T.: Subset method for one-dimensional QCD.  JHEP 10, 140 (2013).  DOI 10.1007/JHEP10(2013)140
  * (9) Bloch, J., Bruckmann, F., Wettig, T.: Sign problem and subsets in one-dimensional QCD.  PoS LATTICE2013, 194 (2014)
  * (10) Durr, S., et al.: Ab-Initio Determination of Light Hadron Masses.  Science 322, 1224–1227 (2008).  DOI 10.1126/science.1163233
  * (11) Fristedt, B., Gray, L.: A Modern Approach to Probability Theory.  Probability and Its Applications. Birkhäuser Basel (1997)
  * (12) M. E. Peskin and D. V. Schroeder, An Introduction to quantum field theory.
  * (13) Gattringer, C., Lang, C.B.: Quantum chromodynamics on the lattice.  Lect. Notes Phys. 788, 1–343 (2010).  DOI 10.1007/978-3-642-01850-3
  * (14) T. DeGrand and C. E. Detar, Lattice methods for quantum chromodynamics  New Jersey, USA: World Scientific (2006) 345 p
  * (15) H. J. Rothe, Lattice gauge theories: An Introduction.  World Sci. Lect. Notes Phys. 43 (1992) 1  [World Sci. Lect. Notes Phys. 59 (1997) 1]  [World Sci. Lect. Notes Phys. 74 (2005) 1]  [World Sci. Lect. Notes Phys. 82 (2012) 1].
  * (16) Genz, A.: Fully Symmetric Interpolatory Rules for Multiple Integrals over Hyper-Spherical Surfaces.  Journal of Computational and Applied Mathematics 157, 187–195 (2003).  DOI 10.1016/S0377-0427(03)00413-8
  * (17) Genz, A., Kahaner, D.K.: The numerical evaluation of certain multivariate normal integrals.  J. Comput. Appl. Math. 16(2), 255–258 (1986).  DOI 10.1016/0377-0427(86)90100-7.  URL http://www.sciencedirect.com/science/article/pii/0377042786901007
  * (18) Hartung, T., Jansen, K., Leövey, H., Volmer, J.: Improving monte carlo integration by symmetrization.  In: A. Böttcher, D. Potts, P. Stollmann, D. Wenzel (eds.) The Diversity and Beauty of Applied Operator Theory, pp. 291–317. Springer International Publishing, Cham (2018)
  * (19) Hayter, A.: Recursive integration methodologies with statistical applications.  Journal of Statistical Planning and Inference 136(7), 2284–2296 (2006).  DOI 10.1016/j.jspi.2005.08.024.  URL http://www.sciencedirect.com/science/article/pii/S0378375805002223.  In Memory of Dr. Shanti Swarup Gupta
  * (20) Kahaner, D., Moler, C., Nash, S., Forsythe, G.: Numerical methods and software.  Prentice-Hall series in computational mathematics. Prentice Hall (1989)
  * (21) Stoer, J., Bartels, R., Gautschi, W., Bulirsch, R., Witzgall, C.: Introduction to Numerical Analysis.  Texts in Applied Mathematics. Springer New York (2013)
  * (22) Volmer, J.L.: New attempts for error reduction in lattice field theory calculations.  Ph.D. thesis, Humboldt-Universität zu Berlin, Mathematisch-Naturwissenschaftliche Fakultät (2018).  DOI http://dx.doi.org/10.18452/19350
  * (23) Wolff, U.: Collective Monte Carlo Updating for Spin Systems.  Phys. Rev. Lett. 62, 361 (1989).  DOI 10.1103/PhysRevLett.62.361
  * (24) T. Hartung and K. Jansen, Quantum computing of zeta-regularized vacuum expectation values  arXiv:1808.06784 [quant-ph].
  * (25) T. Hartung and K. Jansen, Integrating Gauge Fields in the $\zeta$-formulation of Feynman’s path integral  arXiv:1902.09926 [math-ph].


# Wakefield damping in a distributed coupling LINAC

Evan Ericson  CERN, European Organization for Nuclear Research, Geneva 1221,
Switzerland  Department of Physics and Engineering Physics, University of
Saskatchewan, Saskatoon, Canada.  Alexej Grudiev<EMAIL_ADDRESS>CERN,
European Organization for Nuclear Research, Geneva 1221, Switzerland  Drew
Bertwistle  Department of Physics and Engineering Physics, University of
Saskatchewan, Saskatoon, Canada.  Canadian Light Source, Saskatoon, Canada
Mark J. Boland  Department of Physics and Engineering Physics, University of
Saskatchewan, Saskatoon, Canada.  Canadian Light Source, Saskatoon, Canada

###### Abstract

The number of cells in a $\pi$-mode standing wave (SW) accelerating structure
for the Compact linear Collider (CLIC) project is limited by mode overlap with
nearby modes. The distributed coupling scheme avoids mode overlap by treating
each cell as independent. Designs of cells suitable for distributed coupling
with strong wakefield damping have not previously been studied. In this paper
we develop a SW cell to be used in a distributed coupling structure that can
satisfy the CLIC transverse wakepotential limit. From the middle cell of the
CLIC-G* travelling wave (TW) structure, a SW cell is designed. The cell is
adapted to be suitable for distributed coupling. Its wakepotentials in an
ideal case of open boundaries are reduced to satisfy the wakepotential
threshold. An electric boundary is added to the model to simulate total
reflection at the distribution network. A horizontal coupler cell that
connects to the distribution network such that the reflected wakefields remain
similar to the open boundary case is simulated. A triplet module which takes
advantage of cell-to-cell coupling to reduce reflected wakepotential is
presented.

††preprint: APS/123-QED

##  I Introduction

The main linear accelerator (LINAC) of the CLIC project uses accelerating
structures to generate high-energy particles for particle physics experiments
[1]. The current LINAC design consists of TW X-band radio frequency (RF)
accelerating structures. These structures consist of 26 regular accelerating
cells and two coupler cells. The cell irises radii are tapered from
$3.15\text{\,}\mathrm{mm}$ to $2.35\text{\,}\mathrm{mm}$ to maintain a
constant loaded gradient of
$100\text{\,}\mathrm{M}\mathrm{V}\mathrm{/}\mathrm{m}$. All cells feature
waveguides with silicon carbide loads to dissipate power from higher order
modes (HOM) [2]. Figure 1 shows an example of a TW cell.

Figure 1: Accelerating cell from the CLIC-G* TW structure for CLIC with
quarter symmetry.

Recent work has shown converting the $2\pi/3$-mode TW structure to a
$\pi$-mode SW structure of equal length gives an increase in rf-to-beam
efficiency of $\sim$ $3\text{\,}\mathrm{\char 37\relax}$ [3]. A SW structure
operating in the $\pi$-mode is limited to the number of cells, N, given by
Equation 1.

$\sqrt{\frac{Q_{\pi}\pi^{2}k}{4}}>N.$ (1)

For a SW cell based on CLIC-G* cell shape, the maximum number of cells before
mode overlap is 9. With this many cells, the $\pi$-mode is separated from the
next adjacent mode by

$\Delta
f^{\pi}_{N,N-1}=f_{r}\frac{k\pi^{2}}{4N^{2}}=$2\text{\,}\mathrm{MHz}$.$ (2)

Operating a side-coupled structure in the $\pi/2$-mode simultaneously
increases the number of cells in a structure before mode overlap while keeping
the shunt impedance high. A side-coupled geometry is not easily amenable to
cells with four HOM waveguides. Distributed coupling topologies also allow for
a SW structure to increase the number of cells by treating each cell as
independent. In this configuration, each accelerating cell is connected
directly to the power source by a waveguide which runs parallel to the
structure [4, 5]. Because the structure does not rely on coupling between
cells, the cells’ iris aperture can be made small to increase the shunt
impedance. Distributed coupling schemes can also increase high-gradient
performance by limiting the amount of power flow through the cell irises which
are prone to RF breakdown events [6, 7]. As the iris aperture decreases, the
effects of wakefields increase [8]. Existing designs reduce wakefields by use
of detuning such that problematic dipole modes add coherently later, ideally
after the end of the bunch train [9]. Detuning can provide dipole modes with
Q${}_{ext}=1000$ while CLIC requires modes with Q${}_{ext}\sim 10$ [10, 11].
In this work, we present the development of an RF accelerating cell suitable
for a distributed coupling structure that satisfies the CLIC transverse
wakefield damping requirement. Section II outlines the procedure followed to
get a SW cell design. The cell is adapted to include a power input and then
further to reduce the wakepotential of the cell. In Section III the reflected
wakepotential from a distribution network is analyzed and Section IV describes
cell designs capable of reducing the reflected wakepotentials to an acceptable
level. The final findings are summarized in V.

##  II Open boundary wakes

The middle cell of the TW structure was chosen as a starting point for the
design. The cell length was adjusted so the $\pi$-mode was synchronous with a
relativistic beam. Parameters sweeps of the iris thickness and ellipticity
were performed in CST’s eigenmode solver [12] to select an iris shape that
minimized the surface electric field. Sweeps of the outer wall parameters A0
and A2, described in [2], confirmed the values of the TW cell also minimized
the surface magnetic field in this case. The outer wall was adjusted to obtain
the operating frequency. The RF parameters of the resulting cell are listed in
Table 1. An accelerating structure based on this SW cell design having the
same length as the TW structure would have an RF-to-beam efficiency of
$31\text{\,}\mathrm{\char 37\relax}$.

Table 1: RF parameters of SW cell. Quantity | Unit | Value
---|---|---
cell outer radius, b | $\text{\,}\mathrm{mm}$ | 8.5638
iris thickness, d | $\text{\,}\mathrm{mm}$ | 2.6
iris ellipticity, e | $\text{\,}\mathrm{mm}$ | 2.4
f0 | $\text{\,}\mathrm{GHz}$ | 11.9945
Q0 | - | 6,927
R/Q0 | $\text{\,}\mathrm{\SIUnitSymbolOhm}\mathrm{/}\mathrm{m}$ | 10,962
R | $\text{\,}\mathrm{M\SIUnitSymbolOhm}$ | 75.93
Es/Ea | - | 1.80
Hs/Ea | $\text{\,}\mathrm{m}\mathrm{A}\mathrm{/}\mathrm{V}$ | 4.67
Sc/E${}_{a}^{2}$ | $\text{\,}\mathrm{m}\mathrm{A}\mathrm{/}\mathrm{V}$ | 0.35
Figure 2: Long slot cell geometry with half symmetry. The cell is matched by
adjusting the cell outer radius and the length of the HOM waveguide between
the cell and the coupler waveguide.

A cell with four HOM waveguides is fed by cell-to-cell coupling. To make the
cell suitable for a distributed coupling structure, a coupler cell was
generated based on the TW coupling cell. The TW coupler cells have only two
HOM waveguides and therefore the expected wakepotential of a structure of
repeating coupler cells was higher than the four HOM waveguide cell. The SW
coupler cell included only one coupler waveguide and three HOM waveguides. The
cell was matched to the source by adjusting the cell’s outer radius and the
aperture between the coupler waveguide and the cell. The x-polarization of the
wakepotential was found to be lower in amplitude than the y-polarization. The
evaluated transverse wakepotential of the coupler cell wacs above the
threshold of
$6\text{\,}\mathrm{V}\mathrm{(}\mathrm{p}\mathrm{C}\cdot\mathrm{m}\cdot\mathrm{m}\mathrm{m}\mathrm{)}^{-1}$
[13]. A hybrid cell design with a HOM waveguide transition between the coupler
waveguide and the cell shown in Figure 2 had a transverse wakepotential
between the coupler cell and four HOM waveguide cell at the second bunch
position.

Figure 3: Transverse wakepotential envelopes of long slot cell with various
HOM aperture widths, iw (top) and various HOM waveguide widths, w (bottom).

To decrease the wakepotential below the threshold, two approaches were taken.
First, the aperture of the HOM waveguides were increased. Changing the
geometry in this way caused the Qext of the dominant dipole mode to decrease
at the expense of a lower Q0 and a larger expected temperature rise on the
surface of the cell. The cells needed to be tuned for each simulated value.
The second approach was to increase the HOM waveguide widths. The
wakepotentials from these sweeps are shown in Figure 3. The wakepotential at
the second bunch location was more sensitive to changes in w. A value of w =
$10.7\text{\,}\mathrm{mm}$ was selected. Both polarization of the transverse
wakepotentials and their impedances are shown in Figure 4.

Figure 4: Transverse wakepotentials (top) and impedances (bottom) of long slot
cell w = $10.7\text{\,}\mathrm{mm}$ (bottom).

##  III Reflected wakes

To model the effect of a distribution network, parallel waveguide or
otherwise, on the wakepotential of the structure, an electric boundary
condition was placed on the coupler waveguide surfaces. This represents a
worst-case-scenario in which all the wakefields extracted through the coupler
waveguide are reflected back into the cell where they may act on following
bunches. Figure 5 shows the geometry simulated to include the effect of
perfect reflection at the coupler waveguide.

Figure 5: Model simulated to obtain wakepotential of reflected wakefields of
the long slot cell using half symmetry. The HOM waveguide ports (red) are
connected to an open boundary. The coupler waveguides are slightly shorter
than the HOM waveguides and experience an electric boundary. The model was
parameterized to maintain these boundary conditions as the coupler waveguide
length was increased.

An example of a reflected wakepotential is shown in Figure 6. After the drive
bunch passes through the structure, the remaining wakefields are extracted
through the HOM waveguides. The fields propagating out of the cell through the
coupler waveguide encounter the electric boundary and are reflected back into
the cell. At this time, the wakepotential diverges from the open boundary
case, increasing above the wakepotential threshold. The position where the
reflected wakepotential diverges from the open boundary case is determined by
the distance between the electric boundary and the beam axis. We call this
distance the offset. The spectrogram of Figure 7 shows the mode evolution of
the reflected wakepotential. The 17 and $22\text{\,}\mathrm{GHz}$ modes of
Figure 4 are quickly damped. The reflected wakepotential first consists of
modes near $17\text{\,}\mathrm{GHz}$. While these modes are being damped,
modes near $22\text{\,}\mathrm{GHz}$ arrive. The length of the reflected wake
is due to the difference in arrival time of the two sets of modes and because
of the slow damping of these modes. Increasing the offset decreased the
maximum amplitude of the reflected wakepotential. A larger offset accentuated
the difference in arrival time causing the cell to damp one set of modes at a
time. The reflected wakepotential duration increased as a result. The modes
also spread in time due to dispersion in the waveguides.

Figure 6: Reflected transverse wakepotential envelopes of long slot cell with
w = $10.7\text{\,}\mathrm{mm}$ and offset = $150\text{\,}\mathrm{mm}$. Figure
7: Evolution of the mode content of the reflected wakepotential with offset =
$100\text{\,}\mathrm{mm}$.

The electric boundary, position of the distribution network relative to the
beam, can be adjusted such that the reflected wake arrives after the a
particular bunch but because the reflected wakepotential is damped slowly and
lasts longer than the bunch spacing, the reflected wakepotential would be
above the threshold at the next bunch. Increasing the HOM waveguide aperture
and width, which previously had an effect on the wakepotential, proved
ineffective at reducing the reflected wakepotential. The strategy of
evaluating different iris geometries in an attempt to change the mode beating
to form a minimum at the bunch located in the reflected wakepotential wasof
unsuccessful.

##  IV Designs to reduce reflected wakefields

Two configurations were found to reduce the reflected wakepotentials to an
acceptable level. First, the horizontal coupler cell design shown in Figure 8
feeds power to the cell from a coupler waveguide through an aperture in one of
the HOM waveguides. The aperture’s dimensions and location allow power to
enter the cell to establish the accelerating field while being sufficiently
small that the wakefields couple weakly to the coupler waveguide. The
wakefields are damped in the HOM waveguides instead of being reflected at the
coupler waveguide. The reflected wakepotentials of Figure 9 do not depend
strongly on the offset meaning the horizontal coupler is transparent to the
wakefields and a distribution network can be designed separately.

Figure 8: Horizontal coupler cell geometry. Figure 9: Reflected transverse
wakepotential of horizontal coupler cell for various offsets.

A second scheme where modules were formed of coupler cells and four HOM
waveguide cells like the triplet shown in Figure 10 also reduced the reflected
wakepotentials below the CLIC threshold. The reflected wakepotential
amplitudes were lowered by reducing the number of coupler waveguides from
which reflections occur and by increasing the number of HOM waveguides that
damp the wakefields. The module takes advantage of an iris size of
$2.75\text{\,}\mathrm{mm}$, which is large enough to allow for coupling
between cells. Operating a group of modules as a structure reduces the total
number of coupler waveguides required compared to distributed coupling
structure designs where each cell has a feed to a power source. Figure 11
shows the reflected wakepotential amplitude decreases as the number of four
HOM waveguide cells increases. The wakepotential of the triplet is below the
threshold at all bunch positions out to $1\text{\,}\mathrm{m}$. The outer
radii of the cells in the doublet and triplets need to be tuned to achieve a
flat longitudinal field. A triplet module where the coupler cell is a
horizontal coupler cell, previously shown to considerably reduce the reflected
wakepotential, combines a reduction of feed waveguides and effective
mitigation of reflected wakefields.

Figure 10: Triplet module based on the long slot cell with half symmetry.
Figure 11: Reflected transverse wakepotentials of structures composed of
singlets (long slot cell), doublet (one long slot and one four HOM waveguide
cells) and triplet (one long slot and two four HOM waveguide cells).

##  V Conclusion

A SW cell applicable to distributed coupling structures was developed. The HOM
waveguide width, w, was effective at lowering the wakepotential at the second
bunch position. The wakepotentials were found to be generally less sensitive
to changes in the HOM waveguide aperture, iw. From the SW cell, two
configurations were developed capable of supporting
$100\text{\,}\mathrm{M}\mathrm{V}\mathrm{/}\mathrm{m}$ loaded gradient and
gave transverse wakepotentials below the threshold. The wakefields horizontal
coupler cell couple weakly to the high power waveguide and are kept close to
open boundary levels. The triplet module increases the number of HOM
waveguides and decreases the number of coupler waveguides from which
reflections originate to keep the reflected wakepotential low.

## References

  * Aicheler  _et al._ [2012] M. Aicheler, P. Burrows, M. Draper, and T. Garvey,  _A Multi-TeV Linear Collider Based on CLIC Technology_, Tech. Rep. (2012).
  * Zha and Grudiev [2016] H. Zha and A. Grudiev, Design and optimization of Compact Linear Collider main linac accelerating structure, Physical Review Accelerators and Beams 19, 1 (2016).
  * Khan [2013] V. Khan, A Standing Wave Structure Possibility For CLIC Main LINACs, in  _International Workshop on Breakdown Science and High Gradient Technology HG2013, Trieste, Italy_ (2013).
  * Tantawi  _et al._ [2020] S. Tantawi, M. Nasr, Z. Li, C. Limborg, and P. Borchard, Design and demonstration of a distributed-coupling linear accelerator structure, Physical Review Accelerators and Beams 23, 92001 (2020).
  * Jiang  _et al._ [2021] Y. Jiang, J. Shi, H. Zha, J. Liu, X. Lin, and H. Chen, Analysis and design of parallel-coupled high-gradient structure for ultrashort input power pulses, Physical Review Accelerators and Beams 24, 112002 (2021).
  * Dolgashev and Tantawi [2010] V. A. Dolgashev and S. G. Tantawi, Study of Basic RF Breakdown Phenomena in Hight Gradient Vacuum Structures, in  _Proceedings of Linear Accelerator Conference LINAC2010, Tsukuba, Japan_ (2010) pp. 1043–1047.
  * Simakov  _et al._ [2018] E. I. Simakov, V. A. Dolgashev, and S. G. Tantawi, Advances in high gradient normal conducting accelerator structures, Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 907, 221 (2018).
  * Bane [2003] K. L. F. Bane, Short range dipole wakefields in accelerating structures for the NLC,  94309, 1 (2003).
  * Rimmer [1993] R. A. Rimmer, Extraction and Absorption of Higher Order Modes in Room Temperature Accelerators, in  _Worskhop on Microwave-Absorbing Materials for Accelerators, Newport News, VA, February 22-24, 1993_ (1993).
  * Bane  _et al._ [2018] K. L. Bane, T. L. Barklow, M. Breidenbach, C. P. Burkhart, E. A. Fauve, A. R. Gold, V. Heloin, Z. Li, E. A. Nanni, M. Nasr, M. Oriunno, J. M. Paterson, M. E. Peskin, T. O. Raubenheimer, and S. G. Tantawi,  _An Advanced NCRF Linac Concept for a High Energy e +e- Linear Collider_, Tech. Rep. (2018) arXiv:1807.10195 .
  * Grudiev and Wuensch [2010] A. Grudiev and W. Wuensch, Design of the CLIC Main Linac Accelerating Structure for CLIC Conceptual Design Report, in  _Proceedings of Linear Accelerator Conference LINAC2010_ (2010).
  * [12] Computer Simulation Technology, CST Studio Suite.
  * Schulte [2010] D. Schulte, Multi-Bunch Calculations in the CLIC Main Linac, in  _Proceedings, 23rd Conference, PAC’09, Vancouver, Canada, May 4-8, 2009_ (2010).