# Meson spectroscopy from spectral densities in lattice gauge theories

Ed Bennett  0000-0002-1678-6701<EMAIL_ADDRESS>Swansea Academy of
Advanced Computing, Swansea University (Bay Campus), Fabian Way, Swansea SA1
8EN, United Kingdom  Luigi Del Debbio  0000-0003-4246-3305
<EMAIL_ADDRESS>Higgs Centre for Theoretical Physics, School of
Physics and Astronomy, The University of Edinburgh, Peter Guthrie Tait Road,
Edinburgh EH9 3FD, United Kingdom  Niccolò Forzano  0000-0003-0985-8858
<EMAIL_ADDRESS>Department of Physics, Faculty of Science and
Engineering, Swansea University, Singleton Park, SA2 8PP, Swansea, United
Kingdom  Ryan C. Hill  0000-0003-2383-940X<EMAIL_ADDRESS>School of
Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, United
Kingdom  Deog Ki Hong  0000-0002-3923-4184<EMAIL_ADDRESS>Department of
Physics, Pusan National University, Busan 46241, Korea  Ho Hsiao
0000-0002-8522-5190<EMAIL_ADDRESS>Institute of Physics,
National Yang Ming Chiao Tung University, 1001 Ta-Hsueh Road, Hsinchu 30010,
Taiwan  Jong-Wan Lee  0000-0002-4616-2422<EMAIL_ADDRESS>Particle Theory
and Cosmology Group, Center for Theoretical Physics of the Universe, Institute
for Basic Science (IBS), Daejeon, 34126, Korea  C.-J. David Lin
0000-0003-3743-0840<EMAIL_ADDRESS>Institute of Physics, National Yang Ming
Chiao Tung University, 1001 Ta-Hsueh Road, Hsinchu 30010, Taiwan  Centre for
High Energy Physics, Chung-Yuan Christian University, Chung-Li 32023, Taiwan
Centre for Theoretical and Computational Physics, National Yang Ming Chiao
Tung University, 1001 Ta-Hsueh Road, Hsinchu 30010, Taiwan  Biagio Lucini
0000-0001-8974-8266<EMAIL_ADDRESS>Swansea Academy of Advanced
Computing, Swansea University (Bay Campus), Fabian Way, Swansea SA1 8EN,
United Kingdom  Department of Mathematics, Faculty of Science and Engineering,
Swansea University (Bay Campus), Fabian Way, SA1 8EN Swansea, United Kingdom
Alessandro Lupo  0000-0001-9661-7811<EMAIL_ADDRESS>Aix-
Marseille Université, Université de Toulon, CNRS, CPT, Marseille, France
Maurizio Piai  0000-0002-2251-0111<EMAIL_ADDRESS>Department of Physics,
Faculty of Science and Engineering, Swansea University, Singleton Park, SA2
8PP, Swansea, United Kingdom  Davide Vadacchino  0000-0002-5783-5602
<EMAIL_ADDRESS>Centre for Mathematical Sciences, University
of Plymouth, Plymouth, PL4 8AA, United Kingdom  Fabian Zierler
0000-0002-8670-4054<EMAIL_ADDRESS>Department of Physics,
Faculty of Science and Engineering, Swansea University, Singleton Park, SA2
8PP, Swansea, United Kingdom

###### Abstract

Spectral densities encode non-perturbative information that enters the
calculation of a plethora of physical observables in strongly coupled field
theories. Phenomenological applications encompass aspects of standard-model
hadronic physics, observable at current colliders, as well as correlation
functions characterizing new physics proposals, testable in future
experiments. By making use of numerical data produced with lattice gauge
theories, we perform a systematic study to demonstrate the effectiveness of
recent technological progress in the reconstruction of spectral densities.

To this purpose, we write and test new software packages that use energy-
smeared spectral densities to analyze the mass spectrum of mesons. We assess
the effectiveness of different smearing kernels and optimize the smearing
parameters to the characteristics of available lattice ensembles. For
concreteness, we analyze the $Sp(4)$ lattice gauge theory with matter
transforming in an admixture of fundamental and 2-index antisymmetric
representations of the gauge group. We generate new ensembles for this theory,
with lattices that have a longer extent in the time direction with respect to
the spatial ones. We run our tests on these ensembles, obtaining new results
about the spectrum of light mesons and their excitations. We make available
our algorithm and software for the extraction of spectral densities, that can
be applied to theories with other gauge groups, including the theory of strong
interactions (QCD) governing hadronic physics in the standard model.

††preprint: CTPU-PTC-24-11††preprint: PNUTP-24/A02

###### Contents

  1. I Introduction
  2. II Lattice field theory
    1. II.1 Lattice discretization and bare parameters
    2. II.2 Gradient flow, topological charge and autocorrelations
  3. III Correlation functions
    1. III.1  APE and Wuppertal smearing algorithms
    2. III.2 Generalized Eigenvalue Problem
  4. IV Spectral densities
    1. IV.1 The Hansen-Lupo-Tantalo (HLT) method
    2. IV.2 Smearing kernel and smearing radius
    3. IV.3 Spectral density fits
    4. IV.4 On the effect of APE and Wuppertal smearing on the spectral densities
  5. V Spectral energy density reconstruction and lattice time extent
  6. VI Numerical results and comparisons with GEVP
  7. VII Summary and outlook
  8. A More about spectral density reconstruction

##  I Introduction

Recent years have seen the development of new technology aimed at extracting
spectral densities from numerical lattice data obtained in the non-
perturbative study of strongly coupled field theories—see for instance Refs.
Hansen:2019idp ; Hansen:2017mnd ; Bulava:2019kbi ; Bailas:2020qmv ;
Gambino:2020crt ; Bruno:2020kyl ; Lupo:2021nzv ; Gambino:2022dvu ;
DelDebbio:2022qgu ; Bulava:2021fre ; Lupo:2022nuj ; Bulava:2023brj ;
DeSantis:2023rjl ; Lupo:2023qna ; Kades:2019wtd ; Pawlowski:2022zhh ;
DelDebbio:2021whr ; Bergamaschi:2023xzx ; Bonanno:2023thi ; Frezzotti:2023nun
; Frezzotti:2024kqk . Spectral densities are inverse Laplace transforms of
space-averaged two-point functions involving time-separated operators. They
can be used to compute high-precision spectral observables that are otherwise
difficult to access with conventional methodologies. Additionally, spectral
densities encode information about off-shell physics, providing an alternative
framework to compute scattering amplitudes Bulava:2019kbi , inclusive rates
Gambino:2020crt ; Bulava:2021fre , finite-volume energies and matrix elements.
In QCD, spectral densities enter the calculation of hadronic observables, such
as the $R$-ratio ExtendedTwistedMassCollaborationETMC:2022sta  and inclusive
decays of the $\tau$ lepton Evangelista:2023fmt ; Alexandrou:2024gpl . They
can also be used to access the properties of glueballs Smecca:2023hgr . In
analogy to the derivation of the Weinberg sum rules from the properties of the
spectral functions Weinberg:1967kj ; Dash:1967fq ; Bernard:1975cd , the
spectral representation of two-point functions involving both mesons and
baryons enters the effective potential of new physics models Contino:2010rs ;
Banerjee:2023ipb , that implement Higgs compositeness Kaplan:1983fs ;
Georgi:1984af ; Dugan:1984hq  and top partial compositeness Kaplan:1991dc
(see also the reviews in Refs. Panico:2015jxa ; Witzel:2019jbe ;
Cacciapaglia:2020kgq ; Bennett:2023wjw , and the summary tables in Refs.
Ferretti:2013kya ; Ferretti:2016upr ; Cacciapaglia:2019bqz ), and may trigger
electroweak symmetry breaking via vacuum misalignment Das:1967it ;
Peskin:1980gc ; Preskill:1980mz .

This paper has two main objectives. The first is to develop, test, tune,
optimize, and benchmark the effectiveness of a new software package that
allows the computation of spectral densities from correlation functions
measured on the lattice. The spectral densities of interest are smeared in
energy, and thus have a well-defined infinite-volume limit. After such limit
is taken, the energy smearing can in principle be removed. This is however not
necessary in this work, where we perform correlated fits of smeared spectral
densities to extract the finite-volume spectrum of mesons, as proposed in Ref.
DelDebbio:2022qgu . The computation of smeared spectral densities can be
systematically improved by reducing the statistical noise, and by increasing
the number of time slices in the lattice Hansen:2019idp . To this purpose, we
generate and analyze ensembles with long time extent. This work sets the stage
for future applications, both in the context of QCD and of new physics models,
by demonstrating its viability as an analysis tool for the aforementioned
ambitious endeavors. To validate our techniques., we apply these analysis
tools to observables in theories for which alternative ways exist to gain
access to the relevant non-perturbative information.

The second objective is to make significant progress toward understanding the
spectrum of a specific gauge theory that is a prototypical candidate for UV-
completion of Composite Higgs Models Barnard:2013zea . To this purpose, we
generate new ensembles for the $Sp(4)$ theory coupled to two Dirac fermions
transforming in the fundamental and three in the two-index antisymmetric
representation of the gauge group. This theory and its variations have been
recently studied in the literature on $Sp(2N)$ theories Bennett:2017kga ;
Lee:2018ztv ; Bennett:2019jzz ; Bennett:2019cxd ; Bennett:2020hqd ;
Bennett:2020qtj ; Lucini:2021xke ; Bennett:2021mbw ; Bennett:2022yfa ;
Bennett:2022gdz ; Bennett:2022ftz ; Lee:2022elf ; Hsiao:2022kxf ; Maas:2021gbf
; Zierler:2021cfa ; Kulkarni:2022bvh ; Bennett:2023wjw ; Bennett:2023rsl ;
Bennett:2023gbe ; Bennett:2023mhh ; Bennett:2023qwx ; Pomper:2024otb . We use
the Grid software suite Boyle:2015tjk ; Boyle:2016lbp ; Yamaguchi:2022feu ,
including adaptations Bennett:2023gbe  previously made to implement symplectic
groups. The data are analyzed with the HiRep code DelDebbio:2008zf . We write
and make available new software that reconstructs the smeared energy spectral
density. We focus our analysis on the two-point correlation functions of meson
operators, with a general basis of Dirac structures, but restrict our
attention to flavor non-singlet states. We will analyze the spectra of flavor
singlet mesons and of chimera baryons in separate publications.

Our results on the spectrum of mesons show improved statistical accuracy, with
better control over systematics, with respect to existing published results,
and extend over a larger range of parameter space. We also detect new excited
states, not accessible with earlier existing ensembles. We make publicly
available both the software we developed and the data analysis flow
analysis_release , as well as the new data generated for this study
data_release . In the present paper, we defer to future studies of the
extraction of the couplings (decay constants) of the associated mesons.

The paper is organized as follows. We define the lattice theories of interest
in Sec. II. We present the properties of the ensembles we study in the main
body of the paper. We apply the gradient flow as scale-setting procedure, as
well as a smoothening process for topological observables. We monitor the
topological charge, and use it to estimate autocorrelations. In Sec. III we
introduce the flavored, mesonic operators of interest and their correlation
functions. We also describe our implementation of Wuppertal and APE smearing,
and exemplify how the spectra can be extracted from a variational approach
based upon the generalized eigenvalue problem (GEVP). Spectral densities and
(energy) smearing kernels are introduced in Sec. IV. We also discuss how the
signal and its statistical significance depend on the choices of smearing
parameters. We devote Sec. V to a systematic investigation of how the length
of the time extent of the lattice affects the spectral density reconstruction.
Our (new) results on the spectrum of mesons are summarized and critically
discussed in Sec. VI. We conclude by highlighting avenues for further
investigation in Sec. VII. Technical details about the spectral density
reconstruction are relegated to the Appendix.

##  II Lattice field theory

In this section, we present the quantum field theory of interest and the
discretized lattice action we adopt for its study. By doing so, we fix the
notation so that the presentation is self-contained. We also tabulate and
characterize the lattice ensembles generated for the purposes of this paper,
in which we report our results on the spectrum of flavored mesons. We postpone
to future publications the measurement, on these same ensembles, of the
spectra of other bound states—flavor singlet mesons and chimera baryons.

###  II.1 Lattice discretization and bare parameters

The $Sp(4)$ gauge theory of interest has a continuum Lagrangian density that,
in the presence of $N_{\rm f}$ Dirac fermions, $Q^{I}$, transforming in the
fundamental representation, together with $N_{\rm as}$ Dirac fermions,
$\Psi^{k}$, transforming in the two-index antisymmetric representation, is
given by

$\displaystyle\mathcal{L}=-\frac{1}{2}{\rm
Tr}\,G_{\mu\nu}G^{\mu\nu}+\sum_{I,J=1}^{N_{\rm
f}}\bar{Q}^{I}\left(i\delta_{IJ}\,\gamma^{\mu}D_{\mu}-m_{IJ}^{\rm
f}\right)Q^{J}+\sum_{k,\ell=1}^{N_{\rm
as}}\bar{\Psi}^{k}\left(i\delta_{k\ell}\,\gamma^{\mu}D_{\mu}-m_{k\ell}^{\rm
as}\right)\Psi^{\ell}\,,$ (1)

where $\gamma^{\mu}$ are Dirac gamma matrices, $\mu=0,\,1,\,2,\,3$ refers to
the coordinates in Minkowski space, while $I,\,J=1,\,\cdots,\,N_{\rm f}$, and
$k,\,\ell=1,\,\cdots,\,N_{\rm as}$, are flavor indexes—color and spin indexes
are understood. Throughout this paper we assume mass-degenerate fermions, for
which $m_{IJ}^{\rm f}=m^{\rm f}\delta_{IJ}$ and $m_{k\ell}^{\rm as}=m^{\rm
as}\delta_{k\ell}$. The field-strength tensor, $G_{\mu\nu}$, and the covariant
derivative, $D_{\mu}$, acting on the fermions are defined following the
conventions of Ref. Bennett:2019cxd :

$\displaystyle G_{\mu\nu}$
$\displaystyle=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}+ig\left[A_{\mu},A_{\nu}\right]\,,$
(2) $\displaystyle D_{\mu}Q$ $\displaystyle=\partial_{\mu}Q+igA_{\mu}Q\,,$ (3)
$\displaystyle D_{\mu}\Psi$
$\displaystyle=\partial_{\mu}\Psi+igA_{\mu}\Psi+ig\Psi A_{\mu}^{T}\,,$ (4)

where $g$ is the gauge coupling.

An element, $M\in Sp(4)$, of the gauge group acts on the fermion fields with
the gauge transformation $Q\to MQ$ and $\Psi\to M\Psi M^{T}$. Because the
fundamental representation is pseudo-real, while the 2-index antisymmetric one
is real, the global symmetries are enhanced in comparison with a complex
representation (such as that of QCD). The global symmetry of the Lagrangian is
$U(1)_{\rm f}\times U(1)_{\rm as}\times SU(2N_{\rm f})\times SU(2N_{\rm as})$.
One combination of the $U(1)$ factors is broken by the axial anomaly.111Due to
the multi-representation nature of this theory, both $U(1)$ symmetries are
expected to be spontaneously broken which would lead to two additional
(pseudo-) Nambu-Goldstone bosons (PNGBs): one combination of them is broken by
the anomaly, while the other combination is non-anomalous and may have an
implication on the phenomenological studies of composite Higgs models
Belyaev:2015hgo . The bilinear condensate of fundamental fermions breaks
spontaneously the associated $SU(2N_{\rm f})$ symmetry to its $Sp(2N_{\rm f})$
subgroup, while the condensate made of antisymmetric fermions gives rise to
the breaking pattern $SU(2N_{\rm as})\to SO(2N_{\rm as})$ Kosower:1984aw . The
mass terms break explicitly the symmetry along the same pattern, providing
masses for the PNGBs. From hereon, we ignore the $U(1)$ symmetries and specify
the numbers of flavors to $N_{\rm f}=2$ and $N_{\rm as}=3$. We hence have
$5+20$ PNGBs associated with the two non-Abelian cosets.

For the $Sp(2N)$ Euclidean action on the lattice—see Refs. Bennett:2022yfa ;
Bennett:2023wjw ; Bennett:2023gbe  for technical details—we adopt the standard
Wilson plaquette action. We write it in terms of the gauge links,
$U_{\mu}(x)$, as

$\displaystyle S_{g}=\beta\sum_{x}\sum_{\mu<\nu}\left(1-\frac{1}{2N}{\rm
Re}\,{{\rm
Tr}}\,\left(U_{\mu}(x)U_{\nu}(x+\hat{\mu})U_{\mu}^{\dagger}(x+\hat{\nu})U_{\nu}^{\dagger}(x)\right)\right)\,,$
(5)

independently on the representation of the gauge links, and $\mu,\,\nu$ denote
the direction of the link, starting from lattice site $x$, while
$\hat{\mu},\,\hat{\nu}$ are unit displacement on the lattice. The fermions are
described by the Wilson fermion action Wilson:1974sk , for both
representations:

$\displaystyle S_{f}=a^{4}\sum_{J=1}^{N_{\rm
f}}\sum_{x,y}\overline{Q}^{J}(x)D^{\mathrm{(f)}}(x,y)Q^{J}(y)+a^{4}\sum_{k=1}^{N_{\rm
as}}\sum_{x,y}\overline{\Psi}^{k}(x)D^{\mathrm{(as)}}(x,y)\Psi^{k}(y)\,.$ (6)

The lattice spacing is denoted by $a$, while the fundamental and antisymmetric
Dirac operators, $D^{\mathrm{(f)}}$ and $D^{\mathrm{(as)}}$, depending on the
gauge links in the respective representation. We borrow their definitions from
Refs. Bennett:2023gbe :

$\displaystyle
D^{(R)}(x,y)\equiv\left(\frac{4}{a}+m^{R}_{0}\right)\delta(x,y)-\frac{1}{2a}\sum_{\mu=1}^{4}\left\\{(1-\gamma_{\mu})U^{(R)}_{\mu}(x)\delta(x+\hat{\mu},y)+(1+\gamma_{\mu})U^{(R)\dagger}_{\mu}(x)\delta(x-\hat{\mu},y)\right\\}\,,$
(7)

where $R={\rm f},{\rm as}$ denotes the different representations.
Specifically, $m_{0}^{\mathrm{f}}$ and $m_{0}^{\mathrm{as}}$ are the bare
masses of fermions of species ${\rm(f)}$ and ${\rm(as)}$, respectively. For
the link variable, $U^{\mathrm{(f)}}_{\mu}(x)=U_{\mu}(x)$, while a choice of
parameterisation for $U^{\mathrm{(as)}}_{\mu}(x)$ can be found in Ref.
Bennett:2022yfa .

Label | $\beta$ | $am_{0}^{\rm as}$ | $am_{0}^{\rm f}$ | $N_{t}$ | $N_{s}$ | $N_{\rm therm}$ | $n_{\rm skip}$ | $N_{\rm conf}$ | $\langle P\rangle$ | $w_{0}/a$ | $\tau_{\rm int}^{Q}$ | $\bar{Q}$
---|---|---|---|---|---|---|---|---|---|---|---|---
M1 | 6.5 | -1.01 | -0.71 | 48 | 20 | 3006 | 14 | 479 | 0.585172(16) | 2.5200(50) | 6.9(2.4) | 0.38(12)
M2 | 6.5 | -1.01 | -0.71 | 64 | 20 | 1000 | 28 | 698 | 0.585172(12) | 2.5300(40) | 7.1(2.1) | 0.58(14)
M3 | 6.5 | -1.01 | -0.71 | 96 | 20 | 4000 | 26 | 436 | 0.585156(13) | 2.5170(40) | 6.4(3.3) | -0.60(19)
M4 | 6.5 | -1.01 | -0.70 | 64 | 20 | 1000 | 20 | 709 | 0.584228(12) | 2.3557(31) | 10.6(4.8) | -0.31(19)
M5 | 6.5 | -1.01 | -0.72 | 64 | 32 | 3020 | 20 | 295 | 0.5860810(93) | 2.6927(31) | 12.9(8.2) | 0.80(33)
Table 1: Ensembles generated for this study of the $Sp(4)$ theory coupled to
fermions in multiple representations. The inverse coupling is denoted as
$\beta$ and the bare masses of the (Wilson-Dirac) fermions transforming
according to the fundamental and antisymmetric representation of the gauge
group are denoted by $am_{0}^{\mathrm{f}}$ and $am_{0}^{\mathrm{as}}$,
respectively. The lattice volume is $N_{t}N_{s}^{3}a^{4}$. For each ensemble,
we further report the number of gauge configurations, $N_{\rm conf}$, as well
as the average plaquette, $\langle P\rangle$, the gradient flow scale in
lattice units, $w_{0}/a$, the mean of the topological charge, $\bar{Q}$, and
the correlation length of the topological charge, $\tau_{\rm int}^{Q}$. We
also report, for reproducibility purposes, the number of thermalization steps,
$N_{\rm therm}$, discarded from the analysis, as well as the number of
complete sweeps between configurations $n_{\textrm{skip}}$, which have been
discarded to avoid larger autocorrelation. Notice however that $\tau_{\rm
int}^{Q}$ has been computed on the remaining $N_{\textrm{conf}}$
configurations only, and yet $\tau_{\rm int}^{Q}>1$.

Gauge configurations are generated using Grid Boyle:2015tjk ; Boyle:2016lbp ;
Yamaguchi:2022feu , which has the functionality to work with $Sp(2N)$ gauge
groups Bennett:2023gbe . We include dynamical fermions using the Hybrid Monte-
Carlo (HMC) algorithm Duane:1987de  for the two ${\rm(f)}$ fermions and the
rational HMC (RHMC) Clark:2006fx  algorithm for the three ${\rm(as)}$ Dirac
fermions. In principle, the inclusion of an odd number of degenerate fermions
might give rise to a sign problem, but for an odd number of fermions in the
antisymmetric representation the determinant of the Dirac operator remains
positive and real Hands:2000ei ; Bennett:2022yfa . Acceptance rates were tuned
to be around $70\%$ to $80\%$, which corresponds to $27-36$ integrator steps
for a single unit of Monte-Carlo time. Resymplecticization was performed after
every gauge configuration update. We use even-odd preconditioning as in Ref.
Bennett:2022yfa . It was also shown in Ref. Bennett:2023gbe  that the HMC and
RHMC implementations in $Sp(4)$ yield compatible results for even number of
fermions, and that one could equivalently adopt the HMC algorithm for two of
the ${\rm(as)}$ fermions and RHMC for the third, without visible changes to
the results.

Our hypercubic lattices have $N_{s}$ lattice sites in the spatial and
$N_{t}>N_{s}$ lattice sites in the temporal direction, hence the volume is
${V}_{4}=(N_{s}a)^{3}\times(N_{t}a)=L^{3}\times T$. We impose periodic
boundary conditions for the gauge fields. For the fermion fields, we impose
periodic boundary conditions in the spatial dimensions and anti-periodic
boundary conditions along the temporal direction. The lattice action has three
free parameters: the inverse gauge coupling,
$\beta=2N_{c}/g^{2}=4N/g^{2}=8/g^{2}$, and the masses, $m_{0}^{\mathrm{f}}$
and $m_{0}^{\mathrm{as}}$, of the two types of fermions. At strong coupling
(small $\beta$) a transition into an unphysical bulk phase occurs. For this
action it was found that a value of $\beta\gtrsim 6.3$ is sufficient to avoid
this lattice phase Bennett:2022yfa . We choose $\beta=6.5$ for the ensembles
discussed in this paper, and we keep the bare mass of the ${\rm(as)}$ fermions
fixed, to be $am_{0}^{as}=-1.01$. We allow the mass of the ${\rm(f)}$ fermions
to vary over a modest range of values, as listed in Tab. 1. For each ensemble,
we compute the average plaquette, $\langle P\rangle$, defined as

$\displaystyle\langle P\rangle$ $\displaystyle\equiv$
$\displaystyle\frac{1}{6N_{t}N_{s}^{3}}\sum_{x}\sum_{\mu>\nu}{\rm Re}\,{{\rm
Tr}}\,\frac{1}{2N}\left[\frac{}{}U_{\mu}(x)U_{\nu}(x+\hat{\mu})U_{\mu}^{\dagger}(x+\hat{\nu})U_{\nu}^{\dagger}(x)\frac{}{}\right]\,,$
(8)

as this quantity enters the tadpole-improved gauge coupling
$\tilde{g}^{2}\equiv g^{2}/\langle P\rangle$ Martinelli:1982mw ; Lepage:1992xa
.

###  II.2 Gradient flow, topological charge and autocorrelations

In this study, we adopt the gradient flow Luscher:2011bx ; Luscher:2013vga ,
and its lattice implementation, the Wilson flow Luscher:2010iy . This finds
two applications: on one hand, it allows us to set the physical scale for our
ensembles BMW:2012hcm , on the other hand, it will be used as a smoothening
process in the extraction of topological properties. We follow the convention
and processes described in Ref. Bennett:2022ftz . We define a new observable,
$\mathcal{W}(t)$, as a function of a new gradient flow time, $t$, as

$\displaystyle\mathcal{W}(t)$ $\displaystyle\equiv\frac{\rm d}{\rm d\ln
t}\left\\{t^{2}\langle E(t)\rangle\right\\}\,,$ (9)

where in turn222Here we fix a typo in Eq. (23) Ref. Bennett:2022ftz , in which
a minus sign is missing, and which should read as our current Eq. (10), and in
Eq. (25) of the same reference, the right hand side of which should read as
our current Eq. (9).

$\displaystyle E(t,x)\equiv-\frac{1}{2}~{}{\rm
Tr}\,G_{\mu\nu}(t,x)G_{\mu\nu}(t,x).$ (10)

The field strength tensor in Euclidean space-time, $G_{\mu\nu}(t,x)$,
evaluated at non-vanishing flow time $t$, is defined in terms of the five-
dimensional gauge field $A_{\mu}(t,x)$ as

$\displaystyle\frac{{\rm d}A_{\mu}(t,x)}{{\rm
d}t}=D_{\nu}G_{\nu\mu}(t,x),~{}~{}{\rm with}~{}~{}A_{\mu}(t=0,x)=A_{\mu}(x).$
(11)

We define the gradient flow scale, $w_{0}$, as the square root of the flow
time, $t$, for which $\mathcal{W}(t)|_{t=w_{0}^{2}}=\mathcal{W}_{0}=0.35$ and
report all dimensionful quantities in units of $w_{0}$ Fodor:2012td .
Dimensionful quantities expressed in units of $w_{0}$ are denoted by a hat,
for example masses are denoted as $\hat{m}\equiv(ma)(w_{0}/a)$. The measured
values of $w_{0}/a$, obtained on the lattice by discretizing Eq. (11), using
the clover discretization of the field-strength tensor, are tabulated in Tab.
1.

We monitor the topological charge, $Q$, of each configuration, which in the
continuum is defined as

$\displaystyle Q=\frac{1}{32\pi^{2}}\int{\rm
d}^{4}x\,\epsilon^{\mu\nu\rho\sigma}\,{\rm Tr}\,G_{\mu\nu}G_{\rho\sigma}.$
(12)

On the lattice, we measure $Q$ by following the same process as in Ref.
Bennett:2022ftz , to which we refer the reader for details. By applying the
gradient flow to smoothen the gauge fields, UV fluctuations are removed in the
measurement of $Q$. None of the ensembles used in this study show hard
evidence of topological freezing, with the Monte Carlo algorithm sampling
configurations with multiple values of $Q$. In order to quantify this
statement we study the autocorrelations in Monte-Carlo time of both the
topological charge, $Q$, as well as the average plaquette, $\langle P\rangle$.
We perform the measurements of the gradient flow scale, the topological charge
and the hadron correlators using the HiRep code DelDebbio:2008zf ; HiRepSUN ,
which has been extended to symplectic gauge groups HiRepSpN . To this purpose,
we convert the configurations created with the Grid code using the Gauge Link
Utility (GLU) library GLU .

Figure 1: Thermalised trajectory of the average plaquette, $\langle P\rangle$,
in Monte-Carlo time (left) and histogram of the same data (right), for the
lattice ensemble denoted as M$1$ in Table 1. All other ensembles show
qualitatively similar behaviors. We fitted the histogram to a Gaussian
distribution (red line on the right). Configuration sampling and binning of
the raw data are discussed in the main text.

For a generic observable of interest, $X$, we denote as $\tau=1,\,\dots,\,N$,
the Monte-Carlo time, as $X_{i}$ the individual measurement of the observable,
and as $\bar{X}$ the arithmetic mean of $X$. The Madras-Sokal integrated
autocorrelation time, $\tau_{\rm int}$, is defined as follows Madras:1988ei ;
Wolff:2003sm ; Luscher:2004pav :

$\displaystyle\tau_{\rm int}$ $\displaystyle=$ $\displaystyle\max_{\tau_{\rm
max}}\tau_{\rm int}(\tau_{\rm max})\,,$ (13)

where

$\displaystyle\tau_{\rm int}(\tau_{\rm max})$ $\displaystyle=$
$\displaystyle\frac{1}{2}+\sum_{t=1}^{\tau_{\rm max}}\Gamma(\tau)\,,$ (14)

and

$\displaystyle\Gamma(\tau)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N-\tau}\frac{\left(X_{i}-\bar{X}\right)\left(X_{i+\tau}-\bar{X}\right)}{N-\tau}\,.$
(15)

In applying these definitions, we assume that the Monte-Carlo time series,
$X_{i}$, is fully thermalized. In practice, we vary also the thermalization
cut-off and choose the thermalization time so that histograms of the plaquette
and topological charge show a Gaussian behavior. We use this definition by
treating two values of $X_{i}$ as uncorrelated if they are separated by at
least $2\tau_{\rm int}$ in Monte-Carlo time $\tau_{i}$. For comparison, we
additionally compute the exponential autocorrelation time, $\tau_{\rm exp}$,
by fitting the autocorrelation function, $C_{X}(\tau)$, defined as

$\displaystyle C_{X}(\tau)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N-\tau}\left(X_{i}-\bar{X}\right)\left(X_{i+\tau}-\bar{X}\right)\,,$
(16)

to an exponential decay:

$\displaystyle\frac{C_{X}(\tau)}{C_{X}(0)}$ $\displaystyle\sim$
$\displaystyle\exp\left(-\frac{\tau}{\tau_{\rm exp}}\right)\,.$ (17)

Figure 2: Thermalised trajectory of the topological charge, $Q$, in ensemble
M2, restricted to configurations used in the remainder of this paper, after
removing the effect of autocorrelation in the plaquette. We fitted the
histogram to a Gaussian distribution (red line on the right). Gauge
configurations show a modest level of residual autocorrelation, with respect
to the topological charge. The other ensembles exhibit qualitatively similar
behaviors.

In Fig. 1 we display one example of average plaquette trajectory in Monte-
Carlo time. We find the exponential autocorrelation to be always smaller than
the integrated Madras-Sokal autocorrelation time. We keep one gauge
configuration every $n_{\rm skip}$ Monte-Carlo time units, such that $n_{\rm
skip}\gtrsim\tau_{\rm int}^{\langle P\rangle}$ for the autocorrelation time of
the average plaquette. We measure our observables on these gauge
configurations and bin the resulting dataset with bin size of $2$. We retain
$N_{\rm conf}$ measurements, or $N_{\rm conf}/2$ independent samples. We
tabulate $N_{\rm conf}$, in Tab. 1.

Label | Interpolating operator $\mathcal{O}^{A}$ | $J^{P}$ | $Sp(4)$ | $SO(6)$
---|---|---|---|---
PS | $\bar{Q}^{I}\gamma_{5}Q^{J}$ | $0^{-}$ | $5$ | $1$
V | $\bar{Q}^{I}\gamma_{i}Q^{J}$ | $1^{-}$ | $10$ | $1$
T | $\bar{Q}^{I}\gamma_{0}\gamma_{i}Q^{J}$ | $1^{-}$ | $10$ | $1$
AV | $\bar{Q}^{I}\gamma_{5}\gamma_{i}Q^{J}$ | $1^{+}$ | $5$ | $1$
AT | $\bar{Q}^{I}\gamma_{0}\gamma_{5}\gamma_{i}Q^{J}$ | $1^{+}$ | $10$ | $1$
S | $\bar{Q}^{I}Q^{J}$ | $0^{+}$ | $5$ | $1$
ps | $\bar{\Psi}^{k}\gamma_{5}\Psi^{\ell}$ | $0^{-}$ | $1$ | $20$
v | $\bar{\Psi}^{k}\gamma_{i}\Psi^{\ell}$ | $1^{-}$ | $1$ | $15$
t | $\bar{\Psi}^{k}\gamma_{0}\gamma_{i}\Psi^{\ell}$ | $1^{-}$ | $1$ | $15$
av | $\bar{\Psi}^{k}\gamma_{5}\gamma_{i}\Psi^{\ell}$ | $1^{+}$ | $1$ | $20$
at | $\bar{\Psi}^{k}\gamma_{0}\gamma_{5}\gamma_{i}\Psi^{\ell}$ | $1^{+}$ | $1$ | $15$
s | $\bar{\Psi}^{k}\Psi^{\ell}$ | $0^{+}$ | $1$ | $20$
Table 2:  Meson interpolating operators, $\mathcal{O}^{A}(x)$, built using
Dirac fermions transforming in the fundamental, $Q^{I}$, and antisymmetric,
$\Psi^{k}$, representations. Here, $I,\,J=1,\,2$ and $k,\,\ell=1,\,2,\,3$ are
flavor indices, and we restrict our attention to the case $I\neq J$ and
$k\neq\ell$. The index $i=1,\,2,\,3$ refers to spatial directions. Color and
spinor indices are implicit and summed over. $J^{P}$ are space-time quantum
numbers. We report also the representation under the unbroken $Sp(4)$ and
$SO(6)$ global symmetries.

We then compute the topological charge, $Q$, for the $N_{\rm conf}$ gauge
configurations used in the remainder of this paper. As for the plaquette, we
plot the trajectory, compute the autocorrelation time(s), and display the
measurements in a histogram. One example of the results is shown in Fig. 2. We
report the autocorrelation of the topology, $\tau_{\rm int}^{Q}$, and the
average topological charge, $\bar{Q}$, in Tab. 1. Qualitatively, all resulting
distributions are Gaussian. Yet, $\tau_{\rm int}^{Q}>1$ and a non-vanishing
average value of $Q$ is measured for all ensembles. We conclude that the
configurations used in this paper are affected by a moderate amount of
residual autocorrelation in the topological charge.333In pure gauge theory, no
sizeable effect on the glueball spectrum was detected even in ensembles with
complete topological freezing Bonanno:2022yjr .

##  III Correlation functions

Table 2 summarises the properties of the meson operators, $\mathcal{O}^{A}$,
of interest in this paper. They are constituted by two fermions in the
fundamental, ${\rm(f)}$, representation $Q^{I}$, where $I=1,\,2$, or two in
the antisymmetric, ${\rm(as)}$, representation, denoted $\Psi^{k}$, where
$k=1,\,2,\,3$. We label them as pseudoscalar, vector, tensor, axial-vector,
axial-tensor and scalar, both for mesons made of ${\rm(f)}$ fermions and
${\rm(as)}$ ones—in the latter case, we conventionally label them with lower
case acronyms. We also display the quantum numbers, $J^{P}$, and the
irreducible representation in the unbroken $Sp(4)\times SO(6)$ global
symmetry.

Zero-momentum, two-point correlation functions are defined by averaging over
lattice sites as:

$C^{AB}(t)=\sum_{\vec{x}}\,\langle\mathcal{O}^{A}(t,\vec{x})\,\bar{\mathcal{O}}^{B}(0)\rangle\equiv\langle\mathcal{O}^{A}(t)\bar{\mathcal{O}}^{B}(0)\rangle\,.$
(18)

For $A\neq B$, Equation (18) can be recast for a generic $C(t)$ as a summation
over a complete basis, $|n\rangle\langle n|$, with associated energies
$E_{n}$, as follows:

$C^{AB}(t)=\sum_{n}\dfrac{\langle 0|\mathcal{O}^{A}(0)|n\rangle\langle
n|\bar{\mathcal{O}}^{B}(0)|0\rangle}{2E_{n}}\,b(t,E_{n})\,,$ (19)

where $b(t,E_{n})$ becomes an exponential function in the limit of infinite
temporal lattice extent,
$b(t,E_{n})\xrightarrow[N_{t}\to\infty]{}e^{-t\,E_{n}}$. For a finite,
periodic lattice in the time extent, one expects the following functional form

$b(t,E_{n})=b_{+}(E_{n})e^{-t\,E_{n}}+b_{-}(E_{n})e^{-(N_{t}-t)E_{n}}\,,$ (20)

which suggests identifying the interval for which the ground states dominate
by examining the effective mass, defined as

$a\,m_{\rm
eff}(t)\equiv\cosh^{-1}\left[\dfrac{C^{AB}(t+a)+C^{AB}(t-a)}{2C^{AB}(t)}\right]\,.$
(21)

At large Euclidean times, one expects $a\,m_{\rm
eff}(t)\xrightarrow[t\to\infty]{}a\,E_{0}$. To extract the ground state energy
$E_{0}$, therefore, one has the freedom of fitting the plateau according to
Eq. (21) or the correlators as in Eq. (19). In the present paper, we fitted
correlators setting coefficients $b_{+}(E_{n})=b_{-}(E_{n})=1$ when $A=B$, and
$b_{+}(E_{n})=1,\,b_{-}(E_{n})=0$ when $A\neq B$ over the interval that shows
a plateau in the effective mass.

In the remainder of the section, we describe two techniques, APE and Wuppertal
smearing, that we adopt in order to optimize the extraction of physical
information in the fitting procedure. Moreover, we introduce an additional,
well established methodology in lattice field theory, the GEVP, that reduces
the contributions from the excited states to ground states at large $t$.

###  III.1  APE and Wuppertal smearing algorithms

From Eq. (19), one expects contamination from excited states to affect meson
two-point functions at moderate-to-small time separation. Conversely, for any
finite lattice, there is an intrinsic limitation on the maximum length of the
time separation in the two-point functions. The combination of these two
factors results in a limitation on the length of the plateau displaying
(approximately) constant effective mass—see Eq. (21)—that can be used for
spectroscopy. As visible in Eq. (19), the contribution of the $n^{th}$ excited
state, relative to the ground state, is determined by two factors: the
exponential decays proportional to $\exp\left[-(E_{n}-E_{0})t\right]$, and the
overlap functions encoded in the matrix elements $\langle
0|\mathcal{O}(0)|n\rangle$. Increasing the overlap of the interpolating
operators, $\mathcal{O}$, with the ground state, $|0\rangle$, relative to the
excited states, $|n\rangle$, results in suppressed contamination from
excitations, longer plateaux in effective mass plots, and more precise
spectroscopy.

Our implementation of these ideas combines APE smearing APE:1987ehd ;
Falcioni:1984ei  of the gauge configurations with Wuppertal smearing
Gusken:1989qx ; Roberts:2012tp ; Alexandrou:1990dq  of the sink/source
operators. APE smearing acts on gauge links, smoothening their ultraviolet
fluctuations, and hence improving the statistical control over the plateaux
being analyzed. Wuppertal smearing modifies the fermion fields used to source
the two-point functions, and increases the ground state overlap by using
extended interpolators, instead of point-like ones.

APE smearing consists of applying an iterative procedure involving the staple
operator around each gauge link,
$S_{\mu}(x)\equiv\sum_{\pm\nu\neq\mu}U_{\nu}(x)U_{\mu}(x+\hat{\nu})U^{\dagger}_{\nu}(x+\hat{\mu})$,
as follows

$U^{(m)}_{\mu}(x)=\mathcal{P}\left\\{(1-\alpha_{\textrm{APE}})U^{(m-1)}_{\mu}(x)+\frac{\alpha_{\textrm{APE}}}{6}S^{(m-1)}_{\mu}(x)\right\\}\,,$
(22)

with the initial conditions $U_{\mu}^{0}=U_{\mu}$ and $S^{0}_{\mu}=S_{\mu}$.
The iteration number is $m=1,\,\dots,\,N_{\rm APE}$ and
$\alpha_{\textrm{APE}}$ is called APE-smearing step size. In the case of
interest in this study, we use the same process for gauge links transforming
in the fundamental and antisymmetric representations. As the gauge links at
each iteration are summed over their neighboring staples according to Eq.
(22), they do not necessarily lie within the group manifold. Therefore, we use
a projection operator, $\mathcal{P}$, to project the smeared link variable
into the group manifold. The explicit form of the projector depends on the
gauge group and representation considered.

In order to illustrate the Wuppertal smearing of source and sink operators,
one starts at first with the Dirac equation for point-like source and sink:

$\sum_{y,\,\beta,\,b}D^{R}_{a\alpha\,b\beta}(x,y)\,S^{b\beta}_{R\,c\gamma}(y,0)=\delta_{x,0}\,\delta_{\alpha\gamma}\,\delta_{ac}\,,$
(23)

where $D^{R}_{a\alpha\,b\beta}$ is the Wilson-Dirac operator in the
representation, $R$, spinor indices are denoted as $\alpha,\,\beta,\,\gamma$,
and (generalised) color indexes as $a,\,b,\,c$, respectively. The solution,
$S^{b\beta}_{R\,c\gamma}$, is the hyperquark propagator in such
representation, $R$. The two-point correlation function is then written as

$C(t)=\langle\mathcal{O}_{R}(t)\,\bar{\mathcal{O}}_{R}(0)\rangle=\Big{\langle}\sum_{\vec{x}}{\rm
Tr}\,\left[\left(\Gamma\,S_{R}(x,0)\,\bar{\Gamma}\,S_{R}(x,0)\right)\right]\Big{\rangle}\,,$
(24)

where $\Gamma$ and $\bar{\Gamma}$ depend on the spin structure of the
interpolating operator, $\mathcal{O}_{R}$, listed in Tab. 2.

Wuppertal smearing consists of replacing $\delta_{x,0}$, in the right-hand
side of Eq. (23), with a function $q_{R}^{(n+1)}(x)$, defined through an
iterative procedure based on a diffusion process:

$\displaystyle q_{\rm R}^{(n)}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{1+6\varepsilon_{\rm R}}\left[q_{\rm
R}^{(n-1)}(x)+\varepsilon_{\rm R}\sum_{\mu=\pm 1}^{\pm 3}U^{\rm
R}_{\mu}(x)q_{\rm R}^{(n-1)}(x+\hat{\mu})\right]\,,\qquad(n>0)\,,$ (25)
$\displaystyle q_{\rm R}^{(0)}(x)$ $\displaystyle=$
$\displaystyle\delta_{x,0}\,.$ (26)

Here, $\varepsilon_{\rm R}$ is the representation-dependent Wuppertal-smearing
step size. Each source smearing requires an inversion of the Dirac operator,
and the result is a new, source-smeared, propagator denoted as
$S^{(n)}_{R}(y,\,0)$. Sink smearing is obtained by applying the smearing
iteration of Eq. (25) to the source-smeared propagator. Note, that this does
not require further inversions. We denote the propagator with $N_{\rm source}$
iterations of source smearing and $N_{\rm sink}$ iterations of sink-smearing
as $S^{R}_{(N_{\rm source},N_{\rm sink})}(x,\,0)$. The smeared two-point
functions read

$C_{N_{\rm source},\,N_{\rm sink}}(t)=\Big{\langle}\sum_{\vec{x}}{\rm
Tr}\left[\Gamma\,\left(S^{R}_{(N_{\rm source},N_{\rm
sink})}(x,\,0)\right)\,\bar{\Gamma}\,\left(S^{R}_{(N_{\rm source},N_{\rm
sink})}(x,\,0)\right)\right]\Big{\rangle}\,.$ (27)

The measurement of APE and Wuppertal smeared two-point mesonic correlation
functions are performed using the HiRep code DelDebbio:2008zf ; HiRepSUN ;
HiRepSpN . The APE and Wuppertal smearing step-sizes,
$\alpha_{\textrm{APE}},\,\varepsilon_{\rm f},\,\varepsilon_{\rm as}$, and the
number of iterations, $N_{\rm APE},\,N_{\rm source},\,N_{\rm sink}$, are
parameters that are tuned to optimize the signal, both by improving the
effective mass plateaux, as well as the resolution of peaks in the spectral
density reconstruction—as we will discuss later. For our purposes, for
ensembles M1 to M4, we apply Wuppertal smearing step-sizes $\varepsilon_{\rm
f}=0.20$ for the fundamental sector, and $\varepsilon_{\rm as}=0.12$ for the
antisymmetric one. For the ensemble M5, $\varepsilon_{\rm f}=0.24$ for the
fundamental sector and $\varepsilon_{\rm as}=0.12$ for the antisymmetric one.
Typical iteration numbers for sink and sources are $N_{\rm sink},\,N_{\rm
source}\in[0,80]$. For APE smearing, we apply $\alpha_{\textrm{APE}}=0.4$ and
$N_{\rm APE}=50$. We verified explicitly that the process did not lead to
large changes in the mass of the pseudoscalar ground state, and that the
convexity of the effective mass plots is not altered, hence excluding the
possibility of oversmearing.

###  III.2 Generalized Eigenvalue Problem

The ground state in a given channel, $E_{0}$, can be identified by the plateau
in effective mass, at least as long as the ground state is clearly separated
in mass from its excited states. A pragmatic way to isolate excitations
involves performing multi-functional fits by minimizing a correlated chi-
square functional

$\chi^{2}=\sum_{t,\,t^{{}^{\prime}}}\left(h^{(k)}(t)-C(t)\right)\text{Cov}^{-1}_{t\,t^{{}^{\prime}}}[C]\left(h^{(k)}(t^{{}^{\prime}})-C(t^{{}^{\prime}})\right)\,,$
(28)

where the fitting function is the periodic version of a multi-exponential fit,

$h^{(k)}(t)=\sum_{n=1}^{k}\mathcal{B}_{n}\,b(t,E_{n})\,,$ (29)

while $b(t,E)$ is defined in Eq. (20). Yet, the extraction of such
excitations, $E_{n}$, is hindered by the increasing number of degrees of
freedom it requires, and the signal gets exponentially suppressed with growing
$n$.

A variational approach based on solving a _generalized eigenvalue problem_
(GEVP) can overcome this problem Blossier:2009kd . One starts by defining a
matrix-valued correlation function, $\mathcal{C}_{ij}(t)$, encompassing the
Euclidean-space two-point correlation functions of a set of operators, ${\cal
O}_{i}$:

$\mathcal{C}_{ij}(t)=\langle\mathcal{O}_{i}(t)\bar{\mathcal{O}}_{j}(0)\rangle,\quad
i,j=0,\dots,N-1\,.$ (30)

Under the assumption that non-degenerate energy levels exist, they can be
ordered with $E_{n}<E_{n+1}$, where the eigenvalues are labelled as
$n=0,\,\dots,\,N-1$ In this notation, $E_{0}$ denotes the ground state energy.
The GEVP is defined as

$\mathcal{C}(t)v_{n}(t,t_{1})=\lambda_{n}(t,t_{1})\mathcal{C}(t_{1})v_{n}(t,t_{1})\,,$
(31)

in which the number of new functions, $v_{n}$ and $\lambda_{n}$, matches the
dimension of the variational basis used to build the matrix $\mathcal{C}$. By
fixing a reference value, $t_{1}$, Eq. (31) can be solved as an eigenvalue
equation for each lattice-time slice $t>t_{1}$. This procedure results in
determining the eigenvalues, $\lambda_{n}(t,t_{1})$.

In the next step, the eigenvalue functions, $\lambda_{n}(t,t_{1})$, are used
to define $N$ effective-mass plateaux, $am_{n}^{\rm eff}(t,t_{1})$:

$m^{\rm eff}_{n}(t,t_{1})=-\dfrac{\partial\log\lambda_{n}(t,t_{1})}{\partial
t}=-\dfrac{1}{a}\left\\{\log\left[\lambda_{n}(t+a,t_{1})\right]-\log\left[\lambda_{n}(t,t_{1})\right]\right\\}\,,$
(32)

or, in the more realistic case of a periodic lattice:

$m^{\rm
eff}_{n}(t,t_{1})=\dfrac{1}{a}\,\cosh^{-1}\left[\dfrac{{\lambda_{n}(t+a,t_{1})+\lambda_{n}(t-a,t_{1})}}{2\lambda_{n}(t,t_{1})}\right]\,.$
(33)

These quantities are expected to converge to the energy levels,

$aE_{n}=\lim_{t\to\infty}am_{n}^{\rm eff}(t,t_{1})\quad n=0,\,\dots,\,N-1\,.$
(34)

In the case of lattice correlation functions, to solve the GEVP in Eq. (31)
one focuses on late-time slices, so that the corrections due to higher-energy
excitations, with $E>E_{N-1}$, are suppressed. As shown in Ref. LUSCHER1990222
, for fixed $t_{1}$ such contributions are expected to take the form

$am^{\rm
eff}_{n}(t,t_{1})=aE_{n}+\mathcal{O}(e^{-\Delta(aE_{n})t}),\quad\Delta(aE_{n})=\min_{m\neq
n}|aE_{m}-aE_{n}|\,,$ (35)

so that the size of the contaminations depends on the gap between the energy
level of the spectrum, and the size of $t$. As discussed in Blossier:2009kd
there are systematic effects proportional to $\mathcal{A}\cdot\exp(E_{N}\cdot
t_{1})$. On general grounds, then, one expects the quality of the results and
signal in each plateau to improve when considering lattices with larger time
extents, $N_{t}$. In our study, we build a variational basis by considering
different levels of Wuppertal smearing for the interpolating source and sink
operators in each mesonic channel. For the pseudoscalar, axial-vector, axial-
tensor, and scalar channels (and similar for mesons made of ${\rm(as)}$
fermions) we use three different levels of smearing: $N_{\rm source},\,N_{\rm
sink}=0,\,40,\,80$. Therefore, in this case, the correlation matrix
$\mathcal{C}$ is

$\mathcal{C}^{i}(t)=\begin{pmatrix}C^{i}_{0,\,0}(t)&C^{i}_{0,\,40}(t)&C^{i}_{0,\,80}(t)\\\
C^{i}_{40,\,0}(t)&C^{i}_{40,\,40}(t)&C^{i}_{40,\,80}(t)\\\
C^{i}_{80,\,0}(t)&C^{i}_{80,\,40}(t)&C^{i}_{80,\,80}(t)\end{pmatrix}$ (36)

with
$C^{i}_{nm}(t)=\langle\mathcal{O}^{i}(t)\bar{\mathcal{O}}^{i}(0)\rangle_{nm}$,
where
$i=\text{PS},\,\text{AV},\,\text{AT},\,\text{S},\,\text{ps},\,\text{av},\,\text{at},\,\text{s}$
and $n,\,m$ indicate the level of Wuppertal smearing applied to source and
sink operators, $n=N_{\rm source},\,m=N_{\rm sink}$.

In the case of V and T channels, the operators transform in the same way under
the unbroken symmetry groups, and are expected to mix, and source the same
spectrum. Of course, this argument is no longer true for our measurements on a
discretized lattice as the rotational symmetry is broken. Yet, we found no
discernible difference between these two channels for given statistical
errors. Therefore, we extend the correlation matrix to include the cross-
channels $V/T$ and $T/V$, resulting in the following

$\mathcal{C}^{V,\,T}(t)=\begin{pmatrix}C^{V}_{0,\,0}(t)&C^{V}_{0,\,40}(t)&C^{V}_{0,\,80}(t)&C^{V/T}_{0,\,0}(t)&C^{V/T}_{0,\,40}(t)&C^{V/T}_{0,\,80}(t)\\\
C^{V}_{40,\,0}(t)&C^{V}_{40,\,40}(t)&C^{V}_{40,\,80}(t)&C^{V/T}_{40,\,0}(t)&C^{V/T}_{40,\,40}(t)&C^{V/T}_{40,\,80}(t)\\\
C^{V}_{80,\,0}(t)&C^{V}_{80,\,40}(t)&C^{V}_{80,\,80}(t)&C^{V/T}_{80,\,0}(t)&C^{V/T}_{80,\,40}(t)&C^{V/T}_{80,\,80}(t)\\\
C^{T/V}_{0,\,0}(t)&C^{T/V}_{0,\,40}(t)&C^{T/V}_{0,\,80}(t)&C^{T}_{0,\,0}(t)&C^{T}_{0,\,40}(t)&C^{T}_{0,\,80}(t)\\\
C^{T/V}_{40,\,0}(t)&C^{T/V}_{40,\,40}(t)&C^{T/V}_{40,\,80}(t)&C^{T}_{40,\,0}(t)&C^{T}_{40,\,40}(t)&C^{T}_{40,\,80}(t)\\\
C^{T/V}_{80,\,0}(t)&C^{T/V}_{80,\,40}(t)&C^{T/V}_{80,\,80}(t)&C^{T}_{80,\,0}(t)&C^{T}_{80,\,40}(t)&C^{T}_{80,\,80}(t)\\\
\end{pmatrix}$ (37)

where the cross-channels correlators are defined as
$C^{V/T}_{nm}(t)=\langle\mathcal{O}^{V}(t)\bar{\mathcal{O}}^{T}(0)\rangle_{nm}$,
and
$C^{T/V}_{nm}(t)=\langle\mathcal{O}^{T}(t)\bar{\mathcal{O}}^{V}(0)\rangle_{nm}$
The enlarged variational basis can allow for resolving higher excitations.

##  IV Spectral densities

In this section, we introduce and critically appraise the spectral density
reconstruction algorithm and its application to meson spectroscopy. We assess
its dependence on the energy smearing kernel and on the APE and Wuppertal
smearing present in the input data—the two-point correlation function. To this
extent, we implemented this technology in the Python software package
LSDensities, and made it publicly available in Ref. Forzano:2024 . Specific
details of the latter are described in the Appendix. We pay particular
attention to estimating systematic effects in the reconstruction procedure,
and to optimising parameter choices.

###  IV.1 The Hansen-Lupo-Tantalo (HLT) method

Given a generic two-point correlation function, $C(t)$, the spectral density,
$\rho(E)$, is defined as follows:

$C(t)=\int_{E_{\rm min}}^{\infty}dE\,\rho(E)\,b(t,E)\,,$ (38)

where $b(t,E)$ has been introduced in Eq. (20), for periodic time extent. This
definition reduces to an inverse Laplace transform for infinite time extent,
when one chooses vanishing $E_{\rm min}$. This lower bound in the integration
can be chosen between zero and the energy corresponding to the ground state of
the theory in the channel defined by the interpolating operators yielding
$C(t)$.

The input data takes the form of a set of measurements, labeled by
$k=1,\,\dots,\,N_{m}$ (where $N_{m}\leq N_{\textrm{conf}}$ of Tab. 1), of the
given correlation function, $C_{k}(t)$. On the lattice, time is discretized,
so that $t=1,\,\cdots,\,t_{\rm max}\leq N_{t}/2$. The $t_{\rm max}\times
t_{\rm max}$ covariance matrix is

$\hbox{Cov}_{tt{\prime}}[C]\equiv\dfrac{1}{N_{m}}\sum_{k=1}^{N_{m}}\,\Big{(}C_{k}(t)-\langle
C(t)\rangle\Big{)}\Big{(}C_{k}(t^{\prime})-\langle
C(t^{\prime})\rangle\Big{)}\,,$ (39)

where $\langle C(t)\rangle$ is the arithmetic average over the available
measurements at given $t$.

In the literature, several approaches have been applied to reconstruct
spectral densities starting from a finite set of measurements, $C(t)$, which
are affected by noise, circumstances under which the inversion of Eq. (38) is
an ill-posed problem. Among them, the Backus-Gilbert algorithm was originally
devised in Ref. Backus:1968svk , and then modified in Ref. Hansen:2019idp  to
better suit the context of lattice simulations. We refer to this improved
version of the Backus-Gilbert algorithm as Hansen-Lupo-Tantalo (HLT) method.

The starting point is the introduction of smeared spectral
densities444Smearing of spectral densities should not be confused with (APE
and Wuppertal) smearing of two-point correlation functions., via the rewriting

$\rho_{\sigma}(\omega)\equiv\int_{E_{\rm
min}}^{\infty}\,dE\,\Delta_{\sigma}(E-\omega)\,\rho(E)\,,$ (40)

which consists of a convolution of the original spectral density, $\rho(E)$,
with a smearing kernel, $\Delta_{\sigma}(E-\omega)$. The parameter $\sigma$
characterizes the smearing radius around the point $\omega$. At non-zero
smearing radius, the convolution defined in Eq. (40) is such that the smeared
spectral density is always a smooth function.

In a finite volume, the spectral density is a sum of Dirac $\delta$ functions,
corresponding to the discrete eigenvalues of the Hamiltonian:

$\rho(E)=\sum_{n}k_{n}\,\delta(E-\omega_{n})\,,$ (41)

where the sum runs over the eigenvalues $\omega_{n}$, and the coefficients,
$k_{n}$, depends on the lattice spatial volume $N_{s}a$. The choice of the
smearing kernel is guided by the requirement that it is a smooth function,
approaching a Dirac $\delta$ function as $\sigma\to 0$:

$\Delta_{\sigma}(E-\omega)\xrightarrow[\sigma\to 0]{}\delta(E-\omega)\,,$ (42)

and hence in this limit, the regulator disappears:

$\rho_{\sigma}(\omega)\xrightarrow[\sigma\to 0]{}\rho(E=\omega)\,.$ (43)

A key idea of the HLT method is to choose and fix a given smearing kernel at
the beginning of the procedure. If the chosen smearing kernel can be
represented as an infinite sum over the space generated by the basis
functions:

$\Delta_{\sigma}(E-\omega)=\sum_{t=1}^{\infty}g_{t}(\omega,\sigma)\,b(t,E),$
(44)

in the HLT method we can look for the coefficients $g_{t}$ that provide the
best approximation at a finite $t_{\rm max}\leq N_{t}/2$:

$\bar{\Delta}_{\sigma}(E-\omega)=\sum_{t=1}^{t_{\rm
max}}g_{t}(\omega,\sigma)\,b(t,E)\,,$ (45)

In the HLT procedure, the coefficients $g_{t}(\omega,\sigma)$ are defined by
the minimum value of the following functional:

$A[\vec{g}]\equiv\int_{E_{\rm min}}^{\infty}dE\,e^{\alpha
E}\,|\bar{\Delta}_{\sigma}(E-\omega)-\Delta_{\sigma}(E-\omega)|^{2},$ (46)

which measures the difference between the reconstructed smearing kernel,
$\bar{\Delta}_{\sigma}(E-\omega)$, and the target kernel,
${\Delta}_{\sigma}(E-\omega)$. The unphysical $\alpha$ parameterises different
choices of norm DelDebbio:2022qgu , and we shall discuss it later. If the data
are known with infinite precision, the inversion of Eq. (46) is enough to
provide the best values for the coefficients $g_{t}$ defining the smearing
kernel. In any realistic case, with data affected by uncertainties, the
minimization of Eq. (46) amounts to the inversion of a highly ill-conditioned
matrix. As suggested in Ref. Backus:1968svk , the problem can be regularised
by adding a second functional:

$B[\vec{g}]\equiv\sum_{t,\,t^{\prime}=1}^{t_{\rm
max}}\,g_{t}\,\hbox{Cov}_{tt^{\prime}}[C]\,g_{t^{\prime}}\,,$ (47)

where the $t_{\rm max}\times t_{\rm max}$ covariance matrix,
$\hbox{Cov}_{tt^{\prime}}[C]$, is taken from the input data.

With all of the above in place, the third step is to define the functional
$W[\vec{g}]$ Hansen:2019idp :

$W[\vec{g}]\equiv\dfrac{A[\vec{g}]}{A[0]}+\lambda\dfrac{B[\vec{g}]}{B_{\rm
norm}(\omega)}\quad,\quad\lambda\in(0,\infty)\,,$ (48)

where $B_{\rm norm}(\omega)=C^{2}(1)/\omega^{2}$ is written in terms of the
normalization of the correlator at the initial time slice $C(t=1)$, and
$\lambda$ we refer to as the trade-off parameter. The second part of
$W[\vec{g}]$ is called the _statistical error functional_ , and it is
introduced to regularize the problem.

By minimizing $W[\vec{g}]$ for any given value of $\omega$, one can determine
a set of coefficients
$\vec{g}(\omega)=(g_{1}(\omega,\sigma),\,\cdots,\,g_{t_{\rm
max}}(\omega,\sigma))$ which corresponds to the following estimator for the
smeared spectral density:

$\hat{\rho}_{\sigma}(\omega)=\sum_{t=1}^{t_{\rm
max}}g_{t}(\omega,\sigma)\,C(t)\,.$ (49)

We devote the rest of this subsection to a critical discussion of the various
systematics associated with the minimization of $W$. The ratio
$\sqrt{A[\vec{g}]/A_{0}}$ provides an indication of the size of the systematic
error due to the HLT method for extracting $g_{t}(\omega,\sigma)$ with finite
$t_{\rm max}$ and $N_{m}$, as this quantity describes the relative deviation
between the targeted smearing kernel, $\Delta_{\sigma}(E-\omega)$, and the
reconstructed one, $\bar{\Delta}_{\sigma}(E-\omega)$. Conversely, the
$B[\vec{g}]/B_{\rm norm}$ part of the functional $W[\vec{g}]$ provides an
estimate of the statistical uncertainty for the reconstructed spectral density
$\hat{\rho}_{\sigma}(\omega)$. To illustrate this point, in Fig. 3, we display
the ratio $A[\vec{g}]/A_{0}$, evaluated at the minimum of $W$ for a given
choice of $\lambda$ (and $\alpha$). Systematic effects due to the
reconstruction are unsuppressed when $A[\vec{g}]/A_{0}$ is larger, which
corresponds to choices of the parameters $\lambda$ and $\alpha$ that make the
processing of information ineffective. Indeed, in this regime changes of
$\lambda$ result in sizeable changes in $\hat{\rho}_{\sigma}(\omega)$. In this
case, different numerical choices of norm parameter, $\alpha$, yield results
that are not compatible with one another within statistical uncertainties.

Figure 3: Examples of the plateaux in spectral density reconstruction for a
fixed value of energy. In the top panel, the variation of the trade-off
parameter $\lambda$ is shown on the x-axis, whereas $A_{0}=A[\vec{g}=0]$ is
varied in the bottom panel. The correspondent values of reconstructed spectral
densities $\hat{\rho}_{\sigma}(\omega)$ (shown on y-axes of both panels) are
within statistical uncertainties of the unphysical parameters $\lambda$ and
$\alpha$. Some representative values of $\alpha$ and $\lambda$ are chosen in
the figure. The data correspond to the fundamental vector meson channel of
ensemble M1, using $t_{\rm max}=N_{t}/2$, $\alpha=0$,
$\sigma=0.30m_{\mathrm{V}}$.

Conversely, we can reduce the trade-off parameter, $\lambda$, reducing the
size of $A[\vec{g}]/A_{0}$, and hence yielding smaller systematic effects.
This comes at the price of larger statistical uncertainties, as can be seen in
Fig. 3. This is expected on the basis of the very definition of the functional
$W[\vec{g}]$ in Eq. (48), as having a small trade-off parameter, $\lambda$,
corresponds to minimising predominantly the $A[\vec{g}]/A_{0}$ functional.
Therefore, considering smaller $\lambda$ values corresponds to forcing the
systematic error to be smaller, and in principle the reconstruction more
accurate, but at the same time the functional $B[\vec{g}]/B_{\rm norm}$ is
poorly constraining the system. In order to find a value of lambda that sits
on middle grounds, we adopt the procedure described in Ref. Bulava:2021fre ,
where one looks for small enough values of lambda such that fluctuations due
to this unphysical parameter are not relevant compared to the dominating
statistical error. An example of this stability analysis is shown in Fig. 3:
we can see how values for the prediction exist such that the dependence on
unphysical parameters is not significant, and yet the reconstruction is not
lost into statistical noise. Additional tunable parameters, such as $\alpha$,
can be used with the same criterion to better identify this region in
parameter space.

Figure 4: Smearing kernel reconstruction for the vector meson (V) consisting
of fermions transforming in the fundamental representation of $Sp(4)$. The
data analyzed is taken from ensemble M2, and $m_{V}$ has been measured through
the GEVP procedure—see Tab. 6. For both the Gaussian and Cauchy kernels, we
superimpose the target kernel, ${\Delta}_{\sigma}(E-\omega)$, and the
reconstructed one, $\bar{\Delta}_{\sigma}(E-\omega)$. Parameters used are
$t_{\rm max}=N_{t}/2$, $\alpha=0$, and $\sigma=0.40\,m_{\text{V}}$.

While the stability analysis supports the idea that the dependence on the bias
can be absorbed into statistical fluctuation, we take additional steps to
estimate a possible systematic error left in our estimate, by further varying
the algorithmic parameters:

  * •

The first component of the systematic error is

$\sigma_{1,\,\rm
sys}(E)\equiv|\rho_{\lambda_{*}}(E)-\rho_{\lambda_{*}/10}(E)|\,,$ (50)

where $\lambda_{*}$ is determined by the procedure described above, and
$\alpha$ is fixed.

  * •

The second component of the systematic error is

$\sigma_{2,\,\rm
sys}(E)=|\rho_{\lambda_{*},\alpha}(E)-\rho_{\lambda_{*},\alpha^{{}^{\prime}}}(E)|\,,$
(51)

where $\lambda_{*}$ is determined as above and fixed, where we use different
norms, $\alpha\neq\alpha^{{}^{\prime}}$.

Then, these effects are summed up in quadrature with the statistical error to
obtain the error estimate

$\sigma^{2}_{\rho_{\sigma}}=\sigma_{\rm stat}^{2}+\sigma_{1,\,\rm
sys}^{2}+\sigma_{2,\,\rm sys}^{2}\,.$ (52)

###  IV.2 Smearing kernel and smearing radius

Under reasonable conditions, the HLT method is robust, in the sense that
physical results can be extracted effectively. Given the freedom to choose
different smearing kernels, $\Delta_{\sigma}(E,\omega)$, as the arguments
leading to Eq. (42) can be satisfied with a broad class of choices. As a check
of the stability of our results, two particularly interesting choices are the
Gaussian kernels, defined as

$\Delta^{(1)}_{\sigma}(E-\omega)\equiv\frac{1}{Z(\omega)}\exp\left[{\frac{-(E-\omega)^{2}}{2\sigma^{2}}}\right]\;,$
(53)

with $Z(\omega)=\int_{0}^{\infty}dE\,e^{\frac{-(E-\omega)^{2}}{2\sigma^{2}}}$,
and the Cauchy (or Breit-Wigner) kernels, written as

$\Delta^{(2)}_{\sigma}(E-\omega)=\dfrac{\sigma}{(E-\omega)^{2}+\sigma^{2}}\,.$
(54)

Comparing the results obtained with these two choices provides an estimate of
the systematic uncertainty for the fits of spectral densities deriving from
this source. As expected, we anticipate that no relevant dependence (compared
to statistical uncertainties) will be observed in this latter check.

Figure 5: Illustration of kernel reconstruction systematics dependence on
energy. The curves are the target Gaussian kernels,
$\Delta_{\sigma}(E-\omega)$, for $\omega/m_{\mathrm{V}}\approx
1.0,\,1.5,\,1.9$ (top to bottom panels), and $\sigma=0.28\,m_{\mathrm{V}}$.
The points represent the reconstructed kernels,
$\bar{\Delta}_{\sigma}(E-\omega)$. The data is taken from ensemble M2, and
$m_{\mathrm{V}}$ is extracted through the GEVP procedure—see Tab. 6. Other
parameters are $t_{\rm max}=N_{t}/2$, and $\alpha=0$.

An example of a successful reconstruction of such kernels,
$\bar{\Delta}_{\sigma}(E-\omega)$ (Gaussian and Cauchy cases), compared to the
target ones, $\Delta_{\sigma}(E-\omega)$, is shown in Fig. 4. As an additional
precaution, we only consider values of $\lambda$ for which
$A[\vec{g}]/A_{0}<0.1$, which constrains the size of the systematics. The
other parameters are discussed in the caption. Such an outcome corresponds to
what we generally find in our analysis. As expected, given the finiteness of
the lattice and the finite statistics of correlation function measurements,
the quality of the reconstruction deteriorates at the largest energies
considered, as shown in Fig. 5. The goal for a successful reconstruction is to
keep the discrepancy between the targets $\Delta_{\sigma}(E-\omega)$ and the
reconstructed $\bar{\Delta}_{\sigma}(E-\omega)$ as small as possible for each
energy considered.

We find it convenient to express the choice of smearing radius, $\sigma$, in
units of the ground state mass, $m_{0}$, for the mesonic channel of interest,
which can be measured independently, for example by using the GEVP procedure.
The optimal choice of $\sigma$ depends on the quality of the data: the amount
of statistics available for the input two-point correlation functions also
affects the lower bound for the smearing radius $\sigma$. Typically, the
values used for our ensembles are in the range $0.18m_{0}\leq\sigma\leq
0.35m_{0}$. In this interval, the smaller smearing radii are being used to
resolve the details of the spectrum when the discrete energy levels are
tightly packed close together.

Figure 6 displays the reconstruction of a set of synthetic two-point
correlation functions, built to contain a ground state and an excited state.
We have computed $\hat{\rho}_{\sigma}(E)$ using three illustrative values of
the smearing radii ($\sigma=0.1m_{0},\,0.3m_{0},\,0.7m_{0}$). For the smallest
smearing radius, $\sigma=0.1m_{0}$, as one increases the energy considered,
$E/m_{0}$, the smearing kernels appear to be poorly reproduced when
$E/m_{0}\gtrapprox 1.0$, so that excited states are lost below the level of
systematic errors. On the other hand, with the extreme choice
$\sigma=0.7m_{0}$, the two expected peaks merge completely, and the
reconstruction keeps statistical and systematic error moderate. By choosing an
intermediate value, $\sigma=0.3m_{0}$, there is evidence that the
reconstruction resolves well the two peaks, and keeps the systematic and
statistical uncertainties controllable.

Figure 6: Reconstruction of spectral densities from synthetic 2-point
correlation functions, engineered to have ground state mass $am_{0}=0.33$, and
contain one excited state with $m^{*}_{0}=2m_{0}$. The number of
configurations is $N_{\rm conf}=1000$, the lattice time extent is $N_{t}=96$,
and the relative uncertainty is $2\%$. The three panels correspond to Gaussian
kernel reconstructions with $\sigma=0.1m_{0},\,0.3m_{0},\,0.7m_{0}$ (top to
bottom).

###  IV.3 Spectral density fits

In this Section, we describe how smeared spectral densities can be used for
meson spectroscopy. The method involves minimising the correlated functional,
$\chi^{2}$, defined as in Ref. DelDebbio:2022qgu , along the same lines as Eq.
(28), as follows:

$\chi^{2}\equiv\sum_{E,\,E^{{}^{\prime}}}\left(f^{(k)}_{\sigma}(E)-\hat{\rho}_{\sigma}(E)\right)\text{Cov}^{-1}_{E\,E^{{}^{\prime}}}[\rho_{\sigma}]\left(f^{(k)}_{\sigma}(E^{{}^{\prime}})-\hat{\rho}_{\sigma}(E^{{}^{\prime}})\right)\,,$
(55)

where the covariance matrix, $\text{Cov}_{E\,E^{{}^{\prime}}}$, has been
defined in Eq. (39), and the fitting function, $f^{(k)}_{\sigma}(E)$, is
defined as a sum of instances of either the Gaussian ($i=1$, Eq. (53)) or
Cauchy ($i=2$, Eq. (54)) kernel:

$f^{(k)}_{\sigma}(E)=\sum_{n=1}^{k}\mathcal{A}_{n}\,\Delta^{(i)}_{\sigma}(E-E_{n})\,.$
(56)

The fit parameters $E_{n}$ will be identified with the eigenvalues of the
finite-volume Hamiltonian. The freedom in choosing $k$ refers to the fact that
a priori one does not know how many excited states the method can identify;
this is fixed a posteriori.

Figure 7: Spectral density reconstruction in the $\rm v$ channel, in ensemble
M2 ($N_{t}=64$). The black datapoints are the optimal reconstructed spectral
density values $\hat{\rho}_{\sigma}(\omega)$. The yellow line corresponds to
the fitted curve $f^{(2)}_{\sigma}(E)$ in Eq. (56), the red line corresponds
to the deriving ground state curve $\Delta_{\sigma}(E-E_{0})$ from the fit,
the green line is the first excited state $\Delta_{\sigma}(E-E_{1})$. For the
left panel, the reconstruction is done using Gaussian kernels. For the right
panel, it is performed using Cauchy kernels. Both the spectral densities have
smearing radius $\sigma=0.25\,m_{\hbox{\small{v}}}$. The Wuppertal smearing
applied has $N_{\rm source}=0$, $N_{\rm sink}=40$ iterations. For the left
panel: the fit has reduced ${\chi}^{2}/N_{\rm d.o.f.}=1.6$, the energy levels
are $aE_{0}=0.6511(32)$ and $aE_{1}=0.932(11)$, and the associated amplitudes
have $\mathcal{A}_{1}=3.73(23)\times 10^{-6}$ and
$\mathcal{A}_{2}=2.52(52)\times 10^{-6}$, respectively. For the right panel:
${\chi}^{2}/N_{\rm d.o.f.}=1.1$, $aE_{0}=0.6442(37)$ and $aE_{1}=0.929(12)$,
the associated amplitudes are $\mathcal{A}_{1}=2.33(12)\times 10^{-6}$ and
$\mathcal{A}_{2}=1.63(16)\times 10^{-6}$. The value of $m_{\hbox{\small{v}}}$
has been found through the GEVP analysis—see Tab. 6.

As anticipated, we exploit the availability of two different and independent
smearing kernels as an additional sense check for the results of our
spectroscopy study, given that the eigenvalues of the finite-volume
Hamiltonian do not depend on the kernel. For the determination of the $n^{th}$
energy level, we define the systematic error to be

$\sigma_{\rm sys}(aE_{n})=|aE_{n,\,{\rm Gauss}}-aE_{n,\,{\rm Cauchy}}|\,.$
(57)

In Sec. VI, we demonstrate that possible systematics due to the choice of the
kernel never lead to fluctuations outside of statistical errors.

We display in Fig. 7 an example of the comparison between the results of the
reconstruction process performed using Gauss and Cauchy kernels. The plots
show qualitative and quantitative differences in the spectral shapes, due to
the use of different choices of $\Delta_{\sigma}$ functions. Nonetheless, the
position of the corresponding (Gaussian or Cauchy) peaks are compatible within
statistical uncertainties.

Figure 8: Assessment of the fitting systematic error in the ground state
energy level, for ensemble M2. For each channel, the results obtained by
fitting using $k$ and $k+1$ peaks, and Gaussian and Cauchy kernels, are
compared. Further details are presented in Tab. 6. The offset between
datapoints in the same mesonic channel is set for visual clarity. Figure 9:
Assessment of the fitting systematic error in the first excited energy level,
for ensemble M2. For each channel, the results obtained by fitting using $k$
and $k+1$ peaks, and Gaussian and Cauchy kernels, are compared. Further
details are presented in Tab. 7. The offset between datapoints in the same
mesonic channel is set for visual clarity.

Given the finiteness of the smearing radius, and the deterioration of the
reconstruction at high energies, one will only be able to access a certain
number of states in the mesonic spectrum. Even then, fit results can be
contaminated by additional states that are not accounted for in the model. In
order to ensure this does not hinder our prediction for the first $k$ states,
we can repeat the fit by adding a state in the model–performing a
reconstruction with $k+1$ (Gaussian or Cauchy) peaks in Eq. (55), to check the
presence of excited states contamination. If the targetted states are stable
under this procedure, and the chi-square does not improve significantly, we
consider the estimate reliable. As we will show, our measurements appear free
from this effect. An analysis of the systematics due to different kernels, and
contaminations from excited states, is shown in Figs 8 and 9, for the ground
states and first excited states, respectively, in a selection of ensembles and
meson channels.

###  IV.4 On the effect of APE and Wuppertal smearing on the spectral
densities

In analogy with what is done for the correlation functions in Eq. (19), one
can re-express the spectral density

$\rho_{\sigma}(E)=\sum_{n}\dfrac{\langle 0|\mathcal{O}(0)|n\rangle\langle
n|\bar{\mathcal{O}}(0)|0\rangle}{2E_{n}}\,\Delta_{\sigma}\Big{(}E-E_{n}(N_{s}a)\Big{)}\,,$
(58)

using the matrix elements $\langle 0|\mathcal{O}(0)|n\rangle$ and $\langle
n|\bar{\mathcal{O}}(0)|0\rangle$, in agreement with the functional form in Eq.
(56). The matrix elements in Eq. (58) determine the relative contribution of
each energy level to the smeared spectral density.

A major obstacle in performing spectroscopy of lattice gauge theories consists
in the difficulties of disentangling their spectrum, because the energies
contributing to Eq. (58) can be very close to each other. Moreover, certain
states can have large contributions, at the risk of obfuscating others. This
phenomenon may result in considerable discrepancies in the matrix elements
corresponding to different states, $|n\rangle$, posing a challenge to
applications in spectroscopy.

The same potential obstruction would affect the direct study of two-point
correlation functions, as it results in a distortion of the large-$t$ behavior
of two-point correlation function in Eq. (19), and hence the quality of multi-
exponential fits of effective mass plateaux deteriorates. As discussed in Sec.
III, this difficulty is addressed by introducing an appropriately tuned
combination of APE and Wuppertal smearing in the extraction of the correlation
functions. By doing so, one improves the overlap of the states of interest
with the interpolating operator, therefore reducing the importance of other,
undesired states.

Figure 10: Example of the high quality plateau and fit obtained with an appropriate application of APE and Wuppertal smearing. For ensemble M1, the correlation function in the V channel is obtained by applying APE smearing with $N_{\rm APE}=50$ steps (APE smearing step-size $\alpha_{\textrm{APE}}=0.4$), and Wuppertal smearing with $N_{\rm source}=80$ source smearing steps and $N_{\rm sink}=40$ sink smearing steps (Wuppertal smearing step size $\varepsilon_{\rm f}=0.12$). The effective mass is displayed in the left panel. Applying the GEVP procedure, we find $m_{\mathrm{V}}=0.4098(25)$, $m^{*}_{\mathrm{V}}=0.700(26)$. The right panel displays the reconstructed spectral density with Gaussian kernel, for $\sigma=0.30\,m_{V}$, including two states. The black datapoints are the optimal reconstructed spectral density values $\hat{\rho}_{\sigma}(\omega)$. The yellow line corresponds to the fitted curve $f^{(2)}_{\sigma}(E)$ in Eq. (56), the red line corresponds to the deriving ground state curve $\Delta_{\sigma}(E-E_{0})$ from the fit, the green line is the first excited state $\Delta_{\sigma}(E-E_{1})$. The fit results are $\mathcal{A}_{1}=2.89(30)\times 10^{-7}$, $\mathcal{A}_{2}=3.33(37)\times 10^{-7}$, $aE_{0}=0.4139(49)$, and $aE_{1}=0.702(19)$. The reduced chi square is ${\chi}^{2}/N_{\rm d.o.f.}=0.15$. Case | $\alpha_{\textrm{APE}}$ | $\varepsilon_{\rm f}$ | $N_{\rm source}$ | $N_{\rm sink}$ | $\mathcal{A}_{2}/\mathcal{A}_{1}$ | $aE_{0}$ | $aE_{1}$
---|---|---|---|---|---|---|---
A | 0.4 | 0.12 | 80 | 20 | 1.32(19) | 0.4144(50) | 0.692(27)
B | 0.4 | 0.12 | 80 | 40 | 1.15(11) | 0.4139(49) | 0.702(19)
C | 0.4 | 0.12 | 80 | 80 | 0.75(15) | 0.4131(52) | 0.699(22)
D | 0.4 | 0.12 | 40 | 80 | 1.24(18) | 0.4132(43) | 0.694(27)
E | 0.4 | 0.12 | 20 | 80 | 1.80(28) | 0.4148(51) | 0.714(23)
F | 0.4 | 0.24 | 90 | 30 | 1.01(20) | 0.4148(52) | 0.698(19)
G | 0.4 | 0.4 | 170 | 170 | 0.63(11) | 0.4113(82) | 0.717(33)
H | 0.4 | 0.05 | 20 | 20 | 2.28(27) | 0.4136(74) | 0.705(32)
I | 0.0 | 0.12 | 80 | 40 | 1.27(11) | 0.4154(73) | 0.698(32)
Table 3:  Illustration of the effectiveness of APE and Wuppertal smearing in
spectroscopy analysis. For ensemble M1, the correlation function in the V
channel is obtained by applying different levels of sink and source Wuppertal
smearing and APE smearing, and we demonstrate the effect on the output
results, $\mathcal{A}_{2}/\mathcal{A}_{1}$, $aE_{0}$, and $aE_{1}$.

An example of the result of the optimization of the smearing parameters is
shown in Fig. 10. The right panel displays the effect of applying a level of
smearing that optimizes the GEVP extraction of the ground state, chosen so
that the effective mass plateaux are clearly discernible (left panel). The
ground state has comparable amplitude with the first excited state(s). To be
more explicit, in Tab. 3 we report various choices of the smearing parameters
and they affect the output results, for the same channel and ensemble as in
Fig. 10. Cases A to F in Tab. 3 all correspond to reasonable levels of
smearing: the relative amplitudes of the peaks corresponding to the ground and
first excited states are comparable and yield to accurate spectroscopy
results. Conversely, if one applies too much smearing, such as case G in Tab.
3, or too little, as in H, the relative difference between the amplitudes
increases, which results in less precise determinations of the energy levels.
We verified that the amplitudes scale as expected from Eq. (58). By applying
smearing to the operators increases the overlap with the ground state and
decreases the overlap with the excited states. The ratio
$\mathcal{A}_{2}/\mathcal{A}_{1}$ is smaller for larger smearing. For the same
reason, and as expected from Eq. (19), another effect of smearing is the
appearance of longer effective-mass plateaux for the ground state.

Figure 11:  The V channel in ensemble M1, for which the GEVP procedure yields
$m_{\mathrm{V}}=0.4098(25)$ and $m^{*}_{\mathrm{V}}=0.700(26)$. Top panels:
effective mass and spectral density, for $\sigma=0.30\,m_{\mathrm{V}}$, with
APE smearing with $N_{\rm APE}=50$ steps (APE smearing step-size
$\alpha_{\textrm{APE}}=0.4$), and Wuppertal smearing with $N_{\rm source}=20$
source steps and $N_{\rm sink}=20$ sink steps (Wuppertal step size
$\varepsilon_{\rm f}=0.05$). An example of an optimal amount of (APE and
Wuppertal smearing) in this case is shown in Fig. 10. The black datapoints are
the optimal reconstructed spectral density values
$\hat{\rho}_{\sigma}(\omega)$. The yellow line corresponds to the fitted curve
$f^{(2)}_{\sigma}(E)$ in Eq. (56), the red line corresponds to the deriving
ground state curve $\Delta_{\sigma}(E-E_{0})$ from the fit, the green line is
the first excited state $\Delta_{\sigma}(E-E_{1})$. Fit results are
$\mathcal{A}_{1}=1.00(5)\times 10^{-6}$, $\mathcal{A}_{2}=2.28(54)\times
10^{-6}$, $aE_{0}=0.4166(74)$, $aE_{1}=0.705(32)$, and reduced
${\chi}^{2}/N_{\rm d.o.f.}=0.9$. Bottom panels: same as top, except that
Wuppertal smearing has $N_{\rm source}=170$ source steps and $N_{\rm
sink}=170$ sink steps (Wuppertal step size $\varepsilon_{\rm f}=0.90$). Fit
results are $\mathcal{A}_{1}=1.80(26)\times 10^{-8}$,
$\mathcal{A}_{2}=1.13(28)\times 10^{-8}$, $aE_{0}=0.4113(82)$,
$aE_{1}=0.717(33)$, and reduced ${\chi}^{2}/N_{\rm d.o.f.}=0.8$. Figure 12:
The V channel in ensemble M1, for which the GEVP procedure yields
$m_{\mathrm{V}}=0.4098(25)$ and $m^{*}_{\mathrm{V}}=0.700(26)$. With no APE
nor Wuppertal smearing applied, spectral density reconstruction uses
$\sigma=0.60\,m_{\mathrm{V}}$. An example of an optimal amount of (APE and
Wuppertal smearing) in this case is shown in Fig. 10. The same optimal case
without applying APE smearing is shown in Fig. 13. The black datapoints are
the optimal reconstructed spectral density values
$\hat{\rho}_{\sigma}(\omega)$. The yellow line corresponds to the fitted curve
$f^{(2)}_{\sigma}(E)$ in Eq. (56), the red line corresponds to the deriving
ground state curve $\Delta_{\sigma}(E-E_{0})$ from the fit, the green line is
the first excited state $\Delta_{\sigma}(E-E_{1})$. Fit results are
$\mathcal{A}_{1}=1.41(22)\times 10^{-6}$, $\mathcal{A}_{2}=6(1)\times
10^{-6}$, $aE_{0}=0.413(11)$, $aE_{1}=0.839(77)$, and reduced
${\chi}^{2}/N_{\rm d.o.f.}=0.80$. Figure 13:  The V channel in ensemble M1,
for which the GEVP procedure yields $m_{\mathrm{V}}=0.4098(25)$ and
$m^{*}_{\mathrm{V}}=0.700(26)$. No APE smearing has been applied, while
Wuppertal smearing has $N_{\rm source}=80$ source steps and $N_{\rm sink}=40$
sink steps (Wuppertal step size $\varepsilon_{\rm f}=0.12$). An example of an
optimal amount of (APE and Wuppertal smearing) in this case is shown in Fig.
10. The spectral density assumes $\sigma=0.30\,m_{\mathrm{V}}$. The black
datapoints are the optimal reconstructed spectral density values
$\hat{\rho}_{\sigma}(\omega)$. The yellow line corresponds to the fitted curve
$f^{(2)}_{\sigma}(E)$ in Eq. (56), the red line corresponds to the deriving
ground state curve $\Delta_{\sigma}(E-E_{0})$ from the fit, the green line is
the first excited state $\Delta_{\sigma}(E-E_{1})$. The fit results are
$\mathcal{A}_{1}=8.13(32)\times 10^{-13}$, $\mathcal{A}_{2}=1.03(12)\times
10^{-12}$, $aE_{0}=0.4154(73)$, $aE_{1}=0.698(32)$, and reduced
${\chi}^{2}/N_{\rm d.o.f.}=1.3$.

For pedagogical purposes, we find it useful to display, in Fig. 11, also two
cases in which the choice of smearing parameters leads away from optimal
results, in cases G and H of Tab. 3. In producing the top panels, only a tiny
amount of Wuppertal smearing is applied to the two-point correlation
functions, resulting in very short (or absent) effective mass plateaux and
poor resolution of the spectral density. A similar difficulty emerges if one
applies too large amounts of Wuppertal smearing, as depicted in the bottom
panels of Fig. 11: the plateau practically disappears from the effective mass
plot, which appears to be dominated by uncontrolled systematics, and the
resolution of the spectral density deteriorates.

To make the point clearer, in Fig. 12 we depict the case in which no APE nor
Wuppertal smearing has been applied. In this case, the effective-mass plateau
is short and might lead to arguable determinations of effective mass fits.
This is a particularly effective illustration of why the use of point sources
may be problematic. In such a case, to perform a reliable reconstruction it is
necessary to use a larger spectral density smearing radius, $\sigma$, compared
to all the cases considered above. Moreover, the difference in peak heights is
substantial. All these factors lead to a deterioration of the signal.

In Fig. 13, the same reconstruction as in Fig. 10 is considered, but with no
APE smearing applied. By comparing the effective mass plateaux in the left
panels of Figs. 10 and 13 (cases B and I), one sees how APE smearing affects
the effective mass plots and spectroscopy. As APE smearing averages out the
ultraviolet fluctuations in gauge links, without it, the spectral density fits
deteriorate at large values of the energy, as illustrated by the right panels
of the figures. Hence, APE smearing of two-point correlation functions results
in a widening of the energy window available to spectral reconstruction.

##  V Spectral energy density reconstruction and lattice time extent

Figure 14: Examples of reconstructed Gaussian kernels at large energy, for
$\omega/m_{\mathrm{V}}=1.8$, using different lattice temporal extents,
$N_{t}=48$, $N_{t}=64$, and $N_{t}=96$ cases, in ensembles M1, M2, and M3,
respectively. The value of $m_{V}$ is the ground state for the V channel in
ensemble M3, extracted through the GEVP process—see Tab. 9. Figure 15:
Reconstructed spectral density in the PS meson channel, with Gaussian kernel.
Correlation functions are obtained for APE smearing with $N_{\rm APE}=50$
steps and step-size $\alpha_{\textrm{APE}}=0.4$. The ground state mass,
$m_{\rm PS}$, is obtained for the GEVP of the longer time extent—see Tab. 9.
In all the three panels, the black datapoints are the optimal reconstructed
spectral density values $\hat{\rho}_{\sigma}(\omega)$. Top panel: M1 ensemble
($N_{t}=48$), Wuppertal smearing step size $\varepsilon_{\rm f}=0.18$, fit
results $aE_{0}=0.3685(21)$, $aE_{1}=0.649(20)$ and reduced ${\chi}^{2}/N_{\rm
d.o.f}=0.2$. The yellow line corresponds to the fitted curve
$f^{(2)}_{\sigma}(E)$ in Eq. (56), the red line corresponds to the deriving
ground state curve $\Delta_{\sigma}(E-E_{0})$ from the fit, the green line is
the first excited state $\Delta_{\sigma}(E-E_{1})$. Middle panel: M2 ensemble
($N_{t}=64$), Wuppertal smearing step size $\varepsilon_{\rm as}=0.20$, fit
results $aE_{0}=0.3652(18)$, $aE_{1}=0.659(18)$, and reduced
${\chi}^{2}/N_{\rm d.o.f.}=0.9$. The yellow line corresponds to the fitted
curve $f^{(2)}_{\sigma}(E)$ in Eq. (56), the red line corresponds to the
deriving ground state curve $\Delta_{\sigma}(E-E_{0})$ from the fit, the green
line is the first excited state $\Delta_{\sigma}(E-E_{1})$. Bottom panel: M3
ensemble ($N_{t}=96$), Wuppertal smearing step size $\varepsilon_{\rm
as}=0.20$, fit results $aE_{0}=0.3678(10)$, $aE_{1}=0.6693(92)$,
$aE_{2}=0.888(24)$, and reduced ${\chi}^{2}/N_{\rm d.o.f.}=0.5$. The grey line
corresponds to the fitted curve $f^{(3)}_{\sigma}(E)$ in Eq. (56), the red
line corresponds to the deriving ground state curve $\Delta_{\sigma}(E-E_{0})$
from the fit, the green line is the first excited state
$\Delta_{\sigma}(E-E_{1})$, and the yellow line is the second excited state
$\Delta_{\sigma}(E-E_{2})$. Figure 16: Spectroscopy results corresponding to
Fig. 15 for ground and first excited states. The central value of the
corresponding energy levels corresponds to the datapoints, whereas the width
is the uncertainty in their determination. Off-sets between the different
datapoints are kept for visual clarity.

The reconstructed smearing kernel, $\bar{\Delta}_{\sigma}(E,\omega)$, can be
expressed in terms of the function $b(t,E)$—see Eq. (45) in Sec. IV.1. This
finite sum over $t$ has a precision that is limited by the lattice extent, as
the largest possible choice of $t=t_{\rm max}$ is bound by the constraint
$t_{\rm max}\leq N_{t}/2$.

At a vanishing value of the trade-off parameter $\lambda$, this sum is known
to converge to the true kernel in the limit of infinite time extents
Hansen:2019idp . At a non-zero $\lambda$, one expects the reconstruction to
improve, because the basis of functions generating the kernel is larger, and
one can afford smaller values of $\lambda$, thus reducing the bias. To
quantify this effect, we consider lattices with different extents, $N_{t}$, in
the temporal direction, while keeping all other lattice parameters fixed.
Ensembles M1, M2, and M3, have time extents $N_{t}=48,\,64$ and $96$,
respectively, but all other lattice parameters are common. We then compare the
spectral reconstruction, focusing our attention on the region of parameter
space at large energy.

The first result of this comparative study is shown in Fig. 14, which displays
the reconstructed Gaussian kernel, $\bar{\Delta}_{\sigma}(E,\omega)$, with
$\omega/m_{\rm V}=1.8$, for the three aforementioned ensembles. The choice of
$t_{\rm max}$ is the maximum possible value, $t_{\rm max}=N_{t}/2$. The
reconstruction obtained with the shorter available time extent, $N_{t}=48$,
leads to non-negligible deviations from the target kernel,
$\Delta_{\sigma}(E,\omega)$, and hence higher systematic uncertainties, both
in the region away from the maximum of the kernel, but also in the central
region. These indications are compatible with the results shown in Fig. (5).
Conversely, the longest available time extent, $N_{t}=96$, leads to
reconstructed kernels that have visibly smaller deviations from the exact one.

By increasing the time extent of the lattice one can perform accurate spectral
density fits and spectroscopy, reaching progressively higher energies. Hence,
it becomes possible to extend the number of measurable excited states. We
illustrate this phenomenon in Fig. 15, for the $({\rm PS}$) meson channel and
the three time extents $N_{t}=48,64$ and $96$. A stabilization of the spectral
reconstruction with respect to the algorithmic parameters, described in Sec.
IV.1, can be obtained at higher energies, which become accessible by extending
the time extent. As a consequence, we are able to obtain a second peak
structure, which encodes information about excited states. Moreover, as will
be discussed in Sec. VI, the positions of the first two peaks agree, within
statistical errors, and do not show a significant dependence on $N_{t}$. As a
consequence, we can expect our fits to the smeared spectral densities to
improve, possibly capturing more states. We illustrate this phenomenon in Fig.
16. The three ensembles considered, M1, M2, and M3, respectively, yield the
difference $[aE_{n}-\sigma_{aE_{n}},aE_{n}+\sigma_{aE_{n}}]$, where
$\sigma_{aE_{n}}$ is the uncertainty in the determination of the $n^{th}$
energy level.

Figure 17: Plateau in the effective mass, $am_{\rm eff}$, for the PS channel
in ensemble M2 ($N_{t}=64$). The plateaux are identified via the GEVP process,
by using a basis of nine correlators as shown in Eq. (36), obtained by varying
$N=0$, $40$, $80$ sink and source Wuppertal smearing steps (keeping fixed
$\varepsilon_{\rm as}=0.10$, $\alpha_{\textrm{APE}}=0.4$, and $N_{\rm
APE}=50$). Figure 18: Plateaux in effective mass, $am_{\rm eff}$, for the
combined V and T channels, in ensemble M2 ($N_{t}=64$). The plateaux are
identified via the GEVP process, by using a basis of thirty-six correlators as
shown in Eq. (37): besides varying the smearing levels, $N=0$, $40$, $80$, of
both sink and source Wuppertal smearing (and keeping fixed $\varepsilon_{\rm
f}=0.2$, $\alpha_{\textrm{APE}}=0.4$, and $N_{\rm APE}=50$), we combine the
channels with the cross channels. In the case of degeneracy in the plateau
fitting (within statistical uncertainties), we choose the energy level with a
smaller statistical error as a representative value for the considered energy
level. This is the case in this example, where $aE_{1}$ has been chosen over
$aE_{2}$ (table 7) as an estimation for the first excited state.

##  VI Numerical results and comparisons with GEVP

Figure 19: The spectrum, obtained through spectral density fits, for the
pseudoscalar operators made of fundamental representation fermions for all the
ensembles studied. Tower of masses, $\hat{m}\equiv w_{0}\cdot m$ (in Wilson
flow units), corresponding to ground, first and (possibly) second excited
states is shown. The vertical midpoint of each colorblock is the numerical
result, whereas the heights are comprehensive of statistical and systematic
errors summed in quadrature. Horizontal offsets distinguish different
ensembles. Different shadings of the same colors differentiate ensembles that
differ in time extents ($N_{t}=48,\,64,\,96$ for ensembles M1, M2 and M3),
whereas hatched patterns, with no fill, are used to indicate ensembles that
differ also in bare parameters (ensembles M4 and M5). As expected,
progressively smaller uncertainties are obtained in the results for M1, M2 and
M3.  Figure 20: Meson spectrum in all the ensembles (summarised in Tab. 1),
for both fundamental and antisymmetric representation fermions, found through
spectral densities fitting analysis. For each channel, a tower of masses,
$\hat{m}\equiv w_{0}\cdot m$ (in Wilson flow units), corresponding to ground,
first and (where available) second excited states is shown. The vertical
midpoint of each colorblock is the numerical result, whereas the heights are
comprehensive of statistical and systematic errors summed in quadrature.
Horizontal offsets distinguish different ensembles. Different shadings of the
same colors differentiate ensembles that differ in time extents
($N_{t}=48,\,64$ and $96$ for ensembles M1, M2 and M3), whereas no filling
color and patterns are used to indicate ensembles that differ also in bare
parameters (ensembles M4 and M5). Progressively smaller uncertainties are
obtained in the results for M1, M2 and M3. Six colors distinguish the
different meson channels, and the colors match where two different
representations–fundamental and antisymmetric–are used to study meson
operators built with the same gamma-matrix structure.

In this section, we display and discuss our spectroscopy results, obtained
with the HLT method. Firstly, we showcase how our results match the
expectations from the GEVP ones. Figures 17 and 18 display two representative
examples of GEVP computations. The former case is obtained as described in Eq.
(36), by using a variational basis of nine elements, whereas the latter
combines the V and T channels, according to Eq. (37). As shown in the figures,
typically the signal-to-noise ratio is good enough to find both ground state
and first excited state even for the smaller variational basis. By adding the
cross channel one gains access to the second excited state, $E_{2}$.

$C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{0}$ $k$-G |  $aE_{0}$ ($k+1$)-G |  $aE_{0}$ $k$-C |  $aE_{0}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 80 | 40 | 0.3685(21) | 0.3698(28) | 0.3675(32) | 0.3668(38) | 0.3678(17) | 0.33 | 0.32
V | 2 | 80 | 40 | 0.4099(59) | 0.4132(79) | 0.4083(25) | 0.4075(28) | 0.4098(25) | 0.30 | 0.22
T | 2 | 80 | 40 | 0.4086(27) | 0.4093(59) | 0.4032(55) | 0.4094(35) | 0.4098(25) | 0.30 | 0.30
AV | 2 | 80 | 40 | 0.5543(73) | 0.5494(92) | 0.5469(85) | 0.5553(94) | 0.5485(81) | 0.20 | 0.18
AT | 2 | 80 | 40 | 0.5519(86) | 0.5494(78) | 0.5458(75) | 0.5453(89) | 0.5514(73) | 0.18 | 0.20
S | 2 | 80 | 40 | 0.5286(85) | 0.5273(99) | 0.5251(69) | 0.5293(81) | 0.5241(64) | 0.20 | 0.20
ps | 2 | 80 | 40 | 0.5998(32) | 0.6008(28) | 0.5998(35) | 0.6009(32) | 0.60161(91) | 0.18 | 0.20
v | 2 | 80 | 40 | 0.6458(44) | 0.6451(39) | 0.6546(31) | 0.6493(26) | 0.6503(13) | 0.20 | 0.26
t | 2 | 80 | 40 | 0.6475(39) | 0.6478(32) | 0.6553(35) | 0.6583(69) | 0.6503(13) | 0.30 | 0.26
av | 2 | 80 | 40 | 0.8376(72) | 0.8369(94) | 0.8334(71) | 0.8289(63) | 0.8299(81) | 0.20 | 0.18
at | 2 | 80 | 40 | 0.8517(85) | 0.8566(99) | 0.8499(72) | 0.8534(91) | 0.8408(87) | 0.18 | 0.18
s | 2 | 80 | 40 | 0.7993(85) | 0.7997(84) | 0.7954(78) | 0.7930(87) | 0.7957(83) | 0.18 | 0.20
Table 4:  Numerical results for the ground state mass in ensemble M1. $k$-G
stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for
$k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates
the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis,
$\sigma_{G}$ is the energy smearing radius used for the Gaussian fits,
$\sigma_{C}$ for the Cauchy fit.

In Fig. 19, the spectrum for the fundamental pseudoscalar meson is shown for
all the ensembles of Tab. 1. The figure shows masses normalized in Wilson flow
units, $\hat{m}\equiv w_{0}\cdot m$. Different shadings and small horizontal
offsets have been applied to distinguish the ensembles M1, M2 and M3, whereas
larger offsets and patterns display the ensembles M4 and M5. The vertical
extent of each colorblock represents the sum in quadrature of statistical
errors and systematic effects. The latter ones are given by the excited states
contaminations artifacts and by evaluating methodological effects using
different smearing kernels, as the maximal difference between lattice results
in the $k$ and $k+1$ peaks Gauss and Cauchy fits, $k$-G, $(k+1)$-G and $k$-C,
$(k+1)$-C. In Fig. 20, the same logic is applied to the whole meson spectrum
of the theory, and six different colors have been used for different meson
channels, and the matching colors have been used for the same meson channel in
different fermion representations–fundamental and antisymmetric.

The numerical results for the ground states, first excited states, and (where
available) second excited states obtained are reported in Tabs. 4 to 17, for
all the ensembles M1 to M5. We write the masses in units of the lattice
spacing, $a$, but these can be converted to Wilson flow units with the
numerical results reported in Tab. 1. The tables display the results for all
the twelve meson channels of interest, together with details of the smearing
level, and the number of functions, $k$, included in Eq. (56). We tabulate the
results obtained with five alternative analyses of the same correlation
functions: two measurements obtained from spectral density reconstruction with
Gaussian kernels, with $k$ or $k+1$ functions, two measurements obtained with
Cauchy kernels, with $k$ or $k+1$ functions, and one obtained with the GEVP
analysis. We also report the value of the smearing radius in the Gaussian and
Cauchy cases, respectively.

$C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{1}$ $k$-G |  $aE_{1}$ ($k+1$)-G |  $aE_{1}$ $k$-C |  $aE_{1}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 80 | 40 | 0.649(20) | 0.651(26) | 0.652(29) | 0.649(25) | 0.661(25) | 0.33 | 0.32
V | 2 | 80 | 40 | 0.702(19) | 0.697(28) | 0.698(19) | 0.691(11) | 0.700(26) | 0.30 | 0.22
T | 2 | 80 | 40 | 0.706(12) | 0.695(22) | 0.692(21) | 0.697(32) | 0.700(26) | 0.30 | 0.30
AV | 2 | 80 | 40 | 0.760(33) | 0.752(32) | 0.725(23) | 0.721(32) | 0.743(44) | 0.20 | 0.18
AT | 2 | 80 | 40 | 0.789(43) | 0.791(41) | 0.765(31) | 0.785(18) | 0.768(47) | 0.18 | 0.20
S | 2 | 80 | 40 | 0.751(22) | 0.745(28) | 0.754(25) | 0.741(19) | 0.748(14) | 0.20 | 0.20
ps | 2 | 80 | 40 | 0.881(21) | 0.877(31) | 0.864(31) | 0.879(36) | 0.891(19) | 0.18 | 0.20
v | 2 | 80 | 40 | 0.930(22) | 0.924(18) | 0.941(23) | 0.931(21) | 0.955(11) | 0.20 | 0.26
t | 2 | 80 | 40 | 0.959(12) | 0.967(19) | 0.948(22) | 0.968(21) | 0.955(11) | 0.30 | 0.26
av | 2 | 80 | 40 | 1.041(32) | 1.047(34) | 1.027(41) | 1.019(44) | 1.063(42) | 0.20 | 0.18
at | 2 | 80 | 40 | 1.063(34) | 1.068(36) | 1.060(31) | 1.055(35) | - | 0.18 | 0.18
s | 2 | 80 | 40 | 1.037(21) | 1.038(22) | 1.023(19) | 1.038(27) | 1.052(19) | 0.18 | 0.20
Table 5:  Numerical results for the first excited state mass in ensemble M1. $k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for $k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis, $\sigma_{G}$ is the energy smearing radius used for the Gaussian fits, $\sigma_{C}$ for the Cauchy fit. $C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{0}$ $k$-G |  $aE_{0}$ ($k+1$)-G |  $aE_{0}$ $k$-C |  $aE_{0}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 80 | 40 | 0.3632(32) | 0.3652(11) | 0.3643(12) | 0.3623(26) | 0.3656(12) | 0.35 | 0.30
V | 3 | 0 | 40 | 0.4080(30) | 0.4035(32) | 0.4041(22) | 0.4049(19) | 0.4054(19) | 0.28 | 0.33
T | 3 | 0 | 40 | 0.4023(34) | 0.4042(35) | 0.4054(23) | 0.4044(29) | 0.4054(19) | 0.23 | 0.23
AV | 2 | 80 | 40 | 0.5483(82) | 0.5463(84) | 0.5454(91) | 0.5463(89) | 0.5423(90) | 0.30 | 0.18
AT | 2 | 80 | 40 | 0.5523(73) | 0.5512(85) | 0.5475(65) | 0.5427(69) | 0.5477(84) | 0.30 | 0.20
S | 2 | 0 | 40 | 0.5174(76) | 0.5200(74) | 0.5164(94) | 0.5156(98) | 0.5222(78) | 0.30 | 0.20
ps | 3 | 0 | 40 | 0.6020(10) | 0.6048(10) | 0.6068(18) | 0.6025(14) | 0.6007(11) | 0.18 | 0.18
v | 2 | 40 | 80 | 0.6511(32) | 0.6522(42) | 0.6442(37) | 0.6452(28) | 0.6473(12) | 0.20 | 0.20
t | 2 | 40 | 40 | 0.6442(37) | 0.6469(23) | 0.6501(32) | 0.6498(28) | 0.6473(12) | 0.20 | 0.20
av | 3 | 0 | 40 | 0.8213(82) | 0.8187(75) | 0.8217(91) | 0.829(10) | 0.821(11) | 0.18 | 0.20
at | 2 | 0 | 40 | 0.833(12) | 0.836(10) | 0.844(16) | 0.837(18) | 0.834(15) | 0.18 | 0.23
s | 3 | 80 | 40 | 0.7864(87) | 0.7883(94) | 0.7873(99) | 0.7912(91) | 0.7820(96) | 0.18 | 0.18
Table 6:  Numerical results for the ground state mass in ensemble M2. $k$-G
stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for
$k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates
the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis,
$\sigma_{G}$ is the energy smearing radius used for the Gaussian fits,
$\sigma_{C}$ for the Cauchy fit.

As described in Sec. V, we expect to find smaller uncertainties, and achieve
more precise estimates of the energy levels, $aE_{n}$, by considering longer
time extents, $N_{t}$. This is confirmed by inspecting the results in
ensembles M1, M2, and M3, which are characterized by the same lattice bare
parameters and lattice spatial extent, $N_{s}$, while only the time extent,
$N_{t}$, is varied. This is also shown in Fig. 19, where the pseudoscalar
fundamental channel is shown as a representative case of the data shown in
Fig. 20: the numerical results for the masses concerning ensembles M1, M2, M3
progressively present smaller uncertainties, as expected, and their values are
compatible. We consider and show also the cases where the bare fundamental
fermion masses, $am_{0}^{\rm f}$, are varied, in ensembles M4 and M5, for
intermediate lattice time extent, $N_{t}=64$.

$C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{1}$ $k$-G |  $aE_{1}$ ($k+1$)-G |  $aE_{1}$ $k$-C |  $aE_{1}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 80 | 40 | 0.659(23) | 0.664(19) | 0.670(10) | 0.679(13) | 0.688(10) | 0.35 | 0.30
V | 3 | 0 | 40 | 0.677(11) | 0.6721(93) | 0.669(12) | 0.668(11) | 0.6809(65) | 0.28 | 0.33
T | 3 | 0 | 40 | 0.677(10) | 0.671(12) | 0.6784(98) | 0.678(10) | 0.6809(65) | 0.23 | 0.23
AV | 2 | 80 | 40 | 0.793(25) | 0.769(24) | 0.766(22) | 0.741(24) | 0.772(30) | 0.30 | 0.18
AT | 2 | 80 | 40 | 0.803(32) | 0.799(35) | 0.823(29) | 0.805(28) | 0.789(39) | 0.30 | 0.20
S | 2 | 0 | 40 | 0.750(27) | 0.753(20) | 0.760(20) | 0.762(24) | 0.783(25) | 0.30 | 0.20
ps | 3 | 0 | 40 | 0.869(18) | 0.873(19) | 0.865(18) | 0.865(15) | 0.895(16) | 0.18 | 0.18
v | 2 | 40 | 80 | 0.932(11) | 0.928(11) | 0.929(12) | 0.9362(82) | 0.9392(67) | 0.20 | 0.20
t | 2 | 40 | 40 | 0.928(21) | 0.9342(92) | 0.9431(64) | 0.9398(72) | 0.9392(67) | 0.20 | 0.20
av | 3 | 0 | 40 | 1.069(17) | 1.072(16) | 1.080(15) | 1.082(19) | 1.087(14) | 0.18 | 0.20
at | 2 | 0 | 40 | 1.101(20) | 1.102(21) | 1.103(22) | 1.087(23) | 1.076(28) | 0.18 | 0.23
s | 3 | 80 | 40 | 1.037(21) | 1.036(22) | 1.057(11) | 1.060(15) | 1.053(12) | 0.18 | 0.18
Table 7:  Numerical results for the first excited state mass in ensemble M2. $k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for $k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis, $\sigma_{G}$ is the energy smearing radius used for the Gaussian fits, $\sigma_{C}$ for the Cauchy fit. $C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{2}$ $k$-G |  $aE_{2}$ ($k+1$)-G |  $aE_{2}$ $k$-C |  $aE_{2}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 80 | 40 | - | - | - | - | - | 0.35 | 0.30
V | 3 | 0 | 40 | 0.924(18) | 0.909(20) | 0.905(19) | 0.907(20) | 0.913(15) | 0.28 | 0.33
T | 3 | 0 | 40 | 0.878(32) | 0.876(30) | 0.903(17) | 0.900(14) | 0.913(15) | 0.23 | 0.23
AV | 2 | 80 | 40 | - | - | - | - | - | 0.30 | 0.18
AT | 2 | 80 | 40 | - | - | - | - | - | 0.30 | 0.20
S | 2 | 0 | 40 | - | - | - | - | - | 0.30 | 0.20
ps | 3 | 0 | 40 | 1.135(42) | 1.142(43) | 1.132(40) | 1.130(36) | - | 0.18 | 0.18
v | 2 | 40 | 80 | - | - | - | - | 1.028(10) | 0.20 | 0.20
t | 2 | 40 | 40 | - | - | - | - | 1.028(10) | 0.20 | 0.20
av | 3 | 0 | 40 | 1.348(35) | 1.335(34) | 1.337(31) | 1.350(40) | - | 0.18 | 0.20
at | 2 | 0 | 40 | - | - | - | - | - | 0.18 | 0.23
s | 3 | 80 | 40 | 1.389(38) | 1.404(42) | 1.397(37) | 1.394(39) | - | 0.18 | 0.18
Table 8:  Numerical results for the second excited state mass in ensemble M2.
$k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for
$k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates
the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis,
$\sigma_{G}$ is the energy smearing radius used for the Gaussian fits,
$\sigma_{C}$ for the Cauchy fit.

For each energy level, $aE_{n}$, we report the fits obtained with $k$ and
$k+1$ peaks, both for Gaussian and Cauchy kernels, used in Eq. (55). If
contaminations from further excited states are negligible, one expects
compatible results, for both kernels, in going from $k$ to $k+1$; this is
confirmed by the numerical results tabulated. The number of peaks optimizes
the fitting procedure, $k$, is reported in all cases. The number of iterations
for Wuppertal operator smearing, at source and sink, $N_{\rm source}$, and
$N_{\rm sink}$, used for performing the spectral density fits, is reported in
the tables as well.

$C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{0}$ $k$-G |  $aE_{0}$ ($k+1$)-G |  $aE_{0}$ $k$-C |  $aE_{0}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 40 | 0 | 0.3678(10) | 0.36592(80) | 0.36732(92) | 0.36693(79) | 0.36658(81) | 0.30 | 0.27
V | 3 | 40 | 40 | 0.4104(13) | 0.4102(17) | 0.4078(17) | 0.4090(11) | 0.4083(12) | 0.28 | 0.25
T | 3 | 0 | 40 | 0.4072(21) | 0.4074(19) | 0.4081(12) | 0.4091(15) | 0.4083(12) | 0.33 | 0.28
AV | 2 | 40 | 0 | 0.5364(53) | 0.5375(54) | 0.5385(47) | 0.5362(72) | 0.5352(63) | 0.28 | 0.32
AT | 2 | 40 | 40 | 0.5458(55) | 0.5457(65) | 0.5502(54) | 0.5509(49) | 0.5486(47) | 0.30 | 0.18
S | 2 | 40 | 40 | 0.5151(33) | 0.5129(38) | 0.5178(53) | 0.5168(28) | 0.5165(45) | 0.30 | 0.24
ps | 3 | 0 | 40 | 0.60184(89) | 0.60141(81) | 0.60182(62) | 0.60182(51) | 0.60131(58) | 0.23 | 0.22
v | 2 | 0 | 40 | 0.6479(21) | 0.6490(15) | 0.6498(19) | 0.6501(18) | 0.6492(15) | 0.24 | 0.25
t | 2 | 0 | 40 | 0.6482(15) | 0.6480(19) | 0.6518(21) | 0.6519(22) | 0.6492(15) | 0.28 | 0.28
av | 3 | 0 | 40 | 0.8362(58) | 0.8359(61) | 0.8349(48) | 0.8372(59) | 0.8349(65) | 0.18 | 0.25
at | 2 | 0 | 40 | 0.8452(63) | 0.8402(79) | 0.8392(82) | 0.8399(72) | 0.8430(74) | 0.25 | 0.25
s | 3 | 0 | 40 | 0.797(12) | 0.801(14) | 0.798(11) | 0.795(11) | 0.788(12) | 0.23 | 0.23
Table 9:  Numerical results for the ground state mass in ensemble M3. $k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for $k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis, $\sigma_{G}$ is the energy smearing radius used for the Gaussian fits, $\sigma_{C}$ for the Cauchy fit. $C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{1}$ $k$-G |  $aE_{1}$ ($k+1$)-G |  $aE_{1}$ $k$-C |  $aE_{1}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 40 | 0 | 0.6693(92) | 0.6673(96) | 0.6787(84) | 0.6799(81) | 0.6852(87) | 0.30 | 0.27
V | 3 | 40 | 40 | 0.6856(91) | 0.6758(99) | 0.674(12) | 0.678(14) | 0.6819(98) | 0.28 | 0.25
T | 3 | 0 | 40 | 0.688(11) | 0.6848(97) | 0.690(12) | 0.691(13) | 0.6819(98) | 0.33 | 0.28
AV | 2 | 40 | 0 | 0.792(17) | 0.788(16) | 0.794(18) | 0.795(17) | 0.784(12) | 0.28 | 0.32
AT | 2 | 40 | 40 | 0.784(21) | 0.779(23) | 0.778(26) | 0.776(28) | 0.796(23) | 0.30 | 0.18
S | 2 | 40 | 40 | 0.779(22) | 0.785(31) | 0.779(29) | 0.779(39) | 0.793(33) | 0.30 | 0.24
ps | 3 | 0 | 40 | 0.907(12) | 0.906(14) | 0.908(13) | 0.907(15) | 0.9117(94) | 0.23 | 0.22
v | 2 | 0 | 40 | 0.9331(91) | 0.929(13) | 0.9333(59) | 0.9328(74) | 0.9379(69) | 0.24 | 0.25
t | 2 | 0 | 40 | 0.9331(84) | 0.9318(91) | 0.9428(82) | 0.9418(69) | 0.9379(69) | 0.28 | 0.28
av | 3 | 0 | 40 | 1.103(14) | 1.099(14) | 1.099(15) | 1.098(16) | 1.110(15) | 0.18 | 0.25
at | 2 | 0 | 40 | 1.113(22) | 1.113(20) | 1.119(17) | 1.119(17) | 1.115(20) | 0.25 | 0.25
s | 3 | 0 | 40 | 1.032(19) | 1.028(19) | 1.042(32) | 1.043(19) | 1.045(46) | 0.23 | 0.23
Table 10:  Numerical results for the first excited state mass in ensemble M3.
$k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for
$k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates
the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis,
$\sigma_{G}$ is the energy smearing radius used for the Gaussian fits,
$\sigma_{C}$ for the Cauchy fit.

The usage of multiple smearing kernels is an additional safety check against
systematic effects that can potentially be overlooked in the reconstruction.
Where the latter are absent, the spectroscopy should be untainted when
adopting different kernels. This is confirmed by our numerical results. The
fitted energy levels, $aE_{n}$, are compatible with one another, within
statistical uncertainties. In summary, as anticipated in Sec. IV.3, and shown
in an example in Figs. 8 and 9, all the known sources of systematic
uncertainty in the spectral density reconstruction appear to be under control.

$C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{2}$ $k$-G |  $aE_{2}$ ($k+1$)-G |  $aE_{2}$ $k$-C |  $aE_{2}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 40 | 0 | 0.888(24) | 0.878(21) | 0.861(21) | 0.856(24) | - | 0.30 | 0.27
V | 3 | 40 | 40 | 0.924(19) | 0.909(21) | 0.923(17) | 0.911(18) | 0.914(18) | 0.28 | 0.25
T | 3 | 0 | 40 | 0.902(17) | 0.900(19) | 0.921(21) | 0.919(19) | 0.914(18) | 0.33 | 0.28
AV | 2 | 40 | 0 | 1.141(29) | 1.112(32) | 1.132(35) | 1.139(32) | - | 0.28 | 0.32
AT | 2 | 40 | 40 | 1.051(37) | 1.044(34) | 1.021(31) | 1.019(35) | - | 0.30 | 0.18
S | 2 | 40 | 40 | - | - | - | - | - | 0.30 | 0.24
ps | 3 | 0 | 40 | 1.155(38) | 1.162(41) | 1.157(37) | 1.171(39) | - | 0.23 | 0.22
v | 2 | 0 | 40 | - | - | - | - | 1.0191(67) | 0.24 | 0.25
t | 2 | 0 | 40 | - | - | - | - | 1.0191(67) | 0.28 | 0.28
av | 3 | 0 | 40 | 1.374(31) | 1.372(31) | 1.382(27) | 1.376(27) | - | 0.18 | 0.25
at | 2 | 0 | 40 | - | - | - | - | - | 0.25 | 0.25
s | 3 | 0 | 40 | 1.341(31) | 1.339(31) | 1.354(31) | 1.352(32) | - | 0.23 | 0.23
Table 11:  Numerical results for the second excited state mass in ensemble M3. $k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for $k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis, $\sigma_{G}$ is the energy smearing radius used for the Gaussian fits, $\sigma_{C}$ for the Cauchy fit. $C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{0}$ $k$-G |  $aE_{0}$ ($k+1$)-G |  $aE_{0}$ $k$-C |  $aE_{0}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 0 | 40 | 0.4080(30) | 0.4083(34) | 0.4087(23) | 0.4088(22) | 0.4096(11) | 0.30 | 0.30
V | 3 | 40 | 40 | 0.4467(29) | 0.4469(31) | 0.4460(32) | 0.4459(33) | 0.4484(16) | 0.30 | 0.27
T | 3 | 0 | 40 | 0.4473(21) | 0.4475(20) | 0.4488(19) | 0.4491(21) | 0.4484(16) | 0.25 | 0.25
AV | 2 | 80 | 40 | 0.5974(84) | 0.5977(86) | 0.5974(84) | 0.5997(84) | 0.6016(88) | 0.25 | 0.25
AT | 2 | 80 | 40 | 0.6028(89) | 0.6021(85) | 0.6188(84) | 0.6138(83) | 0.6141(95) | 0.25 | 0.25
S | 2 | 0 | 40 | 0.591(10) | 0.590(10) | 0.577(12) | 0.574(13) | 0.580(10) | 0.25 | 0.30
ps | 2 | 0 | 40 | 0.6291(12) | 0.6293(13) | 0.6283(11) | 0.6284(14) | 0.62809(94) | 0.24 | 0.24
v | 2 | 0 | 40 | 0.6692(32) | 0.6690(34) | 0.6722(32) | 0.6725(33) | 0.6716(15) | 0.23 | 0.20
t | 2 | 0 | 40 | 0.6692(30) | 0.6695(33) | 0.6689(32) | 0.6688(31) | 0.6716(15) | 0.23 | 0.25
av | 2 | 40 | 40 | 0.8747(54) | 0.8748(56) | 0.8742(47) | 0.8789(54) | 0.8793(55) | 0.20 | 0.20
at | 2 | 80 | 40 | 0.8745(81) | 0.8765(81) | 0.8765(51) | 0.8755(61) | 0.8797(66) | 0.25 | 0.25
s | 2 | 40 | 40 | 0.8482(99) | 0.8483(96) | 0.8582(87) | 0.8572(82) | 0.8540(90) | 0.24 | 0.20
Table 12:  Numerical results for the ground state mass in ensemble M4. $k$-G
stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for
$k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates
the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis,
$\sigma_{G}$ is the energy smearing radius used for the Gaussian fits,
$\sigma_{C}$ for the Cauchy fit.

Comparison between results obtained by using spectral densities and GEVP show
agreement in correspondence of all mesonic channels of the spectrum of our
theory. The robustness of the results is reassuring and in line with one of
the purposes of this paper, which is to show that measurements obtained with
the former, novel, method match those of the latter, well established, one.
Focussing on the spectrum, in our studies the ground state predictions from
spectral densities tend to have larger uncertainties. This observation agrees
with the findings of Ref. DelDebbio:2022qgu .

$C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{1}$ $k$-G |  $aE_{1}$ ($k+1$)-G |  $aE_{1}$ $k$-C |  $aE_{1}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 0 | 40 | 0.7402(85) | 0.7372(86) | 0.7352(76) | 0.7352(72) | 0.7348(86) | 0.30 | 0.30
V | 3 | 40 | 40 | 0.784(15) | 0.782(14) | 0.764(12) | 0.768(14) | 0.760(15) | 0.30 | 0.27
T | 3 | 0 | 40 | 0.788(24) | 0.789(24) | 0.754(17) | 0.760(18) | 0.760(15) | 0.25 | 0.25
AV | 2 | 80 | 40 | 0.823(23) | 0.820(21) | 0.833(28) | 0.831(27) | 0.809(16) | 0.25 | 0.25
AT | 2 | 80 | 40 | 0.873(13) | 0.869(13) | 0.854(13) | 0.864(13) | 0.857(14) | 0.25 | 0.25
S | 2 | 0 | 40 | 0.765(22) | 0.766(23) | 0.755(18) | 0.757(19) | 0.752(18) | 0.25 | 0.30
ps | 2 | 0 | 40 | 0.923(12) | 0.929(13) | 0.939(11) | 0.923(14) | 0.933(13) | 0.25 | 0.24
v | 2 | 0 | 40 | 0.948(20) | 0.945(20) | 0.944(23) | 0.943(22) | 0.961(19) | 0.23 | 0.20
t | 2 | 0 | 40 | 0.953(25) | 0.951(24) | 0.957(22) | 0.952(21) | 0.961(19) | 0.23 | 0.25
av | 2 | 40 | 40 | 1.120(13) | 1.121(12) | 1.130(15) | 1.123(14) | 1.119(14) | 0.20 | 0.20
at | 2 | 80 | 40 | 1.104(20) | 1.106(19) | 1.146(19) | 1.136(22) | 1.128(20) | 0.25 | 0.25
s | 2 | 40 | 40 | 1.102(19) | 1.104(21) | 1.119(17) | 1.108(17) | 1.110(18) | 0.24 | 0.20
Table 13:  Numerical results for the first excited state mass in ensemble M4. $k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for $k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis, $\sigma_{G}$ is the energy smearing radius used for the Gaussian fits, $\sigma_{C}$ for the Cauchy fit. $C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{2}$ $k$-G |  $aE_{2}$ ($k+1$)-G |  $aE_{2}$ $k$-C |  $aE_{2}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 0 | 40 | - | - | - | - | - | 0.30 | 0.30
V | 3 | 40 | 40 | 0.948(32) | 0.952(31) | 0.953(30) | 0.958(33) | - | 0.30 | 0.27
T | 3 | 0 | 40 | 0.943(30) | 0.948(30) | 0.968(30) | 0.962(30) | - | 0.25 | 0.25
AV | 2 | 80 | 40 | - | - | - | - | - | 0.25 | 0.25
AT | 2 | 80 | 40 | - | - | - | - | - | 0.25 | 0.25
S | 2 | 0 | 40 | - | - | - | - | - | 0.30 | 0.30
ps | 2 | 0 | 40 | - | - | - | - | - | 0.25 | 0.24
v | 2 | 0 | 40 | - | - | - | - | - | 0.23 | 0.20
t | 2 | 40 | 40 | - | - | - | - | - | 0.23 | 0.25
av | 2 | 40 | 40 | - | - | - | - | - | 0.20 | 0.20
at | 2 | 80 | 40 | - | - | - | - | - | 0.25 | 0.25
s | 2 | 40 | 40 | - | - | - | - | - | 0.24 | 0.20
Table 14:  Numerical results for the second excited state mass in ensemble M4.
$k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for
$k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates
the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis,
$\sigma_{G}$ is the energy smearing radius used for the Gaussian fits,
$\sigma_{C}$ for the Cauchy fit.

We also find a general trend towards improvement in the extraction of excited
state masses, for which the quality of results is competitive with the GEVP
ones. This is noteworthy, for the quantity of information given as input for
spectral density measurements is smaller than what used in the GEVP analysis.
The simple case of Eq. (36) makes use of nine two-point correlation functions
in the GEVP analysis, and the cross-channel in Eq. (37) uses thirty-six of
them; by comparison, in both cases, the spectral density method uses only a
single measurement. This finding suggests an opportunity for possible gains in
efficiency in future large-scale studies of multiple excited states.

$C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{0}$ $k$-G |  $aE_{0}$ ($k+1$)-G |  $aE_{0}$ $k$-C |  $aE_{0}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 40 | 40 | 0.3113(12) | 0.3103(13) | 0.3110(10) | 0.3115(16) | 0.31025(63) | 0.25 | 0.25
V | 3 | 40 | 40 | 0.3513(21) | 0.3492(20) | 0.3499(24) | 0.3503(23) | 0.3505(12) | 0.30 | 0.30
T | 3 | 40 | 40 | 0.3510(24) | 0.3462(23) | 0.3469(28) | 0.3513(27) | 0.3505(12) | 0.30 | 0.30
AV | 2 | 40 | 40 | 0.5132(31) | 0.5123(30) | 0.5130(29) | 0.5121(28) | 0.5111(29) | 0.25 | 0.20
AT | 2 | 40 | 40 | 0.5201(42) | 0.5193(32) | 0.5221(52) | 0.5205(32) | 0.5191(33) | 0.20 | 0.20
S | 2 | 40 | 40 | 0.4878(35) | 0.4868(31) | 0.4890(30) | 0.4899(37) | 0.4888(29) | 0.25 | 0.25
ps | 2 | 40 | 40 | 0.57823(39) | 0.57833(43) | 0.57805(41) | 0.57813(46) | 0.57853(41) | 0.20 | 0.20
v | 2 | 40 | 40 | 0.6212(32) | 0.6221(32) | 0.6202(34) | 0.6200(32) | 0.6212(10) | 0.20 | 0.25
t | 2 | 40 | 40 | 0.6232(30) | 0.6213(36) | 0.6213(31) | 0.6212(37) | 0.6212(10) | 0.20 | 0.25
av | 2 | 40 | 40 | 0.7984(61) | 0.7980(66) | 0.7954(53) | 0.7954(41) | 0.7983(41) | 0.20 | 0.20
at | 2 | 40 | 40 | 0.8083(44) | 0.8068(69) | 0.8062(59) | 0.8063(52) | 0.8095(48) | 0.20 | 0.20
s | 2 | 40 | 40 | 0.7635(52) | 0.7651(64) | 0.7660(51) | 0.7648(50) | 0.7674(30) | 0.20 | 0.20
Table 15:  Numerical results for the ground state mass in ensemble M5. $k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for $k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis, $\sigma_{G}$ is the energy smearing radius used for the Gaussian fits, $\sigma_{C}$ for the Cauchy fit. $C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{1}$ $k$-G |  $aE_{1}$ ($k+1$)-G |  $aE_{1}$ $k$-C |  $aE_{1}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 40 | 40 | 0.614(13) | 0.619(14) | 0.624(12) | 0.627(16) | 0.626(14) | 0.25 | 0.25
V | 3 | 40 | 40 | 0.621(16) | 0.611(17) | 0.618(14) | 0.623(19) | 0.619(14) | 0.30 | 0.30
T | 3 | 40 | 40 | 0.624(19) | 0.614(19) | 0.613(17) | 0.621(17) | 0.619(14) | 0.30 | 0.30
AV | 2 | 40 | 40 | 0.779(12) | 0.771(12) | 0.769(14) | 0.770(13) | 0.7733(80) | 0.25 | 0.20
AT | 2 | 40 | 40 | 0.779(16) | 0.775(10) | 0.771(13) | 0.773(13) | 0.7778(91) | 0.20 | 0.20
S | 2 | 40 | 40 | 0.740(12) | 0.737(10) | 0.734(14) | 0.739(15) | 0.7381(93) | 0.25 | 0.25
ps | 2 | 40 | 40 | 0.8773(93) | 0.8753(73) | 0.8793(99) | 0.8753(71) | 0.8788(76) | 0.20 | 0.20
v | 2 | 40 | 40 | 0.8801(99) | 0.8893(99) | 0.8852(99) | 0.8833(95) | 0.8899(85) | 0.20 | 0.25
t | 2 | 40 | 40 | 0.8841(98) | 0.8843(92) | 0.8812(82) | 0.8853(84) | 0.8899(85) | 0.20 | 0.25
av | 2 | 40 | 40 | 1.012(20) | 1.002(19) | 1.022(22) | 1.015(21) | 1.007(20) | 0.20 | 0.20
at | 2 | 40 | 40 | 1.042(16) | 1.040(14) | 1.033(13) | 1.029(12) | 1.033(14) | 0.20 | 0.20
s | 2 | 40 | 40 | 1.024(19) | 1.015(18) | 1.003(16) | 1.011(18) | 1.000(15) | 0.20 | 0.20
Table 16:  Numerical results for the first excited state mass in ensemble M5. $k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for $k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis, $\sigma_{G}$ is the energy smearing radius used for the Gaussian fits, $\sigma_{C}$ for the Cauchy fit. $C$ | $k$ | $N_{\rm source}$ | $N_{\rm sink}$ |  $aE_{2}$ $k$-G |  $aE_{2}$ ($k+1$)-G |  $aE_{2}$ $k$-C |  $aE_{2}$ ($k+1$)-C | $am_{C}$ | $\sigma_{G}/m_{C}$ | $\sigma_{C}/m_{C}$
---|---|---|---|---|---|---|---|---|---|---
PS | 2 | 40 | 40 | - | - | - | - | - | 0.25 | 0.25
V | 3 | 40 | 40 | 0.832(19) | 0.842(23) | 0.830(18) | 0.837(24) | 0.838(21) | 0.30 | 0.30
T | 3 | 40 | 40 | 0.831(23) | 0.839(26) | 0.840(23) | 0.832(20) | 0.838(21) | 0.30 | 0.30
AV | 2 | 40 | 40 | - | - | - | - | - | 0.25 | 0.20
AT | 2 | 40 | 40 | - | - | - | - | - | 0.20 | 0.20
S | 2 | 40 | 40 | - | - | - | - | - | 0.25 | 0.25
ps | 2 | 40 | 40 | - | - | - | - | - | 0.20 | 0.20
v | 2 | 40 | 40 | - | - | - | - | 1.083(50) | 0.20 | 0.25
t | 2 | 40 | 40 | - | - | - | - | 1.083(50) | 0.20 | 0.25
av | 2 | 40 | 40 | - | - | - | - | - | 0.20 | 0.20
at | 2 | 40 | 40 | - | - | - | - | - | 0.20 | 0.20
s | 2 | 40 | 40 | - | - | - | - | - | 0.20 | 0.20
Table 17:  Numerical results for the second excited state mass in ensemble M5.
$k$-G stands for $k$-Gauss fit, ($k+1$)-G is $k+1$-Gauss fit, $k$-C stands for
$k$-Cauchy function fit, ($k+1$)-G $k+1$-Cauchy function fit, $C$ indicates
the mesonic channel considered, $am_{C}$ is the result of the GEVP analysis,
$\sigma_{G}$ is the energy smearing radius used for the Gaussian fits,
$\sigma_{C}$ for the Cauchy fit.

##  VII Summary and outlook

The first purpose of this paper is to report on software development and
testing related to a new analysis package that implements the HLT method as a
spectral density reconstruction tool for two-point functions obtained on the
lattice. We performed a systematic study of the method itself, in order to
optimize the parameter choices entering the spectral reconstruction. We
provide details, and numerical examples, in the main body of the paper. Our
main findings can be summarized as follows.

  * •

For a given smearing kernel, $\Delta_{\sigma}(E,\,\omega)$, and a given input
set of 2-point functions, $C_{ij}(t)$, spectral density reconstruction
requires minimizing a functional, $W[\vec{g}]$, defined in Eq. (48). This
procedure depends on two unphysical parameters, $\lambda$ and $\alpha$. We
identify a range of $\lambda\sim{\cal O}(10)$ for which the smeared spectral
densities do not depend on the parameters $\lambda$ and $\alpha$ within the
statistical error, which is nonetheless not substantial.

  * •

In this work, we were able to afford values of the smearing radius, $\sigma$,
in the range $0.18\,m_{0}\leq\sigma\leq 0.35\,m_{0}$, where $m_{0}$ is the
mass of the ground state appearing in the two-point correlation functions used
as input data.

  * •

Under the conditions identified at the previous points, our fit results are
independent of the specific choice of smearing kernel. We illustrate this
point by repeating our analysis with both Gaussian and Cauchy kernel, and
demonstrating that the results of the two processes are compatible with one
another.

  * •

APE and Wuppertal smearings are essential in the production of correlation
functions to be analyzed. APE smearing is necessary in order to explore the
high-energy behavior of the spectral density. Wuppertal smearing must be tuned
so that contributions to the spectral density of all states of interest have
comparable amplitudes.

  * •

A critical quantity for spectral density reconstruction, in particular in
reference to the identification of excited states, is the time extent of the
lattice, $N_{t}$. We illustrate its impact on the results of the physical
analysis, by comparing lattices with $N_{t}=48$, $N_{t}=64$, and $N_{t}=96$,
while keeping all other parameters fixed. For the longest time extent, it is
possible to reconstruct reliably the ground state, first and second excited
states.

The second purpose of this document is to report on progress in the study of
meson spectroscopy in the $Sp(4)$ lattice gauge theory coupled to $N_{f}=2$
Dirac fermions transforming according to the fundamental representation of the
gauge group, and $n_{f}=3$ on the two-index antisymmetric one. This theory is
a prominent candidate for new physics in the context of composite Higgs models
implementing top partial compositeness Barnard:2013zea . To this purpose, we
generated, using the (R)HMC algorithm, five new gauge ensembles. Allowing for
thermalization, we selected configurations to retain in the ensembles so that
the plaquette does not show significant indications of autocorrelation. We
applied the Wilson flow as a smoothening process and measured the topological
charge of the configurations, to demonstrate the absence of topological
freezing. Some modest level of autocorrelation is visible in the topology, but
we do not expect this phenomenon to affect meson spectroscopy studied in this
work.

We then measured the two-point correlation functions involving twelve flavored
meson operators. We applied both the HLT spectral density method, and the GEVP
analysis, and reported the resulting mass spectra in Tabs. 4 to 17. We
expressed the masses in units of the lattice spacing, $a$, and verified that
the results were independent of the analysis process applied, within
statistical uncertainties. Besides the fact that the analysis involves new
ensembles obtained for different choices of lattice parameters, the main
advance with respect to the literature in this theory is that we measured
excited-state masses, including both first and second excited states. Figure
20 is a way to visualize the results: we display the masses, $\hat{m}$,
expressed in units of the Wilson flow scale, $w_{0}$, of all the states we
could identify in the spectral density analysis, and for all the available
channels and ensembles.

The work reported in this paper sets the stage for a plethora of future
studies, both within the ongoing program of systematic exploration of $Sp(2N)$
gauge theories, but also in contexts of more general interest. Firstly, the
new ensembles can be analyzed in other physical channels, looking at
correlation functions involving the singlet meson operators and the chimera
baryon operators. Such measurements are currently underway. Second, having put
in place all the necessary technology, it is possible to generate and analyze
new ensembles with different couplings, $\beta$, and masses, $am^{\mathrm{f}}$
and $am^{\rm as}$, of the fermions. By combining the results within the
framework of Wilson chiral perturbation theory Sheikholeslami:1985ij ;
Rupak:2002sm  (see also Refs Sharpe:1998xm , and, in the context of
improvement: Refs Symanzik:1983dc ; Luscher:1996sc , as well as, for baryon
chiral perturbation theory in QCD: Refs. Jenkins:1990jv ; Bernard:1995dp ;
Beane:2003xv ), one can perform preliminary extrapolations towards the
continuum and chiral limits. Doing so will allow us to classify the properties
of the theory that are of relevance to model building considerations in the
CHM context, and to compare the spectroscopy results with those descending
from other, complementary approaches Bizot:2016zyu ; Erdmenger:2020lvq ;
Erdmenger:2020flu ; Elander:2020nyd ; Elander:2021bmt ; Erdmenger:2023hkl ;
Erdmenger:2024dxf .

The third further direction of future development reverts back to the original
motivation of spectral density studies: the reconstruction of off-shell
properties of correlation functions. One wants to put in place the technology
needed to measure not just the position of the poles in the two-point
functions, but also their residues, which are related to the decay constants
and the overlap functions between operators and individual states. Generally,
one would like to reconstruct the whole complex-space structure of the
correlators, because these enter observable quantities of theoretical and
phenomenological relevance, as discussed in the introduction.

Last but not least, we wrote the main body of this paper by providing complete
technical details, because our systematic study has general implications for
any lattice gauge theory, including also QCD, the theory of strong nuclear
interactions. We expect the progress we reported here to be useful in a
broader context, as spectral densities have broad ranges of phenomenological
applications. We made all our new software and numerical results available as
open source, so that users can download them, use them, and modify them as
appropriate to their own physics goals. Given the potential that spectral
density analysis has for lattice field theory, we envision this work to
contribute towards the general goal of the lattice community to develop and
master such novel techniques, and turn them into high-precision analysis
tools.

###### Acknowledgements.

We would like to thank Giacomo Cacciapaglia, Thomas Flacke, Anna Hasenfratz,
Chulwoo Jung, Gabriele Ferretti, and Sarada Rajeev, for useful discussions
during the “PNU Workshop on Composite Higgs: Lattice study and all”, at
Haeundae, Busan, in February 2024, where preliminary results of this study
were presented. The work of E.B. and B.L. is supported in part by the EPSRC
ExCALIBUR programme ExaTEPP (project EP/X017168/1). The work of E.B. has also
been supported by the UKRI Science and Technology Facilities Council (STFC)
Research Software Engineering Fellowship EP/V052489/1. The work of E.B., B.L.,
and M.P. has been supported in part by the STFC Consolidated Grant No.
ST/X000648/1. The work of N.F. has been supported by the STFC Consolidated
Grant No. ST/X508834/1. A.L. is funded in part by l’Agence Nationale de la
Recherche (ANR), under grant ANR-22-CE31-0011. The work of D.K.H. was
supported by Basic Science Research Program through the National Research
Foundation of Korea (NRF) funded by the Ministry of Education
(NRF-2017R1D1A1B06033701). The work of D.K.H. was further supported by the
National Research Foundation of Korea (NRF) grant funded by the Korea
government (MSIT) (2021R1A4A5031460). L.D.D. and R.C.H. are supported by the
UK Science and Technology Facility Council (STFC) grant ST/P000630/1. The work
of L.D.D. was supported by the ExaTEPP project EP/X01696X/1. The work of
J.W.L. was supported by IBS under the project code, IBS-R018-D1. The work of
H.H. and C.J.D.L. is supported by the Taiwanese MoST grant
109-2112-M-009-006-MY3 and NSTC grant 112-2112-M-A49-021-MY3. The work of
C.J.D.L. is also supported by Grants No. 112-2639-M-002-006-ASP and No.
113-2119-M-007-013-. The work of B.L. and M.P. has been supported in part by
the STFC Consolidated Grant No. ST/T000813/1. B.L., M.P. and L.D.D. received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program under Grant Agreement No. 813942.
The work of D.V. is supported by STFC under Consolidated Grant No.
ST/X000680/1. The work of F.Z. is supported by the STFC Grant No.
ST/X000648/1. Numerical simulations have been performed on the DiRAC Extreme
Scaling service at the University of Edinburgh, and on the DiRAC Data
Intensive service at Leicester. The DiRAC Extreme Scaling service is operated
by the Edinburgh Parallel Computing Centre on behalf of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). This equipment was funded by BEIS capital funding
via STFC capital grant ST/R00238X/1 and STFC DiRAC Operations grant
ST/R001006/1. DiRAC is part of the National e-Infrastructure. Spectral density
measurements have been performed by using the LSDensities software package in
Ref. Forzano:2024 . Research Data Access Statement—The analysis code and data
generated for this manuscript can be downloaded from Refs. analysis_release
and data_release , respectively. Open Access Statement—For the purpose of open
access, the authors have applied a Creative Commons Attribution (CC BY)
license to any Author Accepted Manuscript version arising.

##  Appendix A More about spectral density reconstruction

In this appendix, we show details of our implementation for the HLT method. As
described in Sec. IV.1, the functional $W[\vec{g}]$ in Eq. (48) consists of
two parts: the systematic-error functional $A[\vec{g}]/A_{0}$ and the
statistical-error functional $B/B_{\rm norm}$. Minimisation of $W[\vec{g}]$ is
tantamount to solving the linear system

$\vec{g}=\left(\mathtt{S}+\frac{\lambda\,A_{0}}{B_{\rm
norm}}\;\mathtt{B}\right)^{-1}\vec{\mathtt{f}}\,.$ (59)

where $B_{\rm norm}$ is defined below Eq. (48). The matrix $\mathtt{S}$
descends from minimising in $A[\vec{g}]/A_{0}$ the contribution quadratic in
$\vec{g}$. Its matrix elements can be always expressed as follows

$\begin{split}\mathtt{S}_{tr}=\int_{E_{\rm min}}^{\infty}&dEe^{\alpha
E}\,b(t+1,E)\,b(r+1,E)=\\\ &=\frac{e^{E_{\rm min}(\alpha-
r-t-2)}}{t+r+2-\alpha}+\frac{e^{E_{\rm