# Possible explanation of the irregular energy dependence of the rapidity
width of $\phi$ mesons observed in Pb+Pb collisions

Tom Reichert  Institut für Theoretische Physik, Goethe-Universität Frankfurt,
Max-von-Laue-Str. 1, D-60438 Frankfurt am Main, Germany  Frankfurt Institute
for Advanced Studies (FIAS), Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main,
Germany  Helmholtz Research Academy Hesse for FAIR (HFHF), GSI
Helmholtzzentrum für Schwerionenforschung GmbH, Campus Frankfurt, Max-von-
Laue-Str. 12, 60438 Frankfurt am Main, Germany  Jan Steinheimer  Frankfurt
Institute for Advanced Studies (FIAS), Ruth-Moufang-Str. 1, D-60438 Frankfurt
am Main, Germany  Marcus Bleicher  Institut für Theoretische Physik, Goethe-
Universität Frankfurt, Max-von-Laue-Str. 1, D-60438 Frankfurt am Main, Germany
GSI Helmholtzzentrum für Schwerionenforschung GmbH, Planckstr. 1, D-64291
Darmstadt, Germany  Helmholtz Research Academy Hesse for FAIR (HFHF), GSI
Helmholtzzentrum für Schwerionenforschung GmbH, Campus Frankfurt, Max-von-
Laue-Str. 12, 60438 Frankfurt am Main, Germany

(August 28, 2024)

###### Abstract

Experimental data from the NA49 collaboration show an unexpectedly steep rise
of the rapidity width of the $\phi$ meson as function of beam energy, which
was suggested as possible interesting signal for novel physics. In this work
we show that the Ultra-relativistic Quantum-Molecular-Dynamics (UrQMD) model
is able to reproduce the shapes of the rapidity distributions of most measured
hadrons and predicts a common linear increase of the width for all hadrons.
Only when following the exact same analysis technique and experimental
acceptance of the NA49 and NA61/SHINE collaborations, we find that the
extracted value of the rapidity width of the $\phi$ increases drastically for
the highest beam energy. We conclude that the observed steep increase of the
$\phi$ rapidity width is a problem of limited detector acceptance and the
simplified Gaussian fit approximation.

##  I Introduction

A promising observable for the understanding of strongly interacting matter
created in heavy ion collisions is the production of $\phi$ mesons. Being a
$\bar{s}s$ state with zero net-strangeness it has a small hadronic cross
section and may therefore transport information directly from hadronisation.
Being at the same time a multi-strange hadron and net strangeness neutral it
is further sensitive to different realizations of strangeness suppression in
peripheral collisions, either via canonical strangeness suppression or via a
$\gamma_{s}$ strangeness suppression factor. It has been further found that
many transport simulations have substantial problems to describe the yields of
$\phi$ mesons, which may hint to hitherto unknown production channels. For the
production of $\phi$ mesons in hadronic and partonic matter many mechanisms
have been proposed, such as thermal production Vovchenko:2015idt , OZI
suppressed reactions Chung:1997mp , kaon coalescence Ko:1991kw , catalytic
reactions Kolomeitsev:2009yn ; Tomasik:2011aa  or from heavy resonance decays
Steinheimer:2015sha .

The production of $\phi$ mesons has been studied extensively in experiments at
the CERN-LHC ALICE:2014jbq , RHIC PHENIX:2004spo ; STAR:2008inc ; STAR:2008bgi
, CERN-SPS NA49:2000jee ; NA49:2008goy ; NA61SHINE:2019gqe , BNL-AGS
E917:2003gec  and at GSI FOPI:2002csf ; HADES:2009lnd ; HADES:2017jgz .
Especially deep sub-threshold production of strangeness in the $\phi$ and
$\Xi$ hadrons has been subject of low energy and large baryon density
experiments e.g. at HADES HADES:2009lnd ; HADES:2009mtu ; HADES:2017jgz
accompanied by theoretical efforts Tomasik:2011aa ; Steinheimer:2015sha .

The NA49 collaboration at the CERN-SPS has reported on the rapidity widths of
several hadrons ranging from $\pi$, $K^{\pm}$ to $\bar{\Lambda}$ and to the
$\phi$ meson for central Pb+Pb collisions in the energy range from
$E_{\mathrm{lab}}=20$AGeV to $E_{\mathrm{lab}}=158$AGeV NA49:2008goy .
Recently, additional measurements in p+p collisions in the same energy range
have been reported by the NA61/SHINE experiment NA61SHINE:2019gqe . The
general observations for $\pi$, $K^{\pm}$ and $\bar{\Lambda}$ was that the
rapidity widths increase linearly as a function of the beam rapidity both in
Pb+Pb and in p+p collisions in the investigated energy regime and that the
rapidity width in Pb+Pb and p+p are very similar to each other. These results
are summarized as open (p+p) and filled (Pb+Pb) symbols in figure 1. The
$\phi$ seems to follow this linear trend in p+p collisions, but shows a stark
difference in Pb+Pb collisions. For Pb+Pb collisions the data suggest a
dramatic deviation from the trend observed for all other hadrons, namely a
much broader rapidity width of the $\phi$ meson for the highest collision
energy.

In this article we aim to explore the reason for this apparent broadening of
the rapidity width of the $\phi$ in Pb+Pb collisions at 158A GeV using
microscopic transport simulations. We will first calculate the rapidity
distributions and widths directly from the UrQMD simulations and then proceed
with the analysis following the experimental procedure to extract $\sigma_{y}$
of the $\phi$ meson.

##  II The UrQMD model and resonance reconstruction

The Ultra-relativistic Quantum Molecular Dynamics model (UrQMD) Bass:1998ca ;
Bleicher:1999xi ; Bleicher:2022kcu  is based on the covariant propagation of
hadrons and their resonances in phase space. As a QMD-type model it keeps
track of all n-body correlations providing a realistic time evolution for the
n-body phase space distribution function. In its current version UrQMD
includes mesonic and baryonic hadrons and resonances up to 4 GeV in mass.
Their binary interactions are modeled via measured or derived cross sections,
while additional QMD potentials can be switched on but will not be used in the
present work.

Experimentally, unstable resonances like the $\phi$ are reconstructed from the
invariant masses of their decay products. This means, from the invariant mass
distributions of the measured $\bar{K}K$-pairs, the yield of the $\phi$ is
obtained by fitting a Breit-Wigner distribution in the mass range of the
$m_{\phi}$ after subtracting the uncorrelated background. Thus, only $\phi$
mesons where both decay daughter particles escape without interaction from the
fireball can be reconstructed experimentally. Therefore, in the simulation,
measurable resonances are identified by tracing their daughter particles and
counting only those $\phi$ mesons for which both daughter particles escape
without further interactions. This method is well established and has been
extensively tested at various collision energies and resonances
Bleicher:2002dm ; Bleicher:2002rx ; Vogel:2005pd ; Knospe:2015nva ;
Reichert:2019lny ; Knospe:2021jgt ; Reichert:2022uha .

Figure 1: [Color online] Comparison between experimental results NA49:2000jee
; NA61SHINE:2019gqe ; NA61SHINE:2013tiv ; Pulawski:thesis ; NA49:2008goy ;
NA49:2002pzu ; NA49:2007stj ; NA49:2003bok  of the rapidity widths of $\pi$
(blue), $K^{+}$ (red), $K^{-}$ (magenta) and
$\overline{\Lambda}+\overline{\Sigma}^{0}$ (green) in p+p (open symbols) and
Pb+Pb (filled symbols) collisions and UrQMD simulations for Pb+Pb reactions
(solid lines) in the energy range from $E_{\mathrm{lab}}=20$AGeV to
$E_{\mathrm{lab}}=158$AGeV.

##  III Results

To establish the baseline for the investigation, we first compare the rapidity
widths of $\pi$, $K^{\pm}$ and $\bar{\Lambda}$ in Pb+Pb reactions with those
presented in the NA61/SHINE article NA61SHINE:2019gqe . For this purpose we
calculate the rapidity distributions of $\pi$, $K^{+}$, $K^{-}$ and
$\overline{\Lambda}+\overline{\Sigma}^{0}$ and extract the width of the
rapidity distribution $\sigma_{y}$ from the simulated data precisely by
$\sigma_{y}^{2}=\langle y^{2}\rangle-\langle y\rangle^{2}$. Fig. 1 shows the
comparison of the rapidity widths calculated from UrQMD (solid lines) to the
experimental data NA49:2000jee ; NA61SHINE:2019gqe ; NA61SHINE:2013tiv ;
Pulawski:thesis ; NA49:2008goy ; NA49:2002pzu ; NA49:2007stj ; NA49:2003bok
of $\pi$ (blue), $K^{+}$ (red), $K^{-}$ (magenta) and
$\overline{\Lambda}+\overline{\Sigma}^{0}$ (green) in Pb+Pb (filled symbols)
collisions. The experimental data in p+p collisions are also shown as open
symbols.

Figure 2: [Color online] Rapidity distributions of reconstructable $\phi$
mesons in 0-5% central Pb+Pb collisions at E${}_{\mathrm{lab}}=158A$ GeV (red,
scaled by 2.1) from UrQMD in cascade mode. The symbols represent experimental
data taken at $E_{\mathrm{lab}}=158A$ GeV ($\sqrt{s_{\mathrm{NN}}}=17.3$ GeV)
by NA49 NA49:2008goy  (black circles), earlier by NA49 NA49:2000jee  (grey
squares) and by CERES CERES:2005wwe  (blue hexagons). The reflected
experimental data is shown by brighter symbols and colors.

One observes that the transport model calculations of $\sigma_{y}$ agree very
well with the experimentally measured data for all investigated hadron
species. Both model and data show a linear increase of the width of the
rapidity distributions with increasing beam rapidity, with a roughly constant
slope across all hadron species.

As reported in NA49:2008goy  the $\sigma_{y}$ of the $\phi$ seems to not
follow this behavior, especially close to $\sqrt{s_{\mathrm{NN}}}\approx 20$
GeV.

Before calculating the width of the rapidity distribution of the $\phi$, the
rapidity density should be compared directly between the model calculation and
the experimental data to confirm the ability of the model to describe the
shape of the rapidity distribution. Figure 2 shows the rapidity distribution
of reconstructable $\phi$ mesons in 0-5% central Pb+Pb collisions at 158 AGeV
(red, scaled by a factor 2.1) from UrQMD. The symbols represent experimental
data by the NA49 collaboration NA49:2008goy  (black circles), earlier data by
the NA49 NA49:2000jee  (grey squares) and by the CERES experiment
CERES:2005wwe  (blue hexagons). The reflected experimental data is shown by
lighter/open symbols.

We note that the UrQMD simulations have been scaled for better comparison by a
factor of 2.1 (obtained through a fit to the available data) such that the
integrated yield of the $\phi$ matches the experimental value. This scaling
does, however, not affect the shape and width of the calculated rapidity
distribution.

The shape of the rapidity density of $\phi$ mesons calculated with UrQMD
describes the measured rapidity density very well, especially in terms of the
width. This is also quantified in the reduced chi-squared
$\tilde{\chi}^{2}\equiv\chi^{2}/d.o.f.$ of the UrQMD-Data comparison which
evaluates to $\tilde{\chi}_{\mathrm{UrQMD}}^{2}=0.16$ for UrQMD and to
$\tilde{\chi}_{\mathrm{Gauss}}^{2}=0.10$ for a single Gaussian fit to the
data. The quality of both fits is hence very good and both can be considered a
good description of the available data.111Although the overall normalization
of the UrQMD curve underestimates the measured $\phi$ yields, the
$\chi^{2}/d.o.f.$ of the scaled curve suggests that the kinematics of the
$\phi$ production are handled properly, hence allowing to reliably calculate
the width of the rapidity distribution.

Figure 3: [Color online] The rapidity width $\sigma_{y}$ as a function of the
beam rapidity of the $\phi$ calculated using the experimental technique (red
filled circles) and calculated directly from the data set (red triangles-up)
in 0-5% central Pb+Pb collisions at different energies from UrQMD. The UrQMD
calculations of the width of the $\pi$ (blue triangles-down) are shown as a
comparison. Experimental results for the $\phi$ NA49:2000jee ; NA49:2008goy
in Pb+Pb are shown as black squares with error bars.

To calculate the width of the $\phi$ rapidity density $\sigma_{y}^{\phi}$ two
different methods can be employed:

  1. 1.

The width is directly calculated, assuming full rapidity acceptance, as
$\sigma_{y}=\sqrt{\langle(\delta y)^{2}\rangle}=\sqrt{\langle
y^{2}\rangle-\langle y\rangle^{2}}$ from the simulated data set as it was done
for the comparison of the $\pi$, $K^{\pm}$ and $\overline{\Lambda}$ widths in
figure 1.

  2. 2.

We amploy the method used by the NA49 and NA61/SHINE collaborations, i.e. we
fit the rapidity density with a Gaussian (1) and a double Gaussian222Note,
that in case of a double Gaussian the width of the distribution is not equal
to the parameter in the denominator, hence we renamed it as $\varpi$ to avoid
confusion. (2) in the acceptance of the NA49 and NA61/SHINE detectors, i.e.
$0\leq y\leq 1.5$. Following the experimental analysis, we then extract a
width from each fit and take the average width from both fits as the result,
i.e. $\sigma_{\mathrm{Exp.}}=0.5\,(\sigma_{\mathrm{Fit\,Gauss}}^{0\leq y\leq
1.5}+\sigma_{\mathrm{Fit\,double\,Gauss}}^{0\leq y\leq 1.5})$

Figure 4: [Color online] The rapidity width $\sigma_{y}$ of $\phi$ meson
calculated as the average width of the fit of a single Gaussian and a double
Gaussian in the rapidity range $0\leq y\leq y_{\mathrm{max}}$ to the $\phi$
rapidity density in 0-5% central Pb+Pb collisions at E${}_{\mathrm{lab}}=158A$
GeV from UrQMD. The experimental results NA49:2000jee ; NA49:2008goy  in Pb+Pb
are shown as a black line with error band.

$\displaystyle\frac{\mathrm{d}N}{\mathrm{d}y}$
$\displaystyle\propto\exp\left[\frac{y^{2}}{2\sigma^{2}}\right]$ (1)
$\displaystyle\frac{\mathrm{d}N}{\mathrm{d}y}$
$\displaystyle\propto\exp\left[\frac{(y+a)^{2}}{2\varpi^{2}}\right]+\exp\left[\frac{(y-a)^{2}}{2\varpi^{2}}\right]$
(2)

Fig. 3 compares the calculated rapidity widths $\sigma_{y}$ (red triangles-up)
and $\sigma_{\mathrm{Exp.}}$ (red filled circles) of the $\phi$ from UrQMD, in
0-5% central Pb+Pb collisions at different energies, as a function of the beam
rapidity. The UrQMD calculations of the width of the $\pi$ (blue triangles-
down) are shown as a comparison. Experimental results for the $\phi$
NA49:2000jee ; NA49:2008goy  in Pb+Pb are shown as black squares with error
bars.

Though the experimental value of the width $\sigma_{y}$ of the $\phi$ is still
larger than the widths of the $\phi$, calculated from the simulated UrQMD
events, one notices that the width of the rapidity density of the $\phi$
calculated by following exactly the experimental technique shows also a much
larger value than the direct calculation of $\sigma_{y}$. One should also keep
in mind that this large variation of the width occurs for three different, but
statistically very likely (i.e. all have low $\chi^{2}$), descriptions of the
data. The UrQMD model calculation of $\sigma_{\mathrm{Exp.}}$ of the $\phi$
meson employing the experimental technique comes very close to the
experimental value though the true value of the width $\sigma_{y}$ shows the
same linear trend with the same slope as the other hadrons.

To better understand the systematics of the strong increase of
$\sigma_{\mathrm{Exp.}}$, Fig. 4 shows the rapidity width
$\sigma_{\mathrm{Exp.}}$ of $\phi$ mesons as a function of the reduced
rapidity range $0\leq y\leq y_{\mathrm{max}}$, which is used as input for the
Gaussian fit, from UrQMD as red circles. The experimental results NA49:2000jee
; NA49:2008goy  in Pb+Pb are shown as a black line with error band.

The rapidity width $\sigma_{\mathrm{Exp.}}$ from the UrQMD calculation shows a
strong dependence on the considered fit range (i.e. the acceptance of the
experimental detector). Smaller acceptance windows artificially increase the
extracted rapidity width. This can be attributed to the extrapolation of the
employed fit function to large forward and backward rapidities, which deviates
from the actual tails of the simulated distributions, simply because the
distribution does not exactly resemble a double-Gaussian. This phenomenon is
not observed at lower collision energies because the acceptance in beam
rapidity allowed the single or double Gaussian fit to capture the whole
distribution, constraining it also at higher forward and backward rapidities.

Again, we want to emphasize that all fits will provide a low
$\tilde{\chi}^{2}$ as the fit usually does not put much emphasis on the tails
of the distributions. In conclusion it is not clear whether the systematic
error, given by the experiment, is sufficient to capture the uncertainties
arising from fitting a double-Gaussian in an incomplete acceptance. Therefore
also the interpretation of the steep increase of the width as a signal for new
physics is questionable.

##  IV Conclusion

In the present article we have employed the Ultra-relativistic Quantum
Molecular Dynamics (UrQMD) model to calculate 0-5% central Pb+Pb collisions
and extract the rapidity distributions of reconstructable $\phi$ mesons. The
NA49 and NA61 collaborations observed a strong increase of the width of the
rapidity density of $\phi$ mesons around $\sqrt{s_{\mathrm{NN}}}\approx 20$
GeV. We have demonstrated that this increase in $\sigma_{y}$ may arise as an
artifact from the limited rapidity coverage of the detector being fitted with
a function which may not be able to capture the tails of the distributions
correctly. The rapidity width calculated directly from the simulated data set
of the $\phi$ follows the same linear trend as seen for the $\pi$, $K^{\pm}$
and $\overline{\Lambda}+\overline{\Sigma}^{0}$ while the width extracted
within the experimental acceptance and employing the same fit function shows a
very different systematic. The interpretation of ’novel’ physics is disfavored
by this investigation and caution is advised.

###### Acknowledgements.

T.R. acknowledges support through the Main-Campus-Doctus fellowship provided
by the Stiftung Polytechnische Gesellschaft Frankfurt am Main (SPTG). The
authors acknowledge for the support the European Union’s Horizon 2020 research
and innovation program under grant agreement No 824093 (STRONG-2020). The
computational resources for this project were provided by the Center for
Scientific Computing of the GU Frankfurt and the Goethe-HLR.

## References

  * (1) V. Vovchenko, V. V. Begun and M. I. Gorenstein, Phys. Rev. C 93, no.6, 064906 (2016) doi:10.1103/PhysRevC.93.064906 [arXiv:1512.08025 [nucl-th]].
  * (2) W. S. Chung, G. Q. Li and C. M. Ko, Nucl. Phys. A 625, 347-371 (1997) doi:10.1016/S0375-9474(97)00198-X [arXiv:nucl-th/9704002 [nucl-th]].
  * (3) C. M. Ko and B. H. Sa, Phys. Lett. B 258, 6-10 (1991) doi:10.1016/0370-2693(91)91199-6
  * (4) E. E. Kolomeitsev and B. Tomasik, J. Phys. G 36, 095104 (2009) doi:10.1088/0954-3899/36/9/095104 [arXiv:0903.4322 [nucl-th]].
  * (5) B. Tomasik and E. E. Kolomeitsev, Acta Phys. Polon. Supp. 5, 201-206 (2012) doi:10.5506/APhysPolBSupp.5.201 [arXiv:1112.1437 [nucl-th]].
  * (6) J. Steinheimer and M. Bleicher, J. Phys. G 43, no.1, 015104 (2016) doi:10.1088/0954-3899/43/1/015104 [arXiv:1503.07305 [nucl-th]].
  * (7) B. B. Abelev et al. [ALICE], Phys. Rev. C 91, 024609 (2015) doi:10.1103/PhysRevC.91.024609 [arXiv:1404.0495 [nucl-ex]].
  * (8) S. S. Adler et al. [PHENIX], Phys. Rev. C 72, 014903 (2005) doi:10.1103/PhysRevC.72.014903 [arXiv:nucl-ex/0410012 [nucl-ex]].
  * (9) B. I. Abelev et al. [STAR], Phys. Lett. B 673, 183-191 (2009) doi:10.1016/j.physletb.2009.02.037 [arXiv:0810.4979 [nucl-ex]].
  * (10) B. I. Abelev et al. [STAR], Phys. Rev. C 79, 064903 (2009) doi:10.1103/PhysRevC.79.064903 [arXiv:0809.4737 [nucl-ex]].
  * (11) S. V. Afanasiev et al. [NA49], Phys. Lett. B 491, 59-66 (2000) doi:10.1016/S0370-2693(00)01023-6
  * (12) C. Alt et al. [NA49], Phys. Rev. C 78, 044907 (2008) doi:10.1103/PhysRevC.78.044907 [arXiv:0806.1937 [nucl-ex]].
  * (13) A. Aduszkiewicz et al. [NA61/SHINE], Eur. Phys. J. C 80, no.3, 199 (2020) doi:10.1140/epjc/s10052-020-7675-6 [arXiv:1908.04601 [nucl-ex]].
  * (14) B. B. Back et al. [E917], Phys. Rev. C 69, 054901 (2004) doi:10.1103/PhysRevC.69.054901 [arXiv:nucl-ex/0304017 [nucl-ex]].
  * (15) A. Mangiarotti et al. [FOPI], Nucl. Phys. A 714, 89-123 (2003) doi:10.1016/S0375-9474(02)01366-0 [arXiv:nucl-ex/0209012 [nucl-ex]].
  * (16) G. Agakishiev et al. [HADES], Phys. Rev. C 80, 025209 (2009) doi:10.1103/PhysRevC.80.025209 [arXiv:0902.3487 [nucl-ex]].
  * (17) J. Adamczewski-Musch et al. [HADES], Phys. Lett. B 778, 403-407 (2018) doi:10.1016/j.physletb.2018.01.048 [arXiv:1703.08418 [nucl-ex]].
  * (18) G. Agakishiev et al. [HADES], Phys. Rev. Lett. 103, 132301 (2009) doi:10.1103/PhysRevLett.103.132301 [arXiv:0907.3582 [nucl-ex]].
  * (19) S. A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravina, C. Ernst, L. Gerland, M. Hofmann, S. Hofmann and J. Konopka, et al. Prog. Part. Nucl. Phys. 41, 255-369 (1998) doi:10.1016/S0146-6410(98)00058-1 [arXiv:nucl-th/9803035 [nucl-th]].
  * (20) M. Bleicher, E. Zabrodin, C. Spieles, S. A. Bass, C. Ernst, S. Soff, L. Bravina, M. Belkacem, H. Weber and H. Stoecker, et al. J. Phys. G 25, 1859-1896 (1999) doi:10.1088/0954-3899/25/9/308 [arXiv:hep-ph/9909407 [hep-ph]].
  * (21) M. Bleicher and E. Bratkovskaya, Prog. Part. Nucl. Phys. 122, 103920 (2022) doi:10.1016/j.ppnp.2021.103920
  * (22) M. Bleicher and J. Aichelin, Phys. Lett. B 530, 81-87 (2002) doi:10.1016/S0370-2693(02)01334-5 [arXiv:hep-ph/0201123 [hep-ph]].
  * (23) M. Bleicher, Nucl. Phys. A 715, 85-94 (2003) doi:10.1016/S0375-9474(02)01416-1 [arXiv:hep-ph/0212378 [hep-ph]].
  * (24) S. Vogel and M. Bleicher, Phys. Rev. C 74, 014902 (2006) doi:10.1103/PhysRevC.74.014902 [arXiv:nucl-th/0509105 [nucl-th]].
  * (25) A. G. Knospe, C. Markert, K. Werner, J. Steinheimer and M. Bleicher, Phys. Rev. C 93, no.1, 014911 (2016) doi:10.1103/PhysRevC.93.014911 [arXiv:1509.07895 [nucl-th]].
  * (26) T. Reichert, P. Hillmann, A. Limphirat, C. Herold and M. Bleicher, J. Phys. G 46, no.10, 105107 (2019) doi:10.1088/1361-6471/ab34fa [arXiv:1903.12032 [nucl-th]].
  * (27) A. G. Knospe, C. Markert, K. Werner, J. Steinheimer and M. Bleicher, Phys. Rev. C 104, no.5, 054907 (2021) doi:10.1103/PhysRevC.104.054907 [arXiv:2102.06797 [nucl-th]].
  * (28) T. Reichert and M. Bleicher, Nucl. Phys. A 1028, 122544 (2022) doi:10.1016/j.nuclphysa.2022.122544 [arXiv:2206.00410 [nucl-th]].
  * (29) N. Abgrall et al. [NA61/SHINE], Eur. Phys. J. C 74, no.3, 2794 (2014) doi:10.1140/epjc/s10052-014-2794-6 [arXiv:1310.2417 [hep-ex]].
  * (30) S. Pulawski, Ph.D. thesis (2015), University of Silesia, Katowice, https://edms.cern.ch/ui/file/1550661/1/main.pdf.
  * (31) S. V. Afanasiev et al. [NA49], Phys. Rev. C 66, 054902 (2002) doi:10.1103/PhysRevC.66.054902 [arXiv:nucl-ex/0205002 [nucl-ex]].
  * (32) C. Alt et al. [NA49], Phys. Rev. C 77, 024903 (2008) doi:10.1103/PhysRevC.77.024903 [arXiv:0710.0118 [nucl-ex]].
  * (33) T. Anticic et al. [NA49], Phys. Rev. Lett. 93, 022302 (2004) doi:10.1103/PhysRevLett.93.022302 [arXiv:nucl-ex/0311024 [nucl-ex]].
  * (34) D. Adamova et al. [CERES], Phys. Rev. Lett. 96, 152301 (2006) doi:10.1103/PhysRevLett.96.152301 [arXiv:nucl-ex/0512007 [nucl-ex]].


[a]Walter Wilcox

# Multipolynomial Monte Carlo Trace Estimation

Paul Lashomb  Ronald B. Morgan  Travis Whyte

###### Abstract

In lattice QCD the calculation of disconnected quark loops from the trace of
the inverse quark matrix has large noise variance. A multilevel Monte Carlo
method is proposed for this problem that uses different degree polynomials on
a multilevel system. The polynomials are developed from the GMRES algorithm
for solving linear equations. To reduce orthogonalization expense, the highest
degree polynomial is a composite or double polynomial found with a polynomial
preconditioned GMRES iteration. Matrix deflation is used in three different
ways: in the Monte Carlo levels, in the main solves, and in the deflation of
the highest level double polynomial. A numerical comparison with optimized
Hutchinson is performed on a quenched $24^{4}$ lattice. The results
demonstrate that the new Multipolynomial Monte Carlo method can significantly
improve the trace computation for matrices that have a difficult spectrum due
to small eigenvalues.

##  1 Introduction

Disconnected quark loop operators on the lattice come in form of a trace,

$\langle\bar{\psi}\Theta\psi\rangle=-Tr(\Theta M^{-1}).$ (1)

This is a challenging problem that is generally solved with Monte Carlo or
noise sampling. We give a new method that uses multilevel Monte Carlo with
different degree polynomials used to form the levels. It also incorporates
eigenvalue deflation in various crucial ways. We show that this approach can
significantly improve computation compared to simple Monte Carlo in which
$Tr(\Theta M^{-1})$ is approximated with a Monte Carlo method using $N$
Hutchinson samples: $\frac{1}{N}\sum_{i=1}^{N}\eta^{\dagger}_{i}\Theta
M^{-1}\eta_{i}$ [1]. We will use $\eta_{i}$ as a vector of $Z(4)$ noise, but
other choices are possible. A large system of linear equations
$Mx_{i}=\eta_{i}$ must be solved for each sampling, and this is the main
expense of the method.

##  2 Subtraction Methods

Subtraction methods are needed in order to reduce the variance of noisy
calculations. The noise subtracted trace estimator is

$\displaystyle Tr(\Theta M^{-1})=$
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}\eta^{(n)\dagger}\Theta\left(M^{-1}-\tilde{M}^{-1}\right)\eta^{(n)}+Tr(\Theta\tilde{M}^{-1}),$
(2)

where $M^{-1}$ is the inverse quark matrix and $\tilde{M}^{-1}$ is the
subtracted matrix which reduces the variance[2]. The second term in Equation
(2) is a trace correction term which does not affect the variance. This forms
the basis for all of our subtraction methods, each of which involve a
different approximation, $\tilde{M}^{-1}$, for the inverse of the Wilson-Dirac
matrix $M^{-1}$.

##  3 The GMRES Polynomial

We use a polynomial, $\pi(M)$, formed from the GMRES algorithm[3]. It is given
by

$\pi(M)=\prod_{i=1}^{d}\Big{(}1-\frac{M}{\tilde{\theta}_{i}}\Big{)},$ (3)

where $\tilde{\theta}_{1},\tilde{\theta}_{2},\ldots,\tilde{\theta}_{d}$ are
the Leja ordered [4], harmonic Ritz values obtained from a single cycle of
GMRES($d$), and $d$ is the desired degree of the GMRES polynomial. This method
of finding the polynomial is much more stable than using the normal
equations[5], allowing for high degree polynomials to be formed.

The GMRES polynomial $p(M)$ of degree $d-1$ related to $\pi(M)$ is given by
$\pi(M)=1-M\,p(M)$. A backwards solve shows it can be expressed as

$p(M)=\sum_{k=1}^{d}u_{k},\,\,u_{k}=\frac{M}{\tilde{\theta}_{k}}\left(1-\frac{M}{\tilde{\theta}_{1}}\right)\left(1-\frac{M}{\tilde{\theta}_{2}}\right)\cdots\left(1-\frac{M}{\tilde{\theta}_{k-1}}\right),$
(4)

where $d$ is the order of the polynomial $\pi(M)$. We only need to run a
single cycle of GMRES($d$) in order to extract $d$ harmonic Ritz values. These
are then Leja ordered for numerical stability. For subtraction, we make the
assignment $\tilde{M}^{-1}\equiv p(M)$.

##  4 Relative Variance

In order to monitor the reduction in variance of each method over that of no
subtraction (NS), we define the relative variance as

$\bar{\sigma}_{R}^{2}\equiv\frac{\bar{\sigma}_{A}^{2}}{\bar{\sigma}_{NS}^{2}},$
(5)

where $\bar{\sigma}_{A}^{2}$ is the variance obtained from a particular
subtraction method $A$ averaged over the configurations, and
$\bar{\sigma}^{2}_{NS}$ is the variance of the trace estimator with no
subtraction averaged over all configurations.

The standard relative variance results are often dominated by a single
configuration and have large fluctuations. Log-averaged relative variance
results, although not unbiased, have significantly reduced errors and are far
easier for other numerical groups to make comparisons to. Thus, we also define
a base-10 log-averaged relative variance using the geometric mean [6] as

$\bar{\sigma}_{R,log}^{2}\equiv 10^{\bar{\rho}^{2}_{A}-\bar{\rho}^{2}_{NS}},$
(6)

where
$\bar{\rho}^{2}_{A}=\frac{1}{N}\sum_{j=1}^{N}\log_{10}[(\sigma_{A}^{2})^{(j)}]$
($\bar{\rho}^{2}_{NS}$ defined similarly), $(\sigma_{A}^{2})^{(j)}$ is the
variance obtained for configuration $j$ for a particular method, and $N$ is
the number of configurations.

We now turn our attention to considering the polynomial degree dependence of
our polynomial noise subtraction methods. We first examined both the standard
and log-averaged relative variance for several polynomial degrees on different
lattice volumes[7]. Here we consider the scalar operator but with four
quenched lattices of volumes $4^{3}\times 4$, $8^{3}\times 8$, $12^{3}\times
16$ and $24^{3}\times 32$ with $\beta=6.0$ configurations and hopping
parameter $\kappa=0.1570\approx\kappa_{crit}$ using a total of 10 thermalized
configurations for each lattice volume. The three smaller lattices used $N=10$
unpartitioned $Z(4)$ noise vectors for the Monte Carlo while the largest used
$N=100$. $N=100$ was sufficient to reveal the relative variance behavior for
the $24^{3}\times 32$ lattices at approximately at the same relative variance
error level as the $12^{3}\times 16$ lattices.

Figure 1: Log-averaged standard relative variance of the scalar operator using
polynomial subtraction against the degree of the subtraction polynomial.
Lattice volumes $4^{3}\times 4$, $8^{3}\times 8$, $12^{3}\times 32$, and
$24^{3}\times 32$ are shown, each averaged over 10 quenched configurations at
$\beta=6.0$ and $\kappa=0.1570$.

Figure 1 shows the log-averaged standard relative variance of the scalar
operator versus the polynomial degree for each of the four lattice volumes
that we used. Error bars are indicated by the shaded regions. Note the
approximate exponential falloff in the relative variances for larger
polynomial degrees for either figure. We also see that the increase in lattice
volume requires a higher polynomial degree to achieve the same variance
reduction performance. This polynomial degree can become impractically high
for larger lattices. To increase the polynomial degree further, one would need
to reduce the expense of computing such high degree polynomials as the single
cycle of GMRES($d+1$) involved in forming the subtraction polynomial of degree
$d$ becomes expensive due to the orthogonalization costs of the Arnoldi
iteration. As will be seen shortly, double polynomials as discussed in Refs.
[8] and [3], and described in the next section can be used to combat this.

##  5 Reducing Costs Through Double Polynomials

The degree of the subtraction polynomials are dependent on the size of the
Krylov subspace used to form them. Since the GMRES algorithm is based on the
Arnoldi iteration, forming large subspaces to produce high degree subtraction
polynomials can increase the orthogonalization and memory costs considerably.
As a cost-saving measure, we may make use of double polynomials to avoid the
large orthogonalization and memory costs of the Arnoldi iteration. We start
this section by first outlining double polynomials in the context of
polynomial preconditioning as they were originally [8] used, and then we state
how they can be used in noise subtraction.

In polynomial preconditioning, we use a polynomial preconditioner $p_{in}(M)$
on the system $M$, yielding the polynomial preconditioned system
$\phi_{in}(M)\equiv Mp_{in}(M)$ as

$\phi_{in}(M)y=b,$ (7)

where the solution is $x=p_{in}(M)y$. If we then polynomial precondition the
system a second time using the preconditioner $p_{out}(\phi_{in}(M))$, we will
have a double polynomial preconditioned system

$\phi_{out}(\phi_{in}(M))z\equiv\phi_{in}(M)p_{out}(\phi_{in}(M))z=b,$ (8)

where the solution is now $x=p_{in}(M)p_{out}(\phi_{in}(M))z$. If we now
consider

$\phi_{out}(\phi_{in}(M))=Mp_{in}(M)p_{out}(\phi_{in}(M)),$ (9)

we see that performing double polynomial preconditioning can also be thought
of as applying a single polynomial preconditioner $p_{double}(M)\equiv
p_{in}(M)p_{out}(\phi_{in}(M))$ to the original system $M$.

To form such a double polynomial, we perform a single cycle of
GMRES($d_{in}+1$) on $M$ to obtain $d_{in}+1$ Leja ordered, harmonic Ritz
values for $M$, and then perform a single cycle of GMRES($d_{out}+1$) on
$\phi_{in}(M)$ to obtain $d_{out}+1$ Leja ordered, harmonic Ritz values for
$\phi_{in}(M)$. Note that the latter single cycle of GMRES($d_{out}+1$) on
$\phi_{in}(M)$ can alternatively be thought of as performing a single cycle of
polynomial preconditioned GMRES or PP-GMRES($d_{out}+1$) on $M$ to form the
$d_{out}+1$ Leja ordered, harmonic Ritz values. From these, one can form
polynomials $p_{in}(M)$ and $p_{out}(\phi_{in}(M))$ of degrees $d_{in}$ and
$d_{out}$, respectively, using the algorithm discussed in Ref. [8] to obtain a
double polynomial of degree $d_{double}=(d_{in}+1)*(d_{out}+1)-1$.

In this way, we may form two low-degree polynomials $p_{in}(M)$ and
$p_{out}(\phi_{in}(M))$ using smaller Krylov subspaces to achieve an effective
high degree polynomial, thus avoiding the large Krylov subspaces formed by the
Arnoldi iteration of GMRES applied solely on $M$. In the context of the
multipolynomial method described in the next section, we make use of this
double polynomial for subtraction purposes, $\tilde{M}^{-1}_{poly}\equiv
p_{double}(M)$, at the highest level in our multilevel Monte Carlo.

##  6 Multipolynomial Monte Carlo

Multilevel Monte Carlo [9, 10] uses different levels of approximation to the
problem being solved in order to reduce the variance. It creates several
independent Monte Carlo estimations with each except the last one having
differences between levels. The last has a Monte Carlo for the least accurate
approximation to the original problem. Essentially, much of the Monte Carlo
sampling is shifted from the original problem to sampling with the
approximations that are cheaper to use. Multilevel Monte Carlo has been
applied in QCD [11, 12]. Also, Ref. [13] has levels with different size
multigrid matrices.

In [7] it is suggested to use a high degree polynomial of $M$ to approximate
$M^{-1}$, and in Ref. [14] to turn this into a multilevel Monte Carlo method.
Consider

$Tr(M^{-1})=Tr(M^{-1}-p_{1}(M))+Tr(p_{1}(M)).$ (10)

The two traces on the right can be computed independently with Monte Carlo
sampling. Or if $p_{1}$ is low degree, the second trace can be computed
exactly with probing [15, 16]. If $p_{1}$ is a good approximation to $M^{-1}$,
the variance is greatly reduced for the $Tr(M^{-1}-p_{1}(M))$ part of the
Monte Carlo. Thus very few samples are needed for this part of the trace
estimate. Each sample is expensive however, needing solution of linear
equations and application of a high degree polynomial. We use the stable high
degree polynomials that are generated by the GMRES algorithm in [8, 3].
However, Equation (10) cannot easily be used with a high degree polynomial
because finding $Tr(p_{1}(M))$ would require a Monte Carlo approach that has
variance similar to that of the original problem of $Tr(M^{-1})$. It does gain
in that implementing $p_{1}(M)$ times a vector is generally less expensive
than solving the linear equations needed to multiply $M^{-1}$ times a vector.

To reduce the Monte Carlo cost for $Tr(p_{1}(M))$ we suggest approximating
$p_{1}$ with a lower degree polynomial, say $p_{2}$. This is done by making
$p_{2}(M)$ another approximation to $M^{-1}$ and thus also an approximation to
$p_{1}(M)$. The next equation shows how $p_{2}$ is used so that there are two
Monte Carlo’s along with the problem of finding the trace of $p_{2}(M)$:

$Tr(M^{-1})=Tr(M^{-1}-p_{1}(M))+Tr(p_{1}(M)-p_{2}(M))+Tr(p_{2}(M)).$

This is a polynomial version of multilevel Monte Carlo [17, 9, 10]. In [18],
multiple Chebyshev polynomials are used in a multilevel Monte Carlo method,
however only for symmetric matrices and not for computing the trace of an
inverse. As mentioned above, we instead use GMRES polynomials for our non-
Hermitian complex QCD matrices. Our approach is perhaps easier to set up than
a multilevel Monte Carlo method that uses multigrid matrices [13] since it
does not need the development of the matrices. However, there are also
complications for our new approach such as the need of a high degree
polynomial $p_{1}$ and the need for eigenvalue deflation.

For each lower order polynomial $p_{i}$ we want $p_{i}(M)\approx M^{-1}$.
These polynomials can be generated from a run of GMRES. For example, for a
degree 100 polynomial $p(\alpha)$, stop GMRES at iteration 101. We call the
GMRES residual polynomial $\pi(\alpha)$. It is degree 101. This polynomial in
$\alpha$ can be written as $\pi(\alpha)=1-\alpha p(\alpha)$, where $p(\alpha)$
is the polynomial that approximates the inverse of $M$. In fact, $\pi(0)=1$
and $\pi(\alpha)\approx 0$ over the spectrum of $M$, once GMRES has run far
enough. So $p(\alpha)\approx 1/\alpha$ over the spectrum and thus $p(M)\approx
M^{-1}$.

There can be any number of levels in a multilevel Monte Carlo scheme, however
there is cost for extra levels because then each individual Monte Carlo
sampling must be estimated more accurately. We will only go as high as using
three polynomials which gives four trace calculations. Let the polynomials
$p_{1}$, $p_{2}$ and $p_{3}$ have degrees be $d_{1}>d_{2}>d_{3}.$ The
multilevel Monte Carlo with polynomials for the different levels, which we
refer to as Multipolynomial Monte Carlo, uses this formula:

$Tr(M^{-1})=Tr(M^{-1}-p_{1}(M))+Tr(p_{1}(M)-p_{2}(M))+\\\
Tr(p_{2}(M)-p_{3}(M))+Tr(p_{3}(M)).$ (11)

The first three trace computations on the right of Equation (11) can be
performed with Hutchinson Monte Carlo trace estimation. Here we do the last
one exactly with probing.

##  7 Deflation for the Monte Carlo

Deflation of eigenvalues is used in three different ways in our method. There
is deflation in the Monte Carlos levels that have differences between
polynomials, deflation in the solution of linear equations, and in the
development of a deflated double polynomial. Here we discuss the Monte Carlo
level deflation. The more complete discussion is in Ref. [14].

We use polynomial preconditioned Arnoldi method (PP-Arnoldi) [8] to deflate
eigenvalues. In some of the tests, we determine the stopping point by
monitoring the relative residual of the corresponding polynomial
preconditioned GMRES. Multipolynomial Monte Carlo can reduce expense due to
moving much of the noise sampling to the cheaper Monte Carlo parts such as
$Tr(p_{2}(M)-p_{3}(M))$. However, the improvement due to having more levels
may not be dramatic due to the increased number of samples needed for more
Monte Carlos with more demanding error tolerances. Needed for a big
improvement is deflation in the Monte Carlos that have the difference between
polynomials. We find eigenvalues $\lambda_{i}$ and corresponding right and
left eigenvectors $z_{i}$ and $u_{i}$, respectively. Then the Monte Carlos for
$Tr(p_{1}(M)-p_{2}(M))$ and $Tr(p_{2}(M)-p_{3}(M))$ have deflation applied.
For the first of these, this is

$Tr\Bigl{(}p_{1}(M)-p_{2}(M)-\sum_{i}(p_{1}(\lambda_{i})-p_{2}(\lambda_{i}))z_{i}u_{i}^{*}\Bigr{)}.$

The deflation reduces the variance for Monte Carlo sampling. We correct for
the deflated part by finding its trace exactly using that the trace of
$z_{i}u_{i}^{*}$ is the inner product between the vectors.

Note that when simply doing a Monte Carlo for $Tr(M^{-1})$, deflating
eigenvalues is not very effective because of non-normality effects. For a non-
normal matrix, there can even be an increase in the norm of the inverse as
some spectral components are removed. The norm of $M^{-1}$ will not increase
if singular value decomposition components are instead removed, hence they
have been previously used [19, 21, 22, 20, 23]. Nevertheless, we remove
eigenvalue components, because the singular values and vectors of a polynomial
of $M$ are not necessarily related to those of $M$. Also, eigenvalues work
well in our context because the polynomials all approximate $M^{-1}$ but
differ mostly at the small eigenvalues that we are deflating. So removing
these components significantly reduces the variance.

##  8 Algorithm and Testing

We test with 10 quenched matrices from $24^{4}$ configurations using the
scalar operator. The desired error is set at $\epsilon=0.0005*24^{4}$. We use
two versions of the multpolynomial approach. Both are with three polynomials.
The first polynomial is from deflated PP(50)-GMRES solved to residual norm
below $10^{-5}$. This gives $p_{1}(M)$ that is a very accurate approximation
to $M^{-1}$, so the Monte Carlo for $Tr(M^{-1}-p_{1}(M))$ requires only two
noises and achieves far more than the requested accuracy. The second and third
polynomials are degree 200 and degree 4. The standard error for the second
Monte Carlo is checked after three noises and then rechecked after every
three. The third Monte Carlo is checked after six noises and then every third
thereafter. The first test computes the eigenvectors by running PP(50)-Arnoldi
until the corresponding GMRES relative residual norm goes below $10^{-12}$.
From the subspace thus generated, 30 right eigenvectors are computed
corresponding to the smallest Ritz vectors, the left eigenvectors are found
and the eigenvalue is accepted if residual norms are below $10^{-3}$. Ritz
values are identified as real if the imaginary parts are below $2*10^{-4}$. An
average of 22.4 small eigenvalues are used. The time and matrix-vector
products for computing $Tr(p_{3}(M))$ with probing are included in the table
(the time averages 0.45 hours). The time for the whole multipolynomial process
varies greatly from 3.9 hours up to 16 hours per configuration. This is
because the Monte Carlo deflation is sometimes more effective than for other
cases. The number of noise vectors needed for the second Monte Carlo varies
from 3 to 66, and every noise vector is expensive because it requires
multiplication by the high-degree polynomial $p_{1}$ (and by $p_{2}$).

Next, we check whether it is worth spending more effort developing the
eigenvalues for deflation. We use a degree 69 polynomial for $p_{in}$, and run
PP(70)-Arnoldi for 150 iterations (stopping at that set point instead of at a
GMRES residual level). This takes a total of 10,500 matrix-vector products
which compares to an average of 6250 matrix-vector products used just above
when going to $10^{-12}$. (This varies considerably for the 10 matrices,
between 4850 and 7000). Next, the 60 smallest Ritz values and vectors are
tested for accuracy. An average of 52.9 approximate eigenvalues are accepted.
Most of them are very accurate with residual norms around $10^{-15}$. This
makes the decisions more clear-cut for determining which eigenvalues to use
for deflation. We call the procedure “extra deflation". The overall time goes
down for all tests even though more time is spent on the initial steps of
finding eigenvalues and polynomials before the actual Monte Carlos.

We show a comparison of the two Multipolynomial Monte Carlo tests along with
Hutchinson in Table 1 with averaging over 10 configurations. Linear equations
are solved for Hutchinson using PP(70)-GMRES with extra deflation. The time is
reduced by a factor of more than six with the extra deflation version of the
multipolynomial method. We note that the Hutchinson method takes an average of
19.6 hours, but this would take much longer if polynomial preconditioning and
deflation were not used.

Table 1:  Trace results for $24^{4}$ lattices with error tolerance $0.0005*24^{4}$. Compare Hutchinson with the Multipolynomial Method. We average with 10 configurations. Method | Noise vectors | Time | MVP’s
---|---|---|---
Hutch., defl. PP(70)-GMRES | 188 | 19.6 hours | $2.02*10^{5}$
3 Poly’s, $10^{-12}$ for deflation | 2, 17.1, 101 | 6.84 hours | $6.68*10^{4}$
3 Poly’s, extra deflation | 2, 3, 24.9 | 3.07 hours | $2.80*10^{4}$

##  9 Conclusions

We have proposed a new approach to computing the trace of the inverse of a
large QCD matrix. It consists of a multilevel Monte Carlo method with
different degree polynomials. The polynomials give approximations of different
accuracy to the inverse of the matrix. The lower degree polynomials are from
GMRES and the highest degree polynomial is from polynomial preconditioned
GMRES.

This is the first use of polynomials for multilevel Monte Carlo with
nonsymmetric matrices. For this method to be effective, it is essential to
deflate eigenvalues from the Monte Carlo steps that have a difference of two
polynomials. This deflation is also new, as is the use of a deflated
polynomial to lower the degree of the highest degree polynomial. Also
important to the method is that this highest degree polynomial is a double
polynomial with much less orthogonalization expense needed to generate it.
Also implemented is a deflated version of polynomial preconditioned GMRES that
is very efficient and is new to QCD calculations. Computation of eigenvalues
and eigenvectors is done with a polynomial preconditioned Arnoldi method that
makes deflation much more practical. Putting the pieces together gives a
multipolynomial Monte Carlo method that is very efficient.

We plan future work on this method, including automatic optimization of the
cost function [10] as well as using nonsymmetric Lanczos to develop the
polynomials and to compute eigenvalues. This will likely give a lower degree
polynomial than the current double polynomial. For QCD matrices, this can be
implemented efficiently; see Ref. [24].

##  10 Acknowledgements

We acknowledge the use of the Baylor High Performance Cluster as well as the
Texas Advanced Computing Center (TACC). We also acknowledge the Baylor Summer
Research Award program, the previous work of Suman Baral, and the QQCD code
from Randy Lewis.

## References

  * [1] M. F. Hutchinson, Commun. Stat. Simul. Comput. 18:1059 (1989).
  * [2] S. Bernardson, P. McCarty, and C. Thron, Comput. Phys. Commun. 78, 256 (1994).
  * [3] J. A. Loe and R. B. Morgan, Numer. Linear Algebra Appl. 29: e2427 (2021).
  * [4] Z. Bai, D. Hu, and L. Reichel, IMA J. Numer. Anal. 14: 563 (1994).
  * [5] Q. Liu, W. Wilcox, and R. Morgan, arXiv: 1405.1763.
  * [6] D. J. Sheskin, Handbook of Parametric and Nonparametric Statistical Procedures, 3rd ed. (Chapman & Hall/CRC, Boca Raton, 2004), p.8.
  * [7] P. Lashomb, R. B. Morgan, T. Whyte, and W. Wilcox. arXiv: 2306.06188 (2023).
  * [8] M. Embree, J. A. Loe, and R. Morgan, SIAM J. Sci. Comput. 43: A1 (2021).
  * [9] M. B. Giles. Oper. Res., 56: 607 (2008).
  * [10] M. B. Giles. Acta Numer., 24: 259 (2015).
  * [11] L. Giusti, T. Harris, A. Nada, and S. Schaefer. Eur. Phys. J. C, 79: 586 (2019).
  * [12] T. Whyte, A. Stathopoulos, E. Romero, and K. Orginos. arXiv: 2212.04430 (2022).
  * [13] A. Frommer, M. N. Khalil, and G. Ramirez-Hidalgo. SIAM J. Sci. Comput., 44: A2536 (2022).
  * [14] P. Lashomb, R. B. Morgan, T. Whyte, and W. Wilcox. arXiv: 2307.08620 (2023).
  * [15] J. Laeuchli and A. Stathopoulos. SIAM J. Sci. Comput., 42: A1459 (2020).
  * [16] A. Stathopoulos, J. Laeuchli, and K. Orginos. SIAM J. Sci. Comput., 35: S299 (2013).
  * [17] S. Heinrich. J. Complexity, 14: 151 (1998).
  * [18] E. Hallman and D. Troester. Linear Algebra Appl., 638: 125 (2022).
  * [19] V. Guerrero, R. B. Morgan, and W. Wilcox. Proceedings of Science, LAT2009, 1 (2009).
  * [20] S. Baral, T. Whyte, W. Wilcox, and R. B. Morgan. Comput. Phys. Commun., 241: 64 (2019).
  * [21] S. Baral, W. Wilcox, and R. B. Morgan. Proceedings of Science, LAT2016, 355 (2016).
  * [22] A. S. Gambhir, A. Stathopoulos, and K. Orginos. SIAM J. Sci. Comput., 39: A532 (2017).
  * [23] E. Romero, A. Stathopoulos, and K. Orginos. J. Comput. Phys., 409: 109356 (2020).
  * [24] A. Frommer, B. Nöckel, S. Güsken, T. Lippert, and K. Schilling. Int. J. Mod. Phys., C6: 627 (1995).