# Inclusive rates from smeared spectral densities
in the two-dimensional O(3) non-linear $\sigma$-model

John Bulava<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron DESY,
Platanenallee 6, 15738 Zeuthen, Germany  Maxwell T. Hansen
<EMAIL_ADDRESS>Higgs Centre for Theoretical Physics, School of
Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK
Michael W. Hansen<EMAIL_ADDRESS>University of Graz, Institute
for Physics, A-8010 Graz, Austria  Agostino Patella
<EMAIL_ADDRESS>Institut für Physik und IRIS Adlershof,
Humboldt-Universität zu Berlin, Zum Großen Windkanal 6, D-12489 Berlin,
Germany  Nazario Tantalo<EMAIL_ADDRESS>University and INFN of
Roma Tor Vergata, Via della Ricerca Scientifica 1, I-00133, Rome, Italy

###### Abstract

This work employs the spectral reconstruction approach of Ref. Hansen:2019idp
to determine an inclusive rate in the $1+1$ dimensional O(3) non-linear
$\sigma$-model, analogous to the QCD part of
${e}^{+}{e}^{-}\rightarrow\rm{hadrons}$. The Euclidean two-point correlation
function of the conserved current $j$ is computed using Monte Carlo lattice
field theory simulations for a variety of spacetime volumes and lattice
spacings. The spectral density of this correlator is related to the inclusive
rate for $j\rightarrow{\rm X}$ in which all final states produced by the
external current are summed. The ill-posed inverse problem of determining the
spectral density from the correlation function is made tractable through the
determination of smeared spectral densities in which the desired density is
convolved with a set of known smearing kernels of finite width $\epsilon$. The
smooth energy dependence of the underlying spectral density enables a
controlled $\epsilon\to 0$ extrapolation in the inelastic region, yielding the
real-time inclusive rate without reference to individual finite-volume
energies or matrix elements. Systematic uncertainties due cutoff effects and
residual finite-volume effects are estimated and taken into account in the
final error budget. After taking the continuum limit, the results are
consistent with the known analytic rate to within the combined statistical and
systematic errors. Above energies where 20-particle states contribute, the
overall precision is sufficient to discern the four-particle contribution to
the spectral density.

††preprint:  DESY 21-201, HU-EP-21/49

###### Contents

  1. I Introduction
  2. II General framework in continuous infinite volume
  3. III Finite-volume estimator
  4. IV Lattice-discretized model
  5. V Numerical Results
    1. V.1 Fixed smearing width
    2. V.2 Extrapolation to zero smearing width
  6. VI Conclusions
  7. A Analytic expressions for the spectral density
  8. B Spectral reconstruction implementation
  9. C Finite-volume effects in a model smeared spectral density
    1. C.1 $1/L$ cancellation
    2. C.2 Exponentially suppressed volume effects
  10. D Simulation algorithm

##  I Introduction

Markov chain Monte Carlo simulations of lattice QCD continue to provide a
quantitative window into the strong nuclear force. However, the Euclidean
metric signature required for Monte Carlo sampling of the QCD path integral
complicates the study of real-time scattering processes. While the spatial
volume dependence of energies and matrix elements has successfully been used
as a probe of few-particle scattering and decay amplitudes Luscher:1990ux ;
Lellouch:2000pv , this approach is inapplicable to energies above arbitrary
multi-particle thresholds111For recent reviews on the status of treating
three-particle states see Refs. Hansen:2019nir ; Rusetsky:2019gyk . and is
restricted to the discrete set of center-of-mass energies that appear in a
given finite-volume geometry. This last disadvantage hampers the continuum
limit of amplitudes at fixed energy, which can only be achieved by maintaining
a constant physical volume as the lattice spacing is decreased. Using the
finite-volume formalism to determine inclusive rates is also difficult as
these can only be calculated by extracting all contributing amplitudes
separately and summing the individual rates.

An alternative approach, following the seminal work of Ref. Barata:1990rn , is
to use Euclidean correlation functions computed in lattice simulations to
determine spectral densities which are definitionally independent of the
metric signature. This alternative is not subject to the same limitations as
finite-volume methods and has been formally developed for inclusive rates
mediated by an external current Liu:2016djw ; Hansen:2017mnd ; Gambino:2020crt
as well as arbitrary scattering and transition amplitudes Bulava:2019kbi ;
Bruno:2020kyl .222This approach is inspired in part by the extraction of
analogous observables in lattice calculations at finite temperature. See Ref.
Meyer:2011gj  for an instructive review. A major obstacle of this program is
the solution of an ill-posed inverse problem to determine continuous spectral
densities from correlation functions sampled on a finite set of time
separations with statistical errors. In Refs. BG1 ; BG2  Backus and Gilbert
deal with this by instead computing spectral densities smeared with a kernel
that is known _a posteriori_. An important advancement by Ref. Hansen:2019idp
(see also Refs. pt_hlt ; Bailas:2020qmv ) enables the smearing kernel to be
specified _a priori_. A review and comparison of alternative spectral
reconstruction techniques goes beyond the scope of this work Maass1990_0e27e ;
Tripolt:2018xeo .

Although these methods are not new, their application to compute inclusive
rates and scattering amplitudes from actual Monte Carlo simulation data, with
a detailed analysis of both statistical and systematic errors, has not yet
been performed. It is therefore worthwhile to pursue such a test in a
controlled setting. This work employs the spectral reconstruction strategy of
Ref. Hansen:2019idp  to determine an inclusive rate in the two-dimensional
O(3) non-linear $\sigma$-model. As in Ref. Luscher:1990ck , comparison is made
with exact continuum results. However, our work overcomes the limitations
mentioned above: the inclusive rate is computed above inelastic thresholds and
the continuum limit taken at fixed center-of-mass energy without the need for
a constant physical volume. Large physical simulation volumes are however
required to control finite-volume effects in the smeared spectral density at
fixed smearing width. After taking the continuum limit and accounting for
systematic errors due to finite-volume effects, the desired spectral density
is obtained by extrapolating the smearing width to zero.

For this first application we treat a process akin to
${e}^{+}{e}^{-}\rightarrow{\rm hadrons}$, the QCD component of which is given
by the spectral density

$\displaystyle\rho_{\mu\nu}(k)=\frac{1}{2\pi}\int d^{4}x\,e^{-ik\cdot
x}\langle\Omega|\hat{j}^{\rm em}_{\mu}(x)\,\hat{j}^{\rm
em}_{\nu}(0)|\Omega\rangle=(g_{\mu\nu}k^{2}-k_{\mu}k_{\nu})\,\rho(k^{2}),$ (1)

where the integral is performed in Minkowski space, $\hat{j}^{\rm em}_{\mu}$
is the quark-level electromagnetic current and the relation to the physical
process is

$\displaystyle\rho(s)=\frac{R(s)}{12\pi^{2}},\quad
R(s)=\frac{\sigma\left[{e}^{+}{e}^{-}\rightarrow{\rm
hadrons}\right](s)}{4\pi\alpha_{\rm em}(s)^{2}/(3s)}.$ (2)

This celebrated ‘$R$-ratio’ has a number of phenomenological applications
including the hadronic vacuum polarization contribution to the anomalous
magnetic moment of the muon. It can be obtained via spectral reconstruction of
the Euclidean correlator in the time-momentum representation

$\displaystyle C(t)=\int d^{3}\bm{x}\,\langle\Omega|\hat{j}^{\rm
em}_{z}(\bm{x})\,e^{-\hat{H}t}\,\hat{j}^{\rm
em}_{z}(0)^{\dagger}|\Omega\rangle=\int_{0}^{\infty}d\omega\,\omega^{2}\rho(\omega^{2})\,e^{-\omega
t}\,,$ (3)

computed in lattice QCD simulations.333As written, this expression holds for
both the anti-hermitian definition of $\hat{j}^{\rm em}_{z}(0)$,
conventionally used with the Euclidean metric, as well as the hermitian
Minkowski current implicitly used in Eq. (1). Everything discussed in this
work is readily transferable to the determination of $\rho(s)$ from $C(t)$ in
lattice QCD. While computations of $\rho(s)$ are naturally compared with
experimental determinations of $R(s)$, numerical results for the spectral
density computed here are compared with the corresponding analytical
predictions in Fig. 9, which is our main result.

This remainder of this work is organized as follows. Sec. II defines the
spectral density of interest in the O(3) model, outlines the reconstruction
approach, and suggests a strategy for extrapolating the smearing width to
zero. Sec. III examines the influence of the finite torus on which the
simulations are performed while Sec. IV defines the lattice regularization and
discusses the continuum limit. Sec. V presents numerical results and Sec. VI
concludes.

##  II General framework in continuous infinite volume

This section introduces the two-dimensional O(3) non-linear $\sigma$-model and
the Euclidean correlator $C(t)$ with spectral density $\rho(\omega)$, which
are both defined via the conserved vector current. We additionally review the
algorithm of Ref. Hansen:2019idp  for systematically determining a smeared
version of $\rho(\omega)$ from numerical estimates of $C(t)$ and detail an
extrapolation procedure for taking the smearing width to zero. In this section
the Euclidean spacetime is assumed to be continuous and infinite in both
directions. Peculiarities due to the finite simulation volume are discussed in
Sec. III and the lattice discretization in Sec. IV.

The continuum Euclidean action of the two-dimensional O(3) model is defined as

$\displaystyle S[\sigma]=\frac{1}{2g^{2}}\int
d^{2}x\,\partial_{\mu}\sigma(x)\cdot\partial_{\mu}\sigma(x)\,,$ (4)

where the $3$-component real field $\sigma^{a}(x)$ has unit length
$\sigma(x)\cdot\sigma(x)=1$. The O(3) model is asymptotically free, has a
dynamically-generated mass gap $m$, and is integrable. It also possesses a
global O(3) symmetry which rotates the field $\sigma(x)$. The corresponding
Noether current is given by

$\displaystyle
j_{\mu}^{c}(x)=\frac{1}{g^{2}}\epsilon^{abc}\sigma^{a}(x)\partial_{\mu}\sigma^{b}(x)\,,$
(5)

where $\epsilon^{abc}$ is the Levi-Civita tensor and repeated indices are
summed. This current transforms irreducibly under the $I=1$ (fundamental)
representation of O(3).

The aim of this work is to reconstruct the spectral density $\rho(\omega)$
associated with $j_{\mu}^{c}$ using Euclidean correlation functions determined
numerically from lattice Monte Carlo simulations. The spectral density is
defined implicitly via the relation

$\displaystyle
2\pi\,\langle\Omega|\hat{j}^{a}_{\mu}(0)\,\delta^{2}(\hat{P}-p)\,\hat{j}^{b}_{\nu}(0)|\Omega\rangle=\frac{\delta^{ab}}{3}\left(\delta_{\mu\nu}-\frac{p_{\mu}p_{\nu}}{p^{2}}\right)\rho\left(\sqrt{p^{2}}\right)\,,$
(6)

where $\hat{P}=(\hat{H},\hat{\bm{P}})$ is the 2-momentum operator and $p$ is
an external 2-momentum coordinate with components denoted $p=(E,\bm{p})$.
Covariance under both the O(3) internal symmetry and the Euclidean spacetime
SO(2) symmetry, as well as current conservation, enable the factorization in
Eq. (6). By specializing $\mu=\nu=1$ and $p=(E,0)$, and contracting the
internal indices, one obtains the more direct definition

$\displaystyle\rho(E)=2\pi\,\langle\Omega|\hat{j}^{a}_{1}(0)\,\delta^{2}(\hat{P}-p)\,\hat{j}^{a}_{1}(0)|\Omega\rangle\,.$
(7)

The utility of the two-dimensional O(3) model for the present study is that
(due to integrability) $\rho(E)$ can be computed analytically, enabling a
comparison between our numerical reconstruction and the exact result. The
spectral density decomposes into sectors defined by the number $n$ of
asymptotic particles444The $n$-particle states interpolated by the current
also transform irreducibly under the $I=1$ representation of O(3). propagating
between the two currents

$\rho(E)=\sum_{\text{even }n\geq 2}\rho^{(n)}(E)\,.$ (8)

The contribution $\rho^{(n)}(E)$ has support for $E>n\,m$, where $m$ is the
mass gap. Even though integrability (together with a number of mild and
generally accepted assumptions) fixes $\rho^{(n)}(E)$ for every even value of
$n$, explicit expressions have been worked out only for $n=2,4,6$. At the
energies considered here the sum over the number of particles is rapidly
convergent and the $n=6$ contribution is at least a couple of orders of
magnitude smaller than the $n=4$ contribution. In the following we therefore
refer to the spectral density summed over $n=2,4,6$ as the exact spectral
density. The $n=2$ contribution can be written in closed form

$\displaystyle\rho^{(2)}(E)=\frac{3\pi^{3}}{8\theta^{2}}\,\frac{\theta^{2}+\pi^{2}}{\theta^{2}+4\pi^{2}}\,\tanh^{3}\frac{\theta}{2}\,\,\bigg{|}_{\theta=2\cosh^{-1}\frac{E}{2m}}\,,$
(9)

while the $n=4,6$ contributions are expressed in terms of phase-space
integrals that require numerical evaluation as discussed in App. A.

This spectral density is related to the Euclidean-signature current-current
correlator at zero spatial momentum via the Laplace transform

$\displaystyle C(t)\equiv\int
d\bm{x}\,\langle\Omega|\,\hat{j}_{1}^{a}(0,\bm{x})\,e^{-\hat{H}t}\,\hat{j}_{1}^{a}(0)\,|\Omega\rangle=\int_{0}^{\infty}d\omega\,e^{-\omega
t}\,\rho(\omega)\,.$ (10)

The demonstration of a systematic method to invert this relation, given
realistic numerical estimates of $C(t)$ at a finite set of time slices with
statistical errors, is the central focus of this work.

A ubiquitous problem, faced across a breadth of scientific disciplines, is the
numerically ill-conditioned nature of the inverse Laplace transform.555See for
example Ref. Maass1990_0e27e  and references therein. A promising way forward
is to recognize that $\rho(\omega)$ is not directly extractable from numerical
data and that one should instead target the smeared spectral density

$\displaystyle\rho_{\epsilon}(E)=\int_{0}^{\infty}d\omega\,\delta_{\epsilon}(E,\omega)\,\rho(\omega)\,,$
(11)

where $\delta_{\epsilon}(E,\omega)$ is any approximation of the Dirac
$\delta$-function satisfying $\lim_{\epsilon\to
0}\delta_{\epsilon}(E,\omega)=\delta(E-\omega)$ and
$\int_{-\infty}^{\infty}d\omega\,\delta_{\epsilon}(E,\omega)=1$. The challenge
of recovering $\rho_{\epsilon}(E)$ from $C(t)$ can then be made arbitrarily
mild (or severe) by varying the specific functional form of
$\delta_{\epsilon}(E,\omega)$ and the values of $\epsilon$ and $E$.

Following Ref. Hansen:2019idp  we consider smearing kernels
$\delta_{\epsilon}(E,\omega)$ that can be represented exactly as

$\displaystyle{\delta}_{\epsilon}(E,\omega)=a\sum_{\tau=1}^{\infty}g_{\tau}^{\rm
target}\,b_{\tau}(\omega)\,,$ (12)

where $\tau$ is a dimensionless integer variable, $a$ an arbitrary scale with
dimensions of inverse energy to be later identified with the lattice spacing,
the $b_{\tau}(\omega)$ are basis functions, and $g_{\tau}^{\rm target}\equiv
g_{\tau}^{\rm target}(\epsilon,E)$ components of the target smearing kernel
${\delta}_{\epsilon}(E,\omega)$ in this basis. If we identify
$b_{\tau}(\omega)={\rm e}^{-a\omega\tau}$, then the smeared spectral density
is given by

${\rho}_{\epsilon}(E)=a\sum_{\tau=1}^{\infty}g^{\rm
target}_{\tau}C(a\tau)=\int_{0}^{\infty}d\omega\,\delta_{\epsilon}(E,\omega)\,\rho(\omega)\,.$
(13)

The choice $b_{\tau}(\omega)={\rm e}^{-a\omega\tau}$ is based on the non-
essential assumption that spacetime is infinite. In practice, as explained in
Sec. III, we use a basis that takes into account the periodicity of the finite
temporal direction.

An estimator for $\rho_{\epsilon}(E)$ is obtained by approximating
$\delta_{\epsilon}(E,\omega)$ with an element of the space spanned by a finite
number of basis functions. The coefficients $g_{\tau}$ corresponding to the
approximation of $\delta_{\epsilon}(E,\omega)$ are determined by minimizing
the functional666Possible additional constraints on the minimization are
discussed in App. B.

$\displaystyle W_{\lambda}[g]=(1-\lambda)\frac{A[g]}{A[0]}+\lambda B[g]\,,$
(14)

where $\lambda\in[0,1]$ is the ‘trade-off’ parameter (discussed shortly) and
the functionals $A[g]$ and $B[g]$ are given by777The relative normalization of
the $A[g]$ and $B[g]$ functionals in Eq. (14) differs from Ref. Hansen:2019idp
. The factor $1/(aC(0))^{2}$ in the normalization of $B[g]$ present there is
not included here since $aC(0)\simeq 1$ while the factor $1/A[0]$ in the first
term of Eq. (14) makes it dimensionless ($A[0]\propto 1/\epsilon$).

$\displaystyle
A[g]=\int_{E_{0}}^{\infty}d\omega\left\\{{\delta}_{\epsilon}(E,\omega)-a\sum_{\tau=1}^{\tau_{\rm
max}}g_{\tau}\,b_{\tau}(\omega)\right\\}^{2}\,,\qquad
B[g]=\sum_{\tau,\tau^{\prime}=1}^{\tau_{\rm
max}}g_{\tau}g_{\tau^{\prime}}\,{\rm
Cov}\left[aC(a\tau),aC(a\tau^{\prime})\right]\,.$ (15)

As discussed in App. B, $A[g]$ and $B[g]$ are both positive quadratic forms so
that the minimum conditions

$\displaystyle\left.\frac{\partial W_{\lambda}[g]}{\partial
g_{\tau}}\right|_{g_{\tau}=g_{\tau}^{\lambda}}=0\,,$ (16)

are a linear system of equations to be solved for the coefficients
$g_{\tau}^{\lambda}\equiv g_{\tau}^{\lambda}(\epsilon,E,\tau_{\rm max})$.
These coefficients define the approximation of $\delta_{\epsilon}(E,\omega)$
and the associated estimator for ${\rho}_{\epsilon}(E)$ according to

$\displaystyle{\delta}^{\lambda}_{\epsilon}(E,\omega)=a\sum_{\tau=1}^{\tau_{\rm
max}}g^{\lambda}_{\tau}\,b_{\tau}(\omega)\;,\qquad{\rho}^{\lambda}_{\epsilon}(E)=a\sum_{\tau=1}^{\tau_{\rm
max}}g^{\lambda}_{\tau}C(a\tau)=\int_{0}^{\infty}d\omega\,\delta^{\lambda}_{\epsilon}(E,\omega)\,\rho(\omega)\;.$
(17)

The functional $B[g^{\lambda}]$ is simply the statistical variance of
$\rho^{\lambda}_{\epsilon}(E)$ and therefore vanishes in the ideal case of
infinitely precise input data. On the other hand, $A[g^{\lambda}]$ measures
the distance between the target kernel ${\delta}_{\epsilon}(E,\omega)$ and its
approximation $\delta^{\lambda}_{\epsilon}(E,\omega)$ in the range888The
parameter $E_{0}$ can be adjusted by exploiting the fact that $\rho(E)$ has
support only for $E>2m$, so that $\rho_{\epsilon}(E)$ in Eq. (11) is
insensitive to $\delta_{\epsilon}(E,\omega)$ for $\omega<2m$. The same holds
for $\rho^{\lambda}_{\epsilon}(E)$ so that the functional form of
$\delta^{\lambda}_{\epsilon}(E,\omega)$ can be left unconstrained for
$\omega<2m$. Any $E_{0}\leq 2m$ is therefore a viable choice in determining
the coefficients $g_{\tau}^{\lambda}$ so $E_{0}$ can be chosen to improve the
numerical stability of the minimization procedure. $E\in[E_{0},\infty]$. It
can only vanish in the limit $\tau_{\rm max}\rightarrow\infty$ when
$g^{\lambda}_{\tau}\rightarrow g^{\rm target}_{\tau}$. The coefficients
$g^{\lambda}_{\tau}$ that minimize $W_{\lambda}[g]$ thus represent a
particular balance between statistical and systematic errors, as dictated by
the $\lambda$ parameter. For small $\lambda$ the estimator
$\rho^{\lambda}_{\epsilon}(E)$ is close to $\rho_{\epsilon}(E)$ but with a
large statistical uncertainty. Conversely, for large $\lambda$ the estimator
$\rho^{\lambda}_{\epsilon}(E)$ has a small statistical error but differs
significantly from $\rho_{\epsilon}(E)$.

When evaluated at the minimum, the functional $W_{\lambda}[g]$ is a function
of $\lambda$ only, thus defining $W(\lambda)\equiv W_{\lambda}[g^{\lambda}]$.
The recipe suggested in Ref. Hansen:2019idp  to choose the optimal value of
the trade-off parameter defines $\lambda_{\star}$ such that

$\displaystyle\left.\frac{\partial
W(\lambda)}{\partial\lambda}\right|_{\lambda=\lambda_{\star}}=0\;.$ (18)

A straightforward application of Eq. (16) demonstrates that at
$\lambda_{\star}$ (the maximum of $W(\lambda)$ where
$g_{\tau}^{\lambda}=g_{\tau}^{\star}$) one has
$A[g^{\star}]=A[0]B[g^{\star}]$. This can be understood as the condition of
‘optimal balance’ between statistical and systematic errors. The numerical
results presented here are obtained using this recipe. In the remainder of
this work, unless explicitly stated, we do not distinguish the theoretical
quantity $\rho_{\epsilon}(E)$ from its numerical estimator
$\rho^{\lambda_{\star}}_{\epsilon}(E)$. Estimates of the systematic error on
$\rho_{\epsilon}(E)$ induced by the residual difference
$\delta^{\lambda_{\star}}_{\epsilon}(E,\omega)-\delta_{\epsilon}(E,\omega)$
are discussed in Sec. V.

A defining feature of the approach in Ref. Hansen:2019idp  is that the
smearing kernel $\delta_{\epsilon}(E,\omega)$ and the associated values of
$\epsilon$ and $E$ are inputs of the algorithm, in contrast to the original
Backus-Gilbert method BG1 ; BG2 . This work exploits this by employing four
functional forms for the smearing kernel, each of which is a function of
$x=E-\omega$:

$\displaystyle\delta_{\epsilon}^{\sf g}(x)$
$\displaystyle=\frac{1}{\sqrt{2\pi}\epsilon}\exp\left[-\frac{x^{2}}{2\epsilon^{2}}\right],$
$\displaystyle\delta_{\epsilon}^{{\sf c}0}(x)$
$\displaystyle=\frac{1}{\pi}\frac{\epsilon}{x^{2}+\epsilon^{2}}\,,$ (19)
$\displaystyle\delta_{\epsilon}^{{\sf c}1}(x)$
$\displaystyle=\frac{2}{\pi}\frac{\epsilon^{3}}{(x^{2}+\epsilon^{2})^{2}}\,,$
$\displaystyle\delta_{\epsilon}^{{\sf c}2}(x)$
$\displaystyle=\frac{8}{3\pi}\frac{\epsilon^{5}}{(x^{2}+\epsilon^{2})^{3}}\,.$
(20)

Here ${\sf g}$ and ${\sf c}$ stand for ‘Gauss’ and ‘Cauchy’ respectively, and
the number following ${\sf c}$ gives the order of the Cauchy-like pole as
shown. All kernels are normalized to unit area.

x |  $w^{{\sf x}}_{k},$ even $k$ |  $w^{{\sf x}}_{k},$ odd $k$ |  $\phantom{+}w^{{\sf x}}_{1}$ |  $\phantom{+}w^{{\sf x}}_{2}$ |  $\phantom{+}w^{{\sf x}}_{3}$ |  $\phantom{+}w^{{\sf x}}_{4}$
---|---|---|---|---|---|---
g | $\displaystyle\frac{k!}{(-2)^{k/2}(k/2)!}$ | $0$ | $\phantom{+}0$ | $-1$ | $\phantom{+}0$ | $\phantom{+}3$
c0 | $\displaystyle\phantom{\frac{1}{2}}1\phantom{\frac{1}{2}}$ | $\displaystyle\phantom{\frac{1}{2}}1\phantom{\frac{1}{2}}$ | $\phantom{+}1$ | $\phantom{+}1$ | $\phantom{+}1$ | $\phantom{+}1$
c1 | $\phantom{\frac{1}{2}}(1-k)\phantom{\frac{1}{2}}$ |  $(1-k)$ | $\phantom{+}0$ | $-1$ | $-2$ | $-3$
c2 |  $\displaystyle\frac{1}{3}(k-3)(k-1)$ |  $\displaystyle\frac{1}{3}(k-3)(k-1)$ | $\phantom{+}0$ | $-1/3$ | $\phantom{+}0$ | $\phantom{+}1$
Table 1: Summary of the kernel-dependent coefficients $w^{{\sf x}}_{k}$
entering the small-$\epsilon$ expansion in Eq. (21). Although not reflected in
the general formulae, $w_{3}^{{\sf c1}}$ and $w_{5}^{{\sf c2}}$ are the lowest
non-zero coefficients with odd $k$ for their respective smearing kernels.

Given estimates of $\rho_{\epsilon}(E)$ over a range of $\epsilon$ for each of
the resolution functions shown above, the final step in determining $\rho(E)$
is to perform an $\epsilon\to 0$ extrapolation. To this end it is useful to
understand the small-$\epsilon$ expansion at fixed $E$. We consider this
expansion only at energies away from singularities, i.e. away from
$2m\mathbb{Z}^{+}$. At such points the smeared spectral density satisfies

$\displaystyle\rho^{{\sf{x}}}_{\epsilon}(E)$
$\displaystyle\equiv\int_{0}^{\infty}d\omega\,\delta_{\epsilon}^{\sf
x}(E-\omega)\,\rho(\omega)=\rho(E)+\sum_{k=1}^{\infty}w^{{\sf
x}}_{k}a_{k}(E)\epsilon^{k}\,,$ (21)

where the superscript $\sf x$ labels a particular smearing kernel. As
indicated in the rightmost expression, the $O(\epsilon^{k})$ contribution to
the expansion factorizes into a geometric kernel-dependent coefficient
$w^{{\sf x}}_{k}$ and a kernel-independent piece $a_{k}(E)$ which depends on
$\rho(E)$ only. Tab. 1 summarizes the values of $w^{{\sf x}}_{k}$ for the four
kernels used in this work. The ambiguity in separating $w^{{\sf x}}_{k}$ and
$a_{k}(E)$ is fixed by setting $w^{{\sf c0}}_{k}=1$ for all $k$. For this
choice, all remaining $w^{{\sf x}}_{k}$ are rational numbers and

$a_{k}(E)=\begin{cases}\frac{(-1)^{k/2}}{k!}\left(\frac{d}{dE}\right)^{k}\rho(E)\,,\qquad&\textrm{$k$
even}\\\ \\\ \lim_{\eta\to
0^{+}}\frac{(-1)^{(k-1)/2}}{2\pi}\int_{-\infty}^{\infty}{\rm
d}\omega\,\frac{\rho(E+\omega)+\rho(E-\omega)}{(\omega+i\eta)^{k+1}}\,,\qquad&\textrm{$k$
odd}\end{cases}\ .$ (22)

##  III Finite-volume estimator

The theory is considered here on an $L\times T$ volume with periodic boundary
conditions in both directions. Finite-$L$ and finite-$T$ effects on the
spectral density are significantly different. On the one hand, finite-$T$ (or
thermal) effects are shown below to be exponentially suppressed and are
therefore reliably small. On the other hand, the spectral density is
dramatically different at finite and infinite $L$.

At infinite $L$ the spectral density $\rho(E)$ defined in Eq. (7) is a
continuous function which is analytic for all $E$ except for $E=n\,m$ where
$n$ is a positive even integer. At finite $L$, by contrast, the spectrum of
the Hamiltonian is discrete and the spectral density is a sum of Dirac
$\delta$-functions. As $L$ increases, the spectrum becomes denser and denser
so that the continuous $L=\infty$ spectral density is recovered as a weak (or
distributional) limit. In no meaningful way can this limit be considered
point-wise and the finite-$L$ effects on $\rho(E)$ treated as small
corrections. Questions such as “Are the finite-$L$ corrections to $\rho(E)$
exponentially suppressed?” are simply ill-posed. By contrast, smeared spectral
densities converge to their infinite-$L$ value in a point-wise sense. An
immediate consequence of this discussion is that the $L\to\infty$ and
$\epsilon\to 0^{+}$ limits (where $\epsilon$ is the width of the smearing
kernel) do not commute and make sense only in a precise order: the
$L\to\infty$ limit must be taken before the $\epsilon\to 0^{+}$ limit. The
importance of the ordered double limit was also stressed in Ref.
Hansen:2017mnd  where the combined $L$ and $\epsilon$ dependence of related
spectral functions was considered in a perturbative context.

We know that individual finite-$L$ matrix elements at energies above the two-
particle threshold approach their infinite-$L$ limit with corrections which
vanish as inverse powers of $1/L$. When smeared spectral densities are
considered, one may hope that finite-$L$ corrections vanish faster than any
inverse power in $1/L$, at least if the smearing kernel is ‘reasonable
enough’. This issue is explored in App. C, which considers a fictitious system
where the infinite-$L$ spectral density is given entirely by the two-particle
contribution in Eq. (9). The finite-$L$ spectral density is then determined
(up to exponential corrections) by the Lellouch-Lüscher formalism
Lellouch:2000pv ; Meyer:2011um . One sees explicitly that the smeared spectral
density has finite-$L$ contributions which are $O({\rm e}^{-mL})$ in the case
of the Gaussian kernel, despite the power-law corrections to the individual
energies and matrix elements. For the Cauchy kernels the effects are also
exponentially suppressed by $O({\rm e}^{-\mu L})$ where $\mu\leq m$ depends on
the values of $E,\epsilon$ and $m$ as specified in Eq. (93).

In the derivation of the exponentially-suppressed finite-$L$ contributions,
the analyticity of the smearing kernel plays a central role. Smooth but non-
analytic kernels produce finite-$L$ corrections that decay faster than any
inverse power of $L$ but generally not exponentially. We also observe that the
finite-volume effects are oscillatory functions of $L$ with an envelope
decaying according to the predicted scaling. We do not claim that these
findings are valid beyond the simple exercise considered here, but it is clear
that the landscape of phenomena related to the $L\to\infty$ limit of smeared
spectral densities is rich.

After this lengthy introduction, we give some explicit formulae. At finite $L$
and $T$, the zero-momentum current-current correlator has the following
Hamiltonian representation

$\displaystyle C_{T,L}(t)=\frac{1}{L}\frac{\text{tr}\,\\{{\rm
e}^{-(T-t)\hat{H}_{L}}\hat{A}{\rm
e}^{-t\hat{H}_{L}}\hat{A}\\}}{\text{tr}\,{\rm e}^{-T\hat{H}_{L}}}\ ,$ (23)

where the definition

$\displaystyle\hat{A}=\int_{0}^{L}d\bm{x}\,\hat{j}_{1}(\bm{x})\,,$ (24)

is employed. At finite $L$, the spectrum of the Hamiltonian is discrete. By
introducing an orthonormal basis of energy eigenstates $|n\rangle_{L}$
satisfying $\hat{H}_{L}|n\rangle_{L}=E_{n}(L)|n\rangle_{L}$, one easily
derives the spectral representation of the correlator

$\displaystyle
C_{T,L}(t)=\int_{-\infty}^{\infty}d\omega\,\tilde{\rho}_{T,L}(\omega)e^{-t\omega}\,,$
(25)

with the definition

$\displaystyle\tilde{\rho}_{T,L}(\omega)=\frac{1}{L}\sum_{n,n^{\prime}}\frac{e^{-TE_{n^{\prime}}(L)}}{\sum_{n^{\prime\prime}}e^{-TE_{n^{\prime\prime}}(L)}}\left|\langle
n^{\prime}|\hat{A}|n\rangle_{L}\right|^{2}\delta(E_{n}(L)-E_{n^{\prime}}(L)-\omega)\,.$
(26)

Notice that $\tilde{\rho}_{T,L}(\omega)$ is non-vanishing also for negative
values of $\omega$. By separating terms with $E_{n}\geq E_{n^{\prime}}$ and
$E_{n}\leq E_{n^{\prime}}$, and taking care to avoid double counting
contributions with $E_{n}=E_{n^{\prime}}$, one can split
$\tilde{\rho}_{T,L}(\omega)$ into

$\displaystyle\tilde{\rho}_{T,L}(\omega)=\rho_{T,L}(\omega)+e^{T\omega}\rho_{T,L}(-\omega)\
,$ (27)

where the spectral density

$\displaystyle\rho_{T,L}(E)=\frac{1}{L}\sum_{E_{n}\geq
E_{n^{\prime}}}\frac{e^{-TE_{n^{\prime}}(L)}}{\sum_{n^{\prime\prime}}e^{-TE_{n^{\prime\prime}}(L)}}\left(1-\frac{\delta_{E_{n},E_{n^{\prime}}}}{2}\right)\left|\langle
n^{\prime}|\hat{A}|n\rangle_{L}\right|^{2}\delta(E_{n}(L)-E_{n^{\prime}}(L)-E)\,,$
(28)

vanishes for $E<0$. By plugging this information in the correlator, we obtain
the spectral representation in the form

$\displaystyle
C_{T,L}(t)=\int_{0^{-}}^{\infty}dE\,\rho_{T,L}(E)\\{e^{-tE}+e^{-(T-t)E}\\}\,.$
(29)

As described in Sec. II, for a target smearing kernel
$\delta_{\epsilon}(E-\omega)$, one now seeks an approximation
$\delta^{\lambda}_{\epsilon}(E,\omega)$ in the space generated by the function
basis

$\displaystyle b_{T,\tau}(\omega)=e^{-a\tau\omega}+e^{-(T-a\tau)\omega}\,,$
(30)

by minimizing the functional $W_{\lambda}[g]$. Then

$\displaystyle{\rho}^{\lambda}_{T,L,\epsilon}(E)=a\sum_{\tau=1}^{\tau_{\rm
max}}g^{\lambda}_{\tau}C_{T,L}(a\tau)=\int_{0^{-}}^{\infty}dE\,\rho_{T,L}(E){\delta}^{\lambda}_{\epsilon}(E,\omega)$
(31)

is an approximation for the smeared spectral density

$\displaystyle{\rho}_{T,L,\epsilon}(E)=\int_{0^{-}}^{\infty}dE\,\rho_{T,L}(E){\delta}_{\epsilon}(E-\omega)\,.$
(32)

A couple of comments on these formulae are in order. (1) The estimator in Eq.
(31) formally depends on $\tau_{\rm max}$, but (as demonstrated in Sec. V) our
results are rather insensitive to its particular value. (2) The coefficients
$g_{\tau}^{\lambda}$ obtained by solving the minimization problem in Eq. (16)
depend on $L$ and $T$ via the covariance matrix appearing in $B[g]$, and on
$T$ via our choice of the basis $b_{T,\tau}(\omega)$. In the limit of infinite
statistics (with $B[g]=0$), the $g_{\tau}^{\lambda}$ depend on $T$ but not on
$L$.

Eq. (28) makes manifest that finite-temperature effects on the spectral
density at fixed $L$ are exponentially suppressed, i.e.

$\displaystyle\rho_{T,L}(E)=\rho_{\infty,L}(E)+O({\rm
e}^{-Tm(L)})=\frac{1}{L}\sum_{n}\left|\langle
0|\hat{A}|n\rangle_{L}\right|^{2}\delta(E_{n}(L)-E)+O({\rm e}^{-Tm(L)})\,,$
(33)

where $m(L)$ is the energy gap in finite volume. Smearing the spectral density
amounts to replacing the Dirac $\delta$-function in Eq. (28) with the smearing
kernel $\delta_{\epsilon}$. At fixed $L$ the smeared spectral density also has
finite-temperature effects that are exponentially suppressed with $O({\rm
e}^{-Tm(L)})$. Even though far from obvious, it is reasonable to assume that
finite-$T$ effects are exponentially suppressed also at $L=\infty$ since the
theory has a mass gap.

##  IV Lattice-discretized model

We calculate the Euclidean current-current correlator in Eq. (29) by numerical
simulation of the lattice-discretized model. The field $\sigma(x)$ is now
defined on the set of $(L/a)\times(T/a)$ equally-spaced points
$x=(a\tau,\bm{x})$ on the two-torus, denoted by $\Lambda$. The standard
discretization of the action is employed here

$\displaystyle
S[\sigma]=\frac{\beta}{2}\sum_{x\in\Lambda}a^{2}\sum_{\mu}\hat{\partial}_{\mu}\sigma(x)\cdot\hat{\partial}_{\mu}\sigma(x)=\beta\sum_{x\in\Lambda}\sum_{\mu}\left[1-\sigma(x)\cdot\sigma(x+a\hat{\mu})\right]\,,$
(34)

where $\hat{\partial}_{\mu}f(x)=\frac{1}{a}[f(x+a\hat{\mu})-f(x)]$ is the
forward-difference operator defined with periodic boundary conditions in both
directions. Expectation values of observables are given by the path integral

$\displaystyle\langle A\rangle=\frac{1}{Z}\int\left[\prod_{x\in\Lambda}{\rm
d}\sigma(x)\right]\,{\rm e}^{-S[\sigma]}A[\sigma]\,,$ (35)

where ${\rm d}\sigma(x)$ denotes the O(3)-invariant integration measure over
the unit sphere. The global O(3) symmetry is preserved by the lattice
discretization, and its associated Noether current is

$\displaystyle
j_{\mu}^{a}(x)=\beta\epsilon^{abc}\sigma^{b}(x)\hat{\partial}_{\mu}\sigma^{c}(x)\,.$
(36)

Exact O(3) Ward identities imply that this current is renormalized. The
discretized version of the zero-momentum Euclidean current-current correlator
is given by

$\displaystyle C(a\tau)=\sum_{\bm{x}}a\langle
j_{1}^{a}(x)\,j_{1}^{a}(0)\rangle\,.$ (37)

The strong dynamics of the O(3) model generate a mass gap $m$, which is often
used to set the scale. However, in this work we determine all dimensionful
quantities in units of the appropriate power of the alternative scale
$m_{\star}$, defined as Francis:2013jfa

$\displaystyle m_{\star}^{-1}C(m_{\star}^{-1})=0.046615\,.$ (38)

While in principle one can choose any number in the right-hand side, we use a
value that gives $m_{\star}\simeq m$. Using the correlator reconstructed from
the $2$-, $4$\- and $6$-particle contributions to the spectral density we find
that the relative difference between $m$ and $m_{\star}$ is of order
$10^{-5}$. In practice the correlator $C(t)$ is known only at values of
$t=a\tau$ that are integer multiples of the lattice spacing $a$, and the
equation for $m_{\star}$ is solved using a piecewise linear interpolation of
$\log tC(t)$. The quantity $m_{\star}$ is determined with higher statistical
precision than $m$ and is also affected by smaller finite-volume effects. It
furthermore does not require any additional correlation functions.

ID | $(L/a)\times(T/a)$ | $\beta$ | $am_{\star}$ | $m_{\star}L$ | $m_{\star}T$ | $N_{\rm th}\times 10^{-6}$ | $N_{\rm rep}$ | $N_{\rm bin}$ | $B\times 10^{-6}$ | $N_{\rm c}\times 10^{-10}$
---|---|---|---|---|---|---|---|---|---|---
A1 | $640\times 320$ | 1.63 | 0.0447989(62) | 29 | 14 | 1.6 | 480 | 5 | 3.2 | 3.84
A2 | $1280\times 640$ | 1.72 | 0.0257695(31) | 33 | 17 | 12.8 | 3840 | 1 | 20 | 7.6
A3 | $1920\times 960$ | 1.78 | 0.0176104(31) | 34 | 17 | 12.8 | 7680 | 1 | 10 | 7.6
A4 | $2880\times 1440$ | 1.85 | 0.0112608(29) | 32 | 16 | 12.8 | 30720 | 1 | 2.5 | 7.68
B1 | $5760\times 1440$ | 1.85 | 0.0112607(73) | 65 | 16 | 12.8 | 7680 | 1 | 1.5 | 1.152
B2 | $2880\times 2880$ | 1.85 | 0.0112462(72) | 32 | 32 | 12.8 | 7680 | 1 | 1.5 | 1.152
Table 2:  Details for the ensembles of field configurations generated for this
work. The dimensionful scale $m_{\star}\simeq m$ is defined in Eq. (38).
Ensembles A1-A4 enable the continuum limit at approximately fixed physical
volume, while B1 and B2 are used to estimate finite size effects. Each
ensemble consists of $N_{\rm rep}$ independent identical replica, each of
which is thermalized for $N_{\rm th}$ cluster updates before $N_{\rm bin}$
measurements are taken by averaging over $B$ subsequent updates. This results
in a total number of measurement cluster updates $N_{\rm c}$ for ensemble.
Only two digits are given for $m_{\star}L$ and $m_{\star}T$.

The numerical simulations are performed with the single-cluster algorithm
described in App. D. Details concerning the generated ensembles of field
configurations are summarized in Tab. 2. The ensembles A1, A2, A3, A4 have
different values of $\beta$ and therefore different values of the lattice
spacing, but similar $m_{\star}L\geq 29$ and $m_{\star}T\geq 14$. The
ensembles B1 and B2 have been generated with the same lattice spacing as A4
but with doubled spatial and temporal extent, respectively. While the
ensembles A1, A2, A3, and A4 are used to perform a continuum extrapolation, B1
and B2 enable estimates of the residual finite-$L$ and finite-$T$ effects, as
explained in Sec. V.

The standard discretization of the two-dimensional O(3) $\sigma$-model
employed here is known to approach the continuum limit rather slowly. Ref.
Hasenbusch:2001ht  observed that lattice artifacts behave like $O(a)$ over a
large range of lattice spacings, in apparent contradiction with Symanzik’s
effective theory, which predicts an asymptotic $O(a^{2})$ behavior up to
logarithms. The puzzle was solved in Refs. Balog:2009yj ; Balog:2009np : the
asymptotic behavior is correctly described by Symanzik’s effective theory, but
the logarithmic corrections turn out to be large and must be included in
fitting formulae used to extrapolate to the continuum limit. On-shell
quantities, such as the mass gap or energy levels in finite volume, have the
asymptotic expansion

$\displaystyle
Q(a)=Q(0)+Ca^{2}\beta^{3}\left[1+\sum_{k=1}^{\infty}c_{k}\beta^{-k}\right]+O(a^{4})\,,$
(39)

where $c_{1}$ and $c_{2}$ are universal and calculable in terms of two-loop
and three-loop integrals respectively ($c_{1}=-1.1386\dots$ is analytically
known and $c_{2}=-0.4881$ is estimated numerically in Ref. Balog:2009np ),
while the other constants are non-perturbative. The leading logarithm
$\beta^{3}$ is generated by the dimension-four operator with largest one-loop
anomalous dimension appearing in the Symanzik expansion of the action. The
current two-point function, and consequently its spectral density, get extra
contributions from dimension-three operators appearing in the Symanzik
expansion of the current. These one-loop anomalous dimensions are unknown, and
their calculation is well beyond the scope of this work. For the continuum
extrapolation of the smeared spectral densities we use a fit function of the
type

$\displaystyle Q(a)=Q(0)+Ca^{2}\beta^{r}\,,$ (40)

where the exponent $r$ is fixed to $0,3,6$, and $Q(0)$ and $C$ are fit
parameters. We take $Q(0)$ at $r=3$ as our continuum extrapolation and the
spread generated by the three values of $r$ as an estimate of the systematic
error. One may argue that, since a $\beta^{3}$ term exists in the Symanzik
expansion, the correct asymptotic formula should have $r\geq 3$. However it is
conceivable that the coefficients of the various logarithms conspire in such a
way that an effective power $r<3$ is generated in some intermediate regime. We
therefore choose to include also $r=0$ in our analysis. As is evident in Fig
1, the variation of $r$ across $0$, $3$, and $6$ has little effect on the
continuum extrapolation of $tC(t)$ with $t=0.5m_{\star}^{-1}$ using our three
finest lattice spacings. These conclusions are supported by additional
extrapolations at several values of $t$ in the range
$t\in[0.5m_{\star},4m_{\star}]$. This strategy is thus adopted in Sec. V for
continuum extrapolations of $\rho_{\epsilon}(E)$, where the fourth lattice
spacing is included for stability in the presence of larger statistical errors
on $\rho_{\epsilon}(E)$.

Figure 1:  Illustration of the continuum limit procedure for
$tC(t)|_{m_{\star}t=0.5}$. Three different extrapolation forms are employed
using Eq. (40) with $r=0,3,6$. Shown is the difference between the numerical
estimate and the exact result divided by the exact result, denoted $\Delta$.
The horizontal band indicates our estimate for the continuum extrapolated
$\Delta$, given by the intercept from the fit with $r=3$. The statistical
error on this estimate is similar to the variation across the three values of
$r$, which are taken as an estimate of the systematic error. This validates
the use of the ‘phenomenological’ extrapolation forms in Eq. (40) for
$\rho_{\epsilon}(E)$ where the statistical errors are larger.

##  V Numerical Results

After detailing estimates of the systematic errors due to the reconstruction
and finite volume, this section presents the numerical verification of the
spectral reconstruction procedure using two different tests. The first test
(discussed in Sec. V.1) compares $\rho_{\epsilon}(E)$ with exact results at
fixed smearing width $\epsilon$ for each of the four smearing kernels in Eq.
(19). As anticipated, $\rho_{\epsilon}(E)$ is more difficult to determine with
increasing $E$ and decreasing $\epsilon$. The second test, which is detailed
in Sec. V.2, uses $\rho_{\epsilon}(E)$ at finite $\epsilon$ to extrapolate to
the $\epsilon\rightarrow 0$ limit and obtain the unsmeared spectral density
$\rho(E)$. The extrapolation procedure (as applied here) requires smearing
widths that are small compared to the scale at which the unsmeared spectral
density varies, and is therefore less effective in the elastic region where
$\rho(E)$ increases rapidly with the onset of two-particle phase space. Since
$\rho(E)$ in the elastic region is accessible to the finite-volume formalism,
the focus is instead on energies in the inelastic region $E>4m$. Due to the
absence of any sharp ‘resonance peaks’, $\rho(E)$ varies increasing slowly
with increasing $E$, suggesting that the smearing width should be scaled
$\epsilon\propto(E-2m)$, i.e. with the distance to the rapid variation from
the two-particle threshold. Apart from this scaling of the smearing width, the
analysis for these two tests proceeds identically.

All statistical errors on the results presented here are estimated using the
bootstrap procedure with $N_{\rm boot}=800$ samples, which are generated
independently on each ensemble from the $N_{\rm rep}\times N_{\rm bin}$
measurements listed in Tab. 2. For all reconstructions the lower bound of the
integration range $E_{0}$ defined in Eq. (15) is set using the scale
$m_{\star}$ via $aE_{0}=2\,am_{\star}$ to (approximately) coincide with the
two-particle threshold. This choice minimizes the range in energy over which
the functional $A[g]$ forces the reconstructed and target kernels to be
similar. The values of $am_{\star}$ given in Tab. 2 are also used to fix the
dimensionful parameters $\epsilon$ and $E$ at different lattice spacings.
Properly including the statistical error on $am_{\star}$ requires the
determination of the coefficients $g_{\tau}^{\lambda}$ on each bootstrap
sample, which significantly increases the computational cost. After confirming
that this has no observable effect on a selection of sample reconstructions,
the statistical error on the values of $am_{\star}$ in Tab. 2 is subsequently
ignored.

###  V.1 Fixed smearing width

Before presenting the main results of this section in Fig. 6, several
systematic errors must be estimated. Consider first systematic errors due to
the reconstruction procedure. Due to the finite number of input values
$C(a\tau)$, the reconstructed smearing kernel will never be exactly equal to
the desired one. This source of systematic error is quantified by the
functional $A[g^{\lambda}]/A[0]$ defined in Eq. (15). While this measure does
not take into account the role of the unsmeared spectral density $\rho(E)$ on
the systematic error of $\rho_{\epsilon}(E)$, this error certainly vanishes in
the $A[g^{\lambda}]\rightarrow 0$ limit. The value of $A[g^{\lambda}]$ is
therefore a useful diagnostic for determining the onset of the statistics-
limited regime: if $\lambda$ is lowered such that $A[g^{\lambda}]$ changes
significantly and no significant change is observed in
$\rho^{\lambda}_{\epsilon}(E)$, then this source of systematic error is likely
smaller than the statistical error. This type of ‘plateau’ analysis is well-
known to lattice field theory practitioners and is exemplified in Fig. 2. As
discussed in Sec. II, we follow here the recipe of Ref. Hansen:2019idp  and
quote as the best estimate for $\rho_{\epsilon}(E)$ the result from the
unconstrained reconstruction obtained at $\lambda=\lambda_{\star}$ defined in
Eq. (18).

Figure 2:  For the A1 ensemble detailed in Tab. 2, an illustration of the
compromise between statistical and systematic errors which results from
choosing $\lambda=\lambda_{\star}$. Left: the estimator
$\rho^{\lambda}_{\epsilon}(E)$ for the smeared spectral function at different
values of $\lambda$, plotted against the ‘resolution’ functional
$A[g^{\lambda}]/A[0]$ defined in Eq. (15). The horizontal band indicates the
result at $\lambda_{\star}$ (where $A[g^{\star}]=A[0]B[g^{\star}]$) without
any constraints on the determination of the optimal coefficients. This choice
is consistent with the one obtained by imposing the equal area constraint
(denoted ‘equal area’) as well as with the one demanding that reconstructed
smearing kernel take the correct value at the peak, denoted ‘equal peak
value’. Center: the reconstructed smearing kernel
$\delta_{\epsilon}^{\lambda}(E,\omega)$ for the different constraints together
with the target kernel $\delta_{\epsilon}^{\sf g}(E-\omega)$ shown as a solid
line. The value of $\lambda$ for each of the three reconstructions results in
$A[g^{\lambda}]/A[0]\approx 0.1$ and the vertical shaded band indicates
energies $\omega<E_{0}$ where the reconstructed kernels are left
unconstrained. Right: same as the center panel but with
$\lambda=\lambda_{\star}$.

Another probe of the systematic error due to the reconstruction is the
comparison of $\rho_{\epsilon}(E)$ determined from coefficients subject to
various constraints on the minimization, implemented as explained in App. B.
Fig. 2 shows results for the unconstrained reconstruction, the constraint that
the reconstructed smearing kernel has a signed area equal to the desired one,
defined in Eq. (60), and the constraint that the reconstructed kernel exactly
coincides with the desired one at the peak $\omega=E$, defined in Eq. (61).
For large $\lambda$, these different reconstructions differ
significantly999The equal area constraint is apparently somewhat weaker than
demanding $\delta_{\epsilon}^{\lambda}(E,\omega)=\delta_{\epsilon}^{\sf
x}(E-\omega)$ at $\omega=E$, as evidenced by the plots of
$\delta_{\epsilon}^{\lambda}(E,\omega)$ in the center and left panels of Fig.
2. but coincide well within the statistical error at $\lambda_{\star}$,
lending additional confidence that $\lambda_{\star}$ is indeed in the
statistics-dominated regime. Put more precisely, reconstructions subject to
different constraints result in different reconstructed smearing kernels.
Their consistency when compared at similar values of $A[g^{\lambda}]/A[0]$
demonstrates that these differences are insignificant and suggests that
deviations from the exact kernel are as well. Overall, the approach employed
for choosing $\lambda_{\star}$ discussed in Eq. (18) to balance statistical
and systematic errors is somewhat conservative. Fig 2 suggests the less
restrictive alternative strategy of merely demanding consistency between the
three different reconstructions. This could potentially result in smaller
statistical errors, but requires further investigation beyond the scope of
this work.

Figure 3:  Left: test of the sensitivity of the reconstruction procedure to
the range of time slices $[\tau_{\rm min},\tau_{\rm max}]$ by varying
$\tau_{\sf max}$, plotted using $t_{\rm max}m_{\star}=\tau_{\sf
max}(am_{\star})$ with fixed $\tau_{\rm min}=1$. Right: test of the
sensitivity to the spacing between time slices included in the reconstruction
$\sigma_{t}$. Both tests are performed on the A1 ensemble detailed in Tab. 2.
Evidently the reconstruction procedure is relatively stable under the
variation of $\sigma_{t}$ and $\tau_{\rm max}$.

The input correlator time slices $\\{C(a\tau)\\}$ range over $[\tau_{\rm
min},\tau_{\rm max}]$ with $\tau_{\rm min}=1$ and $\tau_{\rm max}=160$ used as
the best estimate of $\rho_{\epsilon}(E)$ everywhere. They are furthermore
correlated and suffer from an exponential degradation of the signal-to-noise
ratio with increasing $\tau$ with a rate phenomenologically similar to $m$.
The sensitivity to the input correlator data is probed in Fig. 3 which
demonstrates that, for a sample reconstruction, $\rho_{\epsilon}(E)$ is
relatively insensitive to $\tau_{\rm max}$ and to a ‘thinning’ of the input
data where the time slices are separated by
$\sigma_{t}=a(\tau_{n+1}-\tau_{n})$. These observations can be plausibly
explained by the signal-to-noise degradation and correlation of the input
data. The role of $B[g]$ is to penalize coefficients $g_{\tau}^{\lambda}$
which would result in a large statistical error. This penalty naturally
disfavors input data at large $\tau$, effectively ‘turning off’ these time
slices and resulting in an insensitivity to $\tau_{\rm max}$. Similarly, the
correlation between the input time slices may be responsible for the
robustness to changing $\sigma_{t}$.

The preceding discussion demonstrates that the reconstruction procedure on a
particular ensemble of field configurations results in a statistics-dominated
estimator for $\rho_{\epsilon}(E)$. No further systematic error due to the
reconstruction is assigned. We turn now to systematic errors associated with
the finite extent of the lattice. Along the (approximate) line of constant
physics defined by the ensembles A1-A4, these errors are crudely estimated by
considering separately the differences between B1 and A4, and B2 and A4. These
differences

$\Delta_{\rm
L}(\epsilon,E)=\rho_{T,L,\epsilon}(E)-\rho_{T,2L,\epsilon}(E)\qquad{\rm
and}\qquad\Delta_{\rm
T}(\epsilon,E)=\rho_{T,L,\epsilon}(E)-\rho_{2T,L,\epsilon}(E)\,,$ (41)

are taken as estimators for the systematic errors due to the finite $L$ and
$T$ and are added (in absolute value) independently for each $\epsilon$, $E$,
and smearing kernel. A selection of these deviations are shown in Fig. 4,
illustrating that they are at most marginally significant. Nonetheless, this
systematic error is subsequently taken into account and represents the largest
of the systematic error estimates. Although Eq. (33) ensures that finite-$T$
effects are ‘reliably’ ${\rm O}({\rm e}^{-mT})$, we conservatively account for
both finite $L$ and $T$ with additional systematic errors as discussed above.
An increase of statistics on the B1 and B2 ensembles would more accurately
estimate these systematic errors and provide a more stringent test of the
error estimates due to the reconstruction, but is left for future work.

Figure 4: The differences $\Delta_{\rm L}$ and $\Delta_{\rm T}$ from Eq. (41).
The vertical axis is normalized by the error on the estimated difference and
can be interpreted as the statistical significance of each of the finite size
effects. The sum of the mean values $|\Delta_{\rm L}|+|\Delta_{\rm T}|$ is
taken as an estimate of the systematic error due to finite size effects
independently for each $E$, $\epsilon$, and smearing kernel.

In addition to finite size effects, systematic errors due to the lattice
spacing must be estimated as outlined in Sec. IV. The continuum limit is taken
with four lattice spacings using the ensembles A1-A4, but (as detailed in Sec.
IV) the asymptotic $a^{2}$ behavior is affected by large logarithmic
corrections and the ‘phenomenological’ extrapolation form of Eq. (40) is
employed with $r=0$, $3$, and $6$. The value at $r=3$ is taken as the best
estimate and the largest deviation between any two as an estimate of the
systematic error. A selection of the extrapolations for the three different
values of $r$ is shown in Fig. 5. Evidently, the difference in the
extrapolated value varies little across these three extrapolation forms.

Figure 5:  Illustration of the continuum extrapolation procedure from Eq. (40)
with $r=0$, $3$, and $6$ at $E=3m_{\star}$, $\epsilon=m_{\star}$, and the
Gaussian smearing kernel. The three panels show the $r=0$ ansatz (left), the
$r=3$ (center), and the $r=6$ (right). The points at $a=0$ are the continuum
extrapolated values and the horizontal band is the extrapolated value from
$\beta^{3}a^{2}$, which is taken as the best estimate. As in Fig. 1, there is
evidently little variation across the different extrapolation forms.

After estimating systematic errors due to the reconstruction, finite lattice
size, and finite lattice spacing, the results for $\rho_{\epsilon}(E)$ are at
last confronted with the exact values. While no further systematic error is
assigned to the reconstruction procedure, estimates of the remaining three
sources of systematic error (finite $L$, finite $T$, continuum limit) are
added naively. The total systematic error is then combined in quadrature with
the statistical error on the continuum limit estimator. The finite-$L$ and $T$
errors are typically the dominant source of systematic error, and are similar
in magnitude to the statistical error.

Figure 6:  Summary of the continuum extrapolated values with statistical and
systematic errors combined in quadrature as described in the text. Top left:
all four smearing kernels are plotted for a large smearing radius
$\epsilon=2m_{\star}$ together with the exact results including 2-, 4-, and
6-particle contributions shown as solid lines. Bottom left: the ‘pull’ between
the numerical and exact results for the same data. The variation of the points
indicates that the numerical and exact results are consistent at the
2-$\sigma$ level. Top and bottom right: same type of plots but for a variety
of smearing widths with the Gaussian smearing kernel.

A selection of these continuum extrapolated values are shown in Fig. 6. In the
left two panels all four kernels are compared at a large fixed
$\epsilon=2m_{\star}$ to accentuate the difference between them. However, all
of these results are consistent with the exact values (including 2-, 4-, and
6-particle contributions) within the combined statistical and systematic
errors. The right two panels show various $\epsilon$ for the Gaussian kernel
and are similarly consistent with the exact numbers. The increasing difficulty
of the reconstruction problem with increasing $E$ and decreasing $\epsilon$ is
apparent.

###  V.2 Extrapolation to zero smearing width

The continuum extrapolated results for $\rho_{\epsilon}(E)$ displayed in Fig.
6 demonstrate that the spectral reconstruction procedure yields quantitatively
accurate results within the quoted statistical and systematic errors. However,
the difficulty in reconstructing $\rho_{\epsilon}(E)$ at large $E$ and small
$\epsilon$ naively suggests that the approach is of little use in the
inelastic region, precisely where the finite-volume formalism is not yet
developed. However, the increasing smoothness of $\rho(E)$ with increasing $E$
can be exploited by scaling $\epsilon\propto(E-2m_{\star})$ in order to probe
larger energies. The scaled values of $\rho_{\epsilon}(E)$ are used to
extrapolate $\epsilon\rightarrow 0$ and determine the unsmeared spectral
density $\rho(E)$ using the small-$\epsilon$ expansion discussed in Sec. II.
The rapid variation of $\rho(E)$ in the elastic region inhibits the
application of this approach there, so we focus on $E\geq 4m_{\star}$101010The
opening of subsequent $2n$-particle thresholds introduces non-analyticities in
$\rho(E)$, but for the range of energies considered here $\rho(E)$ is
dominated by fewer-particle contributions and therefore approximately smooth..

Figure 7: The same as Fig. 5 for $E=40m_{\star}$,
$\epsilon=1.2(E-2m_{\star})$, and the Gaussian smearing kernel.

A natural concern is the growth of cutoff effects with increasing $E$ and
$\epsilon$, since $\rho_{\epsilon}(E)$ has increasingly significant
contributions from more energetic states. Fig. 7 shows the continuum
extrapolations of $\rho_{\epsilon}(E)$ at the largest energy considered here,
which is $E=40m_{\star}$, and $\epsilon=1.2(E-2m_{\star})$, which is the
largest $\epsilon$ used in the $\epsilon\rightarrow 0$ extrapolation. Although
the continuum limit is somewhat steeper than those in Fig. 5, there is still
little variation across the different values of $r$. However, it appears that
the statistical error on the estimator $\rho_{\epsilon}^{\lambda}(E)$
decreases more rapidly with $am_{\star}$ in Fig. 7 than in Fig. 5. This
phenomenon requires further study, but could be influenced by our choice of
basis vectors, namely the inclusion of all $\tau\in[1,160]$ at each lattice
spacing. The edge of the Brillouin zone may also play a role since
$1/(am_{\star})\simeq E/m_{\star}$ when $E\simeq 20m_{\star}$ for the A1
ensemble and $E\simeq 90m_{\star}$ for the A4. Although not displayed here,
the continuum-extrapolated values for the energies and smearing widths used in
this section are consistent with the exact results at a level similar to Fig.
6.

The small-$\epsilon$ behavior of $\rho_{\epsilon}(E)$ depends on the smearing
kernel111111For a fair comparison of the effectiveness of the four kernels in
Eq. (19), $\epsilon$ for the c2 is enlarged by the rescaling
$\epsilon\rightarrow\sqrt{3}\epsilon$ so that $w_{2}^{\sf c2}=1$ in Tab. 1.
This rescaling is not employed for the c2 data in Sec. V.1. As discussed in
Sec. II, for the four smearing kernels considered here, this dependence is
encoded in the coefficients $w_{k}^{\sf x}$, which are known analytically and
given in Tab. 1. The coefficients $a_{k}(E)$ in Eq. (21) however depend on
$\rho(E)$ only and are therefore identical across the different smearing
kernels. This means that an extrapolation ansatz which includes terms up to
$O(\epsilon^{n})$ contains $n+1$ fit parameters for each of the $a_{k}$, which
are constrained by data from all the smearing kernels. In this work we target
$a_{0}(E)=\rho(E)$ only: no additional information from the other $a_{k}(E)$
is used to further constrain the desired spectral density.

Figure 8:  Left: a sample constrained extrapolation of the g, c1, and c2
kernels at $E=14m_{\star}$ with the (constrained) extrapolation form including
terms up to $O(\epsilon^{4})$. Right: summary of the extrapolated values for
different extrapolations forms and different extrapolation ranges
$[\epsilon_{\rm min},\epsilon_{\rm max}]$ with $\epsilon_{\rm
min}=0.3(E-2m_{\star})$ fixed and varying $\epsilon_{\rm max}$. For each
$\epsilon_{\rm max}$, the extrapolated result is plotted with the correlated
$\chi^{2}/{\rm d.o.f.}$ as the horizontal coordinate. The horizontal band
shows the final estimate for the extrapolated values, with statistical and
systematic errors combined in quadrature. The systematic errors are estimated
as discussed in the text.

The clear sources of systematic error in this approach are the extrapolation
form and the extrapolation range $[\epsilon_{\rm min},\epsilon_{\rm max}]$.
Due to the difficulty in reconstructing small $\epsilon$, the statistical
errors increase with decreasing $\epsilon$ so that in practice varying
$\epsilon_{\rm min}$ has little effect on the extrapolated value. We therefore
fix $\epsilon_{\rm min}=0.3(E-2m_{\star})$ for all smearing kernels and
extrapolations. The effect of varying $\epsilon_{\rm max}$ must be monitored
more closely, however. In order to fairly include data from different smearing
kernels, $\epsilon_{\rm max}$ for each kernel is re-scaled $\epsilon_{\rm
max}\rightarrow\epsilon_{\rm max}/\alpha_{\sf x}$ with $\int_{-\alpha_{\sf
x}}^{\alpha_{\sf x}}{\rm d}x\,\delta^{\sf x}_{\epsilon=1}(x)={\rm
erf}(2^{-1/2})\approx 0.68269$ chosen to approximately coincide with the
Gaussian. This treats the data from different smearing kernels on equal
footing: all smearing widths are included up to those with the same amount of
‘leakage’ down to the two-particle threshold, where $\rho(E)$ varies rapidly.
For the c0 kernel the large value $\alpha_{\sf c0}\approx 1.84$ together with
the leading $O(\epsilon)$ behavior effectively renders it useless in the
extrapolations. For the Gaussian kernel $\alpha_{\sf g}=1$ by definition,
while $\alpha_{\sf c1}=0.7$ and $\alpha_{\sf c2}=0.855$ increase the fit range
for these two kernels relative to the Gaussian. The extrapolations are
therefore performed with the three kernels g, c1, and c2. An example of such
an extrapolation is shown in the left panel of Fig. 8.

The systematic error estimate due to the $\epsilon\rightarrow 0$ extrapolation
proceeds by employing extrapolation forms up to and including
$O(\epsilon^{p})$ with $p=2$, $3$, and $4$ and varying $\epsilon_{\rm max}$
keeping below $\epsilon_{\rm max}<1.3(E-2m_{\star})$. For each value of $p$,
the largest $\epsilon_{\rm max}$ resulting in a (correlated) $\chi^{2}/{\rm
d.o.f.}<1$ is identified, with the $p=4$ value taken as the best estimate for
$\rho(E)$. The spread of these three values is then added in quadrature as a
systematic error. An example plot showing the consistency between different
extrapolation forms and extrapolation ranges is shown in the right panel of
Fig. 8. This procedure is performed for all energies for which $E/m_{\star}$
takes integral values from $4$ to $40$. A summary of the extrapolated results
is shown in Fig. 9, which are consistent with the exact results including 2,
4, and 6 particle contributions. Interestingly, the data exhibits some
sensitivity to $\rho^{(4)}(E)$ for $E\gtrsim 20m_{\star}$. Four-particle
scattering amplitudes are currently beyond the reach of the finite-volume
formalism.

Figure 9: A summary of the continuum, $\epsilon\rightarrow 0$ extrapolated
results for $\rho(E)$, together with the exact two-particle contribution
(light dashed line), the two-, four-, and six-particle contributions (dark
solid line), and the 2-loop perturbative result (dark dotted line).
Statistical and systematic errors due to the finite volume, continuum limit,
and $\epsilon\rightarrow 0$ extrapolation are combined in quadrature as
described in the text.

##  VI Conclusions

The aim of the preceding sections is to verify the procedure of Ref.
Hansen:2019idp  for numerically computing smeared spectral densities (with an
_a priori_ specified smearing kernel) from lattice field theory correlation
functions. In this regard the two-dimensional O(3) model usefully provides
exact results against which the estimates can be checked. These checks, which
are presented in Figs. 6 and 9, are satisfied and compare both
$\rho_{\epsilon}(E)$ at finite $\epsilon$ and the results from
$\epsilon\rightarrow 0$ extrapolations to determine $\rho(E)$ deep into the
inelastic region where finite-volume methods have not yet been developed. The
highest energy considered here is $E=40m_{\star}$, at which $\rho(E)$ is
determined with a relative accuracy of $5\%$ and differs significantly from
the exact two-particle contribution $\rho^{(2)}(E)$ given in Eq. (9).

Apart from the ‘usual’ sources of systematic error due to the finite lattice
spacing and finite-volume spacetime, we must also consider the imperfect
reconstruction of the smearing kernel due to the finite number of input time
slices and their associated statistical errors. All sources of systematic
error have been estimated and included in Figs. 6 and 9 where the statistical
and systematic errors are added in quadrature. Generally the errors due to the
finite lattice extent are the largest source of systematic uncertainty, and
are typically less than or comparable to the statistical errors.

The determination of $\rho_{\epsilon}(E)$ becomes increasingly difficult for
smaller smearing widths $\epsilon$ at fixed energy $E$, and increasing $E$
with fixed $\epsilon$. As is evident from the right two panels of Fig. 6, it
is difficult to achieve precise results outside of the elastic region for
$\epsilon\lesssim m/2$ with the Gaussian smearing kernel. Better is to exploit
the smoothness of $\rho(E)$ and scale $\epsilon\propto(E-2m)$, so that an
equal proportion of the smearing kernel ‘leaks’ down to the two particle
threshold at each energy. This enables the determination of $\rho(E)$ in Fig.
9, which is the main result of this work.

The analysis performed here is (in principle) directly applicable to the
vector-vector correlator in QCD to compute the $R$-ratio, provided a few key
points are addressed. First, the presence of narrow resonances in QCD
seemingly invalidates our approach to the $\epsilon\rightarrow 0$
extrapolation by scaling $\epsilon\propto(E-2m)$. Nonetheless, the finite
smearing width may be used to probe the onset of the perturbative regime
Poggio:1975af . Furthermore the spatial extent achieved here of $mL=30$ is
currently beyond the state of the art for simulations at physical quark
masses. To this point it is important to note that the density of states is
higher in three spatial dimensions compared to one so that smaller values of
$mL$ may be acceptable, though still likely larger than the typical lattice
QCD volumes. In this respect, the masterfield simulation paradigm
Francis:2019muy ; Giusti:2018cmp  may enable larger volumes in the near future
Ce:2021akh .

Apart from this, the difficulty in the reconstruction problem faced here may
in fact be similar to QCD. Although the use of $\tau_{\rm max}=160$ may seem
daunting, Fig. 3 indicates that reducing the range and density of input
correlator time slices likely has little effect. Although the input correlator
data for this work is generated using the two-cluster algorithm, it still
possesses an exponentially degrading signal-to-noise ratio which decays at a
rate empirically similar to $m$, as expected for the vector-vector correlator
in QCD. Also, as is typical with scalar field theory, the statistical errors
achieved at $C(m_{\pi}^{-1})$ are roughly comparable with state-of-the-art
determinations of the vector-vector correlator in QCD.

Finally, the validation of the spectral reconstruction procedure presented
here is intended as a stepping stone to other applications. Among them is the
determination of exclusive scattering amplitudes by exploiting an asymptotic
formalism as in Refs. Bulava:2019kbi ; Bruno:2020kyl . There are several
additional challenges there compared to this work, including the computation
and spectral reconstruction of at least three-point temporal correlation
functions with two associated time separations. Nonetheless, this may provide
an alternative to the finite-volume formalism for determining few-body
scattering amplitudes where applicable, and extend the reach of such
computations by enabling (in principle) arbitrary center-of-mass energies.

###### Acknowledgements.

We thank Martin R. Hansen for his participation at an early stage of this work
and acknowledge helpful conversations with Rainer Sommer and Mattia Bruno, and
with Peter Weisz, who also provided comments on a previous version of this
manuscript. We also thank the organizers of the CERN Workshop “Advances in
Lattice Gauge Theory” where work on this project began. The work of MTH is
supported by UK Research and Innovation Future Leader Fellowship MR/T019956/1.
The research of AP is funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) - Projektnummer 417533893/GRK2575 “Rethinking
Quantum Field Theory”. This work was supported by CINECA that granted
computing resources on the Marconi supercomputer to the LQCD123 INFN
theoretical initiative under the CINECA-INFN agreement.

##  Appendix A Analytic expressions for the spectral density

Thanks to the integrability of the O(3) model, the matrix elements of the
Noether current $j_{\mu}$ between the vacuum and a generic $n$-particle state
(form factors) are completely determined by a known set of recursive
functional equations, derived by Karowski and Weisz Karowski:1978vz . These
functional equations are written in terms of the S-matrix of the system, which
has been analytically derived by Zamolodchikov and Zamolodchikov
Zamolodchikov:1977nu  on the basis of the integrability of the model and some
fairly weak assumptions. The $n$-particle contribution to the spectral density
of the Noether current can be reconstructed from the form factors. This
program has indeed been carried out by Balog and Niedermaier Balog:1996ey  for
$n=2,4,6$. Interestingly, the original motivation for this calculation was the
calculation of the Euclidean two-point function of the O(3) model, and the
comparison with lattice simulations and perturbative calculations Balog:1996mn
.

For convenience we summarize here the most important formulae taken from
Balog:1996ey . A different normalization for the spectral density with respect
to Balog and Niedermaier (BN) is employed here

$\displaystyle\rho(\mu)=\frac{3\mu}{2}\,\rho_{\text{BN}}(\mu)=\frac{3\mu}{2}\sum_{n=2,4,\dots}\rho_{\text{BN}}^{(n)}(\mu)\
.$ (42)

With this normalization, the large-$\mu$ behaviour of the spectral density is
given by $\rho(\mu)=1/(2\pi)+O(g^{2}(\mu))$.

The form factors are usually written in terms of the rapidities
$\theta_{i=1,\dots,n}$. The (spatial) momentum and energy of the $i$-th
particle are given by

$\displaystyle p_{i}=m\sinh\theta_{i}\ ,\qquad
E_{i}=E(p_{i})=m\cosh\theta_{i}\ .$ (43)

The invariant mass of a system of $n$ particles is then

$\displaystyle
M^{(n)}=\left\\{\left(\sum_{i=1}^{n}E_{i}\right)^{2}-\left(\sum_{i=1}^{n}p_{i}\right)^{2}\right\\}^{1/2}=m\left\\{n+2\sum_{1\leq
i<j\leq n}\cosh(\theta_{i}-\theta_{j})\right\\}^{1/2}\ ,$ (44)

where trivial hyperbolic function identities are applied. Notably the
invariant mass depends only on the rapidity differences, so that it is
convenient to define the new variables

$\displaystyle u_{i}=\theta_{i}-\theta_{i+1},\qquad{\rm for}\qquad
i=1,\dots,n-1\ .$ (45)

The kinematics of the $n$-particle system are completely specified by the set
of variables $\theta_{i=1,\dots,n}$, or equivalently by $\theta_{n}$ and
$u_{i=1,\dots,n-1}$. In terms of the new variables

$\displaystyle\theta_{i}-\theta_{j}=\sum_{k=i}^{j-1}u_{k}\ ,\qquad\text{for
}1\leq i<j\leq n\ .$ (46)

Once the phase-space integral appearing in the definition of the $n$-particle
contribution to the spectral density has been written in terms of the new
variables, one can use the $\delta$-function over the spatial momentum to
eliminate the integral over $\theta_{n}$. One is left with the integral over
the variables $u_{i=1,\dots,n-1}$ of a known function times the
$\delta$-function over the invariant mass

$\displaystyle\rho_{\text{BN}}^{(n)}(\mu)=\int_{0}^{\infty}\frac{du_{1}\cdots
du_{n-1}}{(4\pi)^{n-1}}F^{(n)}(u)\delta\left(\mu-M^{(n)}(u)\right)\ .$ (47)

The form factor is parametrized as

$\displaystyle F^{(n)}(u)=\frac{\pi^{3n-2}}{4}G^{(n)}(u)\prod_{1\leq i<j\leq
n}\left|\frac{\theta_{ij}-i\pi}{\theta_{ij}(2\pi
i-\theta_{ij})}\tanh^{2}\frac{\theta_{ij}}{2}\right|_{\theta_{ij}=\theta_{i}-\theta_{j}}\
,$ (48)

where the r.h.s. can be written as a function of $u$ by means of eq. (46), and
the functions $G^{(n)}(u)$ are polynomials in $u$ which are given explicitly
in Ref. Balog:1996ey  for $n=2,4,6$. We report here only

$\displaystyle G^{(2)}(u)=2\ ,$ (49) $\displaystyle
G^{(4)}(u)=-4[6\tau_{2}^{3}+9\tau_{3}^{2}+40\tau_{2}\tau_{4}]+8\pi^{2}[25\tau_{2}^{2}+44\tau_{4}]-448\pi^{2}\tau_{2}+272\pi^{6}\
,$ (50)

where, for the $n=4$ case, the following auxiliary variables have been
introduced

$\displaystyle\tau_{k}=\sum_{1\leq i_{1}<\dots<i_{k}\leq
4}\prod_{\ell=1}^{k}\left(\theta_{i_{\ell}}-\frac{1}{4}\sum_{j=1}^{4}\theta_{j}\right)\quad{\rm
for}\quad k=1,2,3,4.$ (51)

The corresponding expressions for $n=6$ are considerably longer, but present
no logical complication. In order to calculate the integral (47) numerically
one can solve the equation $M^{(n)}(u)=\mu>0$ for e.g. $u_{n-1}$. The only
positive solution $\bar{u}_{n-1}(u_{1},\dots,u_{n-2})$ can be written
explicitly by means of standard algebraic manipulations. The integral over
$u_{n-1}$ is then explicitly given by

$\displaystyle\rho_{\text{BN}}^{(n)}(\mu)=\int_{M^{(n)}(u_{1},\dots,u_{n-2},0)<\mu}\frac{du_{1}\cdots
du_{n-2}}{(4\pi)^{n-1}}\left\\{\left[\frac{\partial M^{(n)}}{\partial
u_{n-1}}(u)\right]^{-1}F^{(n)}(u)\right\\}_{u_{n-1}=\bar{u}_{n-1}(u_{1},\dots,u_{n-2})}\
.$ (52)

Notice that for $n=2$ there is no remaining integral, and one obtains the
closed-form expression

$\displaystyle\rho_{\text{BN}}^{(2)}(\mu)=\theta(\mu-2m)\frac{8\pi^{3}m^{5}}{\mu^{6}}\left(\frac{\mu^{2}}{4m^{2}}-1\right)^{5/2}\left[\frac{\bar{u}_{1}^{2}+\pi^{2}}{\bar{u}_{1}^{2}(4\pi^{2}+\bar{u}_{1}^{2})}\right]_{\bar{u}_{1}=2\cosh^{-1}\frac{\mu}{2m}}\
.$ (53)

For $n=4,6$, the integral in (52) is equivalently replaced by

$\displaystyle\int_{M^{(n)}(u_{1},\dots,u_{n-2},0)<\mu}\frac{du_{1}\cdots
du_{n-2}}{(4\pi)^{n-1}}\ \longrightarrow\
\int_{0}^{\cosh^{-1}\frac{\mu^{2}-n\,m}{2}}\frac{du_{1}\cdots
du_{n-2}}{(4\pi)^{n-1}}\,\theta\left(\mu-M^{(n)}(u_{1},\dots,u_{n-2},0)\right)\
,$ (54)

and is calculated numerically for $m=1$ (the $m$ dependence can be trivially
reintroduced by dimensional analysis) and a selection of values of $\mu$ by
means of the GNU Scientific Library implementation of Press and Farrar’s MISER
Monte Carlo algorithm Press:1989vk . In the region $\mu\leq 40m$, the
$6$-particle contribution to the spectral density is always smaller than the
$0.2\%$ of the sum of the $2$\- and $4$-particle contributions.

##  Appendix B Spectral reconstruction implementation

Explicit expressions and numerical implementation details of the spectral
reconstruction algorithm are provided here for completeness. The presentation
follows Ref. Hansen:2019idp  with a slightly different notation. Additional
information concerning the constraints used in this work is also provided.

By using a matrix-vector notation, the functionals in Eq. (15) can be
expressed as

$\displaystyle
A[g]=a\left\\{\bm{g}^{T}\cdot\bm{A}\cdot\bm{g}-2\bm{g}^{T}\cdot\bm{f}\right\\}+A[0]\;,\qquad
B[g]=\bm{g}^{T}\cdot\bm{B}\cdot\bm{g}\;,$ (55)

where $\bm{g}^{T}=(g_{1},\cdots,g_{\tau_{\rm max}})$ is the vector collecting
the coefficients $g_{\tau}$, $\bm{f}$ the vector with entries

$\displaystyle
f_{\tau}=\int_{E_{0}}^{\infty}d\omega\,{\delta}_{\epsilon}(E,\omega)\,b_{\tau}(\omega)\;,$
(56)

and $\bm{A}$ and $\bm{B}$ matrices with elements

$\displaystyle
A_{\tau\tau^{\prime}}=a\int_{E_{0}}^{\infty}d\omega\,b_{\tau}(\omega)\,b_{\tau^{\prime}}(\omega)\;,\qquad
B_{\tau\tau^{\prime}}={\rm Cov}\left[aC(a\tau),aC(a\tau^{\prime})\right]\;.$
(57)

Here $C(a\tau)$ is the correlator at time $\tau$ in lattice units and

$\displaystyle
A[0]=\int_{E_{0}}^{\infty}d\omega\,\left\\{{\delta}_{\epsilon}(E,\omega)\right\\}^{2}\;.$
(58)

With these definitions the vector $\bm{g^{\lambda}}$ that solves the minimum
conditions of Eq. (16) is given by

$\displaystyle\bm{g^{\lambda}}=(1-\lambda)\bm{W^{-1}_{\lambda}}\cdot\bm{f}\;,\qquad\bm{W_{\lambda}}=(1-\lambda)\bm{A}+\lambda\frac{A[0]}{a}\bm{B}\;.$
(59)

In Fig. 2 the minimization problem is subject to two different constraints on
the reconstructed kernel. The first is detailed in Ref. Hansen:2019idp  and
forces the reconstructed and target kernels to have the same area

$\displaystyle\int_{E_{0}}^{\infty}d\omega\,a\sum_{\tau=1}^{\tau_{\rm
max}}g_{\tau}\,b_{\tau}(\omega)=\int_{E_{0}}^{\infty}d\omega\,{\delta}_{\epsilon}(E,\omega)\;.$
(60)

The second constraint requires the reconstructed and target kernels to
coincide at $\omega=E$

$\displaystyle a\sum_{\tau=1}^{\tau_{\rm
max}}g_{\tau}\,b_{\tau}(E)={\delta}_{\epsilon}(E,E)\;.$ (61)

Both constraints are represented in vector notation as

$\displaystyle\bm{R}^{T}\cdot\bm{g}=r\;,$ (62)

where for the constraint in Eq. (60) the entries of the vector $\bm{R}$ and
the constant $r$ are given by

$\displaystyle
R_{\tau}=a\int_{E_{0}}^{\infty}d\omega\,b_{\tau}(\omega)\;,\qquad
r=\int_{E_{0}}^{\infty}d\omega\,{\delta}_{\epsilon}(E,\omega)\;,$ (63)

while for the constraint in Eq. (61) they are given by

$\displaystyle R_{\tau}=b_{\tau}(E)\;,\qquad
r=\frac{{\delta}_{\epsilon}(E,E)}{a}\;.$ (64)

The solution of the minimization problem subject to Eq. (62) is

$\displaystyle\bm{g^{\lambda}}=(1-\lambda)\bm{W^{-1}_{\lambda}}\cdot\bm{f}+\bm{W^{-1}_{\lambda}}\cdot\bm{R}\,\frac{r-(1-\lambda)\bm{R}^{T}\cdot\bm{W^{-1}_{\lambda}}\cdot\bm{f}}{\bm{R}^{T}\cdot\bm{W^{-1}_{\lambda}}\cdot\bm{R}}\;.$
(65)

For the basis functions actually used in this work, namely

$\displaystyle b_{T,\tau}(E)=e^{-aE\tau}+e^{-E(T-a\tau)}\;,$ (66)

explicit expressions for the entries of the matrix $\bm{A}$ are

$\displaystyle
A_{\tau\tau^{\prime}}=\frac{e^{-E_{0}(a\tau+a\tau^{\prime})}}{\tau+\tau^{\prime}}+\frac{e^{-E_{0}(T-a\tau+a\tau^{\prime})}}{T/a-\tau+\tau^{\prime}}+\frac{e^{-E_{0}(T-a\tau^{\prime}+a\tau)}}{T/a-\tau^{\prime}+\tau}+\frac{e^{-E_{0}(2T-a\tau-a\tau^{\prime})}}{2T/a-\tau-\tau^{\prime}}\;.$
(67)

When $\tau_{\rm max}$ is large this matrix is poorly conditioned.
Consequently, when both $\lambda$ and the errors on the input data are small,
the coefficients $g_{\tau}^{\lambda}$ are large in magnitude and oscillate in
sign at different values of $\tau$. In these situations the signal on the
reconstructed spectral density, given in vector notation by

$\displaystyle\rho^{\lambda}_{\epsilon}(E)=\bm{c}^{T}\cdot\bm{g^{\lambda}}\;,\qquad
c_{\tau}=aC(a\tau)\;,$ (68)

results from large numerical cancellations. For this reason (as in Ref.
Hansen:2019idp ), an implementation using extended-precision arithmetic is
advocated. This is relatively straightforward for the kernels considered in
this work using the libraries of Refs. Johansson2017arb ; MPFR  because the
integrals appearing in Eqs. (56), (58), and (63) can be expressed in terms of
standard special functions already implemented in these packages121212See for
example the publicly available code at https://github.com/mrlhansen/rmsd. In
the case of more generic kernels a numerical integration implemented in
extended-precision arithmetic may be required. An implementation based on the
double-exponential quadrature algorithm MORI2001287  is available upon
request.

##  Appendix C Finite-volume effects in a model smeared spectral density

In this appendix we present results concerning finite-volume effects in an
idealized system where the spectral density is given entirely by
$\rho^{(2)}(E)$ in Eq. (9) and its finite-volume analog contains only two-
particle states.131313Stated more precisely, this condition means that the
finite-volume smeared spectral density is equal to the one in Eq. (76). The
results are therefore relevant for $\rho_{\epsilon}(E)$ in the O(3) model for
values of $E$ and $\epsilon$ for which $n\geq 4$ particle states do not
significantly contribute. Using the notation of Sec. III, we first define
$\rho_{\infty,L,\epsilon}(E)$ as the full spectral density (with
$n=2,4,6,\cdots$ contributions) at finite $L$ and infinite $T$ smeared with a
resolution function. This can be written as

$\rho_{\infty,L,\epsilon}(E)=\sum_{n}c_{n}(L)\,\delta_{\epsilon}(E-E_{n}(L))\,,$
(69)

where $\delta_{\epsilon}(x)$ a generic smearing kernel, and

$c_{n}(L)=L\sum_{a}\big{|}\langle
0|\hat{j}^{a}_{1}(0)|n\rangle_{L}\big{|}^{2}\,.$ (70)

In absence of $n\geq 4$ particle states in the sum in Eq. (69) the $E_{n}(L)$
and $c_{n}(L)$ are given by the Lüscher quantization condition Luscher:1990ck
and the Lellouch-Lüscher relation Lellouch:2000pv , respectively141414See also
Ref. Briceno:2020rar  for explicit expressions relevant to a two-dimensional
spacetime.

$\displaystyle E_{n}(L)$ $\displaystyle=Q^{-1}(\pi n)+O({\rm e}^{-mL})\,,$
(71) $\displaystyle c_{n}(L)$ $\displaystyle=\pi\bigg{(}\frac{\partial
Q(E)}{\partial E}\bigg{)}^{-1}\rho^{(2)}(E)\bigg{|}_{E=E_{n}(L)}+O({\rm
e}^{-mL})\,,$ (72)

with

$Q(E)\equiv\delta_{I=1}(E)+\frac{L\sqrt{E^{2}/4-m^{2}}}{2}\,.$ (73)

This means that (up to corrections exponentially suppressed in $L$) both the
energies and overlaps are dictated by kinematic factors and the two-particle
scattering phase shift $\delta_{I=1}(E)$. The $I=1$ label emphasizes that this
scattering phase shift corresponds to two-particle states transforming
irreducibly under the fundamental representation of the global O(3) symmetry
and distinguishes the phase shift from the resolution function
$\delta_{\epsilon}(x)$ and the Dirac $\delta$-function. The phaseshift
$\delta_{I=1}$ can be written explicitly via the analytically-known
$S$-matrix, denoted $S(k)$

$\displaystyle\delta_{I=1}(E)$
$\displaystyle\equiv\frac{1}{2i}\log[S(k)]\bigg{|}_{k=\sqrt{E^{2}/4-m^{2}}},\qquad
S(k)=\frac{\theta+2i\pi}{\theta-2i\pi}\frac{\theta-i\pi}{\theta+i\pi}\bigg{|}_{\theta=2\sinh^{-1}\frac{k}{m}}\,,$
(74)

where the log is defined such that $0<\delta_{I=1}(E)<\pi$. Note that the
$\theta$ used here matches that of Eq. (9), since

$\frac{E^{2}}{4m^{2}}=1+\frac{k^{2}}{m^{2}}\ \Longrightarrow\
\sinh^{-1}\frac{k}{m}=\cosh^{-1}\frac{E}{2m}\,.$ (75)

The two-particle, finite-volume, smeared spectral density is therefore defined
as

$\rho^{(2)}_{L,\epsilon}(E)\equiv\sum_{n}\pi\bigg{(}\frac{\partial
Q(\omega)}{\partial\omega}\bigg{)}^{-1}\delta_{\epsilon}(E-\omega)\rho^{(2)}(\omega)\bigg{|}_{\omega=Q^{-1}(\pi
n)}\,,$ (76)

where we have dropped the $T=\infty$ label. If $\delta_{\epsilon}(x)$ has
compact support and $\epsilon$ and $E$ are chosen so that only $E_{n}(L)<4m$
states contribute, then this corresponds exactly to the smeared finite-volume
spectral density in the full theory. By contrast, for non-compact smearing or
for $E>4m$, $\rho^{(2)}_{L,\epsilon}(E)$ is an approximation in which the role
of infinite-volume states with more than two-particles is ignored.

In the remainder of this appendix, we show that terms scaling as $1/L$ cancel
in $\rho^{(2)}_{L,\epsilon}(E)$ provided that $\delta_{\epsilon}(x)$ is
differentiable and falls off fast enough as $x\to\infty$. The precise
condition is given after Eq. (86) below. We additionally show the stronger
result that $\rho^{(2)}_{L,\epsilon}(E)$ has exponentially suppressed volume
effects provided $\delta_{\epsilon}(x)$ is analytic in some finite-width strip
about the real axis. Before investigating the nature of the scaling at
asymptotically large $L$, we formally evaluate the $L\to\infty$ limit of Eq.
(76)

$\displaystyle\rho^{(2)}_{\infty,\epsilon}(E)$
$\displaystyle\equiv\lim_{L\to\infty}\rho^{(2)}_{L,\epsilon}(E)\,,$ (77)
$\displaystyle=\lim_{L\to\infty}\frac{1}{L}\sum_{n}\pi\bigg{(}\frac{1}{2}\frac{\partial\sqrt{\omega^{2}/4-m^{2}}}{\partial\omega}+\frac{1}{L}\frac{\partial\delta_{I=1}(\omega)}{\partial\omega}\bigg{)}^{-1}\delta_{\epsilon}(E-\omega)\rho^{(2)}(\omega)\bigg{|}_{\omega=Q^{-1}(\pi
n)}\,,$ (78)
$\displaystyle=\int_{0}^{\infty}\frac{dk}{2\pi}\,\pi\bigg{(}\frac{1}{2}\frac{\partial\sqrt{\omega^{2}/4-m^{2}}}{\partial\omega}\bigg{)}^{-1}\delta_{\epsilon}(E-\omega)\rho^{(2)}(\omega)\bigg{|}_{\omega=2\sqrt{m^{2}+k^{2}}}\,,$
(79)

where in the final line we have dropped contributions to the finite-volume
energy and the Lellouch-Lüscher factor that are $1/L$ suppressed, replaced the
sum over $n$ with an integral and performed the change of variables $k=2\pi
n/L$. These results hold for any smearing kernel $\delta_{\epsilon}(x)$ for
which the integral converges. Finally, changing integration variables to
$\omega$ yields

$\displaystyle\rho^{(2)}_{\infty,\epsilon}(E)$
$\displaystyle=\int_{2m}^{\infty}d\omega\,\delta_{\epsilon}(E-\omega)\rho^{(2)}(\omega)\,,$
(80)

as expected.

###  C.1 $1/L$ cancellation

We now show that $1/L$ contributions to the finite-volume energies and
Lellouch-Lüscher factors cancel in the definition of
$\rho^{(2)}_{L,\epsilon}(E)$. First write

$\rho^{(2)}_{L,\epsilon}(E)=\rho^{(2)}_{\infty,\epsilon}(E)+\frac{c^{(1)}(E)}{L}+\mathcal{O}(1/L^{2})\,,$
(81)

implying

$c^{(1)}(E)\equiv\lim_{L\to\infty}L\big{[}\rho^{(2)}_{L,\epsilon}(E)-\rho^{(2)}_{\infty,\epsilon}(E)\big{]}\,.$
(82)

Contributions to this coefficient arise from the leading-order shift to both
the matrix element and the energy. The former can be written as

$\bigg{(}\frac{\partial
Q(\omega)}{\partial\omega}\bigg{)}^{-1}=\frac{1}{L}\frac{8k}{\omega}\bigg{[}1-\frac{1}{L}\frac{8k}{\omega}\frac{\partial\delta_{I=1}(\omega)}{\partial\omega}\bigg{]}+\mathcal{O}(1/L^{3})\,,$
(83)

where we have substituted $\partial k/\partial\omega=\omega/(4k)$ for
$k=\sqrt{\omega^{2}/4-m^{2}}$. The corresponding result for the energy is
obtained by expanding Eq. (71)

$E_{n}(L)=2\sqrt{m^{2}+k^{2}}-\frac{1}{L}\frac{8k}{\omega}\delta_{I=1}(\omega)+\mathcal{O}(1/L^{2})\,,$
(84)

where here $k=2\pi n/L$.

Putting everything together yields

$c^{(1)}(E)=-\lim_{L\to\infty}\frac{1}{L}\sum_{n}\pi\frac{8k}{\omega}\bigg{(}\delta_{I=1}(\omega)\frac{\partial}{\partial\omega}+\frac{\partial\delta_{I=1}(\omega)}{\partial\omega}\bigg{)}\frac{8k}{\omega}\delta_{\epsilon}(E-\omega)\rho^{(2)}(\omega)\bigg{|}_{\omega=2\sqrt{m^{2}+k^{2}}}\,,$
(85)

where the first term in parenthesis arises from the energy shift and the
second from the matrix element. Evaluating the limit gives

$c^{(1)}(E)=-\int_{2m}^{\infty}d\omega\,\frac{d}{d\omega}\bigg{[}\delta_{I=1}(\omega)\frac{8\sqrt{\omega^{2}/4-m^{2}}}{\omega}\delta_{\epsilon}(E-\omega)\rho^{(2)}(\omega)\bigg{]}\,.$
(86)

The precise condition on $\delta_{\epsilon}(x)$ is therefore that the integral
above should vanish. Denoting the quantity in square brackets by $f(\omega)$,
this holds if $df(\omega)/d\omega$ is integrable, $f(\omega)$ is
differentiable everywhere, and $f(\omega)$ vanishes at both $\omega=2m$ and
infinity. All the smearing kernels considered in this work satisfy this
condition as do much more aggressive choices, such as $\delta_{\epsilon}(x)$
with compact support in the region $[x-\epsilon,x+\epsilon]$.

###  C.2 Exponentially suppressed volume effects

We now demonstrate that for a more restrictive class of smearing kernels,
finite-volume effects are indeed exponentially suppressed in $L$. This result
also holds for all of the smearing kernels used in this work. The argument is
based on an elegant identity that is closely related to the derivation of the
Lellouch-Lüscher formalism. One can combine the appearance of $(\partial
Q(\omega)/\partial\omega)$ with the definition of the finite-volume energies
to re-write Eq. (76) as

$\rho^{(2)}_{L,\epsilon}(E)\equiv 2\pi\int_{2m}^{\infty}\\!d\omega\
\delta_{\epsilon}(E-\omega)\ \sum_{n=0}^{\infty}\delta\big{(}2Q(\omega)-2\pi
n\big{)}\rho^{(2)}(\omega)\,,$ (87)

where the ’$\delta$’ in the sum over $n$ denotes the Dirac $\delta$-function.
The relation $\delta\big{(}f(x)\big{)}=|f^{\prime}(x)|^{-1}\delta(x)$ can be
used to show that this expression is equivalent to Eq. (76).

Next we use that $Q(E)$ is non-negative151515This holds for all $L$ since the
phase shift $\delta_{I=1}$ is positive. For systems with negative a scattering
phase shift $Q(E)>0$ is only guaranteed above some sufficiently large value of
$L$. to extend the sum over $n$ to include all integers. Then applying the
Poisson summation formula

$\sum_{n=-\infty}^{\infty}\delta(x-2\pi
n)=\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}e^{inx}\,$ (88)

gives

$\displaystyle\rho^{(2)}_{L,\epsilon}(E)$
$\displaystyle=\sum_{n=-\infty}^{\infty}\int_{2m}^{\infty}\\!d\omega\,\delta_{\epsilon}(E-\omega)\,e^{i2nQ(\omega)}\rho^{(2)}(\omega)\,,$
(89)
$\displaystyle=\sum_{n=-\infty}^{\infty}\int_{2m}^{\infty}\\!d\omega\,\delta_{\epsilon}(E-\omega)\,e^{inLk}[S(k)]^{n}\rho^{(2)}(\omega)\,,$
(90)

where in the second line we have substituted the definition of $Q(E)$ from Eq.
(73).

Taking the expressions for the spectral density and $S$-matrix (Eqs. (9) and
(74) respectively) and changing the variable of integration to
$k=\sqrt{\omega^{2}/4-m^{2}}$ now results in

$\displaystyle\rho^{(2)}_{L,\epsilon}(E)=\frac{3\pi^{3}}{8}\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}dk\,\delta_{\epsilon}\Big{(}E-2\sqrt{m^{2}+k^{2}}\Big{)}\,e^{inLk}\,\mathcal{I}_{n}(k)\,,$
(91)
$\displaystyle\mathcal{I}_{n}(k)\equiv\left.\bigg{[}\frac{\theta+2i\pi}{\theta-2i\pi}\frac{\theta-i\pi}{\theta+i\pi}\bigg{]}^{n}\frac{1}{\theta^{2}}\,\frac{\theta^{2}+\pi^{2}}{\theta^{2}+4\pi^{2}}\tanh^{4}\frac{\theta}{2}\right|_{\theta=2\sinh^{-1}(k/m)}\,.$
(92)

We have also used the invariance of the integrand under the simultaneous
replacements $n\to-n$ and $k\to-k$ ($\theta\to-\theta$) to ‘unfold’ the
integral to the whole real axis. The function $\mathcal{I}_{n}(k)$ is analytic
in the whole complex plane except for two cuts on $\pm i[m,\infty)$.

Assuming that $\delta_{\epsilon}\left(E-2\sqrt{m^{2}+k^{2}}\right)$ is
analytic in the complex strip $\text{Im}\,k<\mu$ for some $\mu\leq m$, one can
shift the integration contour to $\mathbb{R}+i\mu$. Concerning the kernels
used in this work, $\mu=m$ always holds for the Gaussian kernel ${\sf g}$ but
for the Cauchy kernels, ${\sf c0}$, ${\sf c1}$ and ${\sf c2}$, the singularity
at $(E-2\sqrt{m^{2}+k^{2}})^{2}=-\epsilon^{2}$ can be a distance less than $m$
from the real axis. In this case one has

$\mu=\text{Im}\sqrt{\frac{E^{2}-4m^{2}-\epsilon^{2}}{4}+\frac{iE\epsilon}{2}}=\frac{E\epsilon}{2\sqrt{E^{2}-4m^{2}}}+O(\epsilon^{2})\,,$
(93)

where the small $\epsilon$ expansion is included for illustration. The
following formulae are also only valid if the segments at infinity do not
contribute. We have confirmed this is the case for the kernels used in this
work.

Extending the integration range and shifting the contour in Eq. (91) as
described, we find that the difference between finite- and infinite-volume
smeared spectral densities can be written as

$\Delta\rho^{(2)}_{L,\epsilon}(E)\equiv\rho^{(2)}_{L,\epsilon}(E)-\rho^{(2)}_{\infty,\epsilon}(E)=\sum_{n=1}^{\infty}\mathcal{C}_{n}(E,L,\mu)e^{-n\mu
L}\,.$ (94)

The coefficient $\mathcal{C}_{n}$ is oscillatory with an envelope that is
either constant or falling with a power of $L$. It is defined as

$\mathcal{C}_{n}(E,L,\mu)\equiv\frac{3\pi^{3}}{8}\lim_{\eta\to\mu^{-}}\text{Re}\int_{-\infty}^{\infty}dx\,\delta_{\epsilon}(E-\omega_{\eta}(x))\,e^{inLx}\,\mathcal{I}_{n}(\theta_{\eta}(x))\,,$
(95)

with $\theta_{\eta}(x)=2\sinh^{-1}[(x+i\eta)/m]$ and
$\omega_{\eta}(x)=2\sqrt{(x+i\eta)^{2}+m^{2}}$.

Figure 10: Theoretically predicted two-particle spectral densities for three
sets of $\epsilon/m$ and $mL$ values (as indicated) with the Gaussian smearing
kernel. Each panel includes the smeared infinite-volume results (orange), as
well as the finite-volume smeared spectral density defined in Eq. (76) (blue).
The points arise from numerically evaluating the $n=1$ term of Eq. (94) and
combining with $\rho_{\infty,\epsilon}(E)$. Figure 11: Plots of the finite-
volume residue vs. $L$ for the Gaussian kernel (${\sf g}$, left) and Cauchy
kernel (${\sf c0}$, right), determined using Eq. (94). The dashed lines give
the predicted $e^{-\mu L}$ scaling for each (with a pre-factor chosen to
position the curve for ease of comparison); $\mu=m$ for the Gaussian and
$\mu=0.265m$ for the Cauchy, slightly exceeding $\epsilon/2$ as given by Eq.
(93). The Gaussian data includes only the $n=1$ term of Eq. (94) while the
Cauchy includes the sum over the $n=1,2,3$ terms, i.e. all terms that fall-off
slower than $e^{-mL}$.

We have also studied these results numerically and confirmed that the
straightforward expression based on a sum over finite-volume states in Eq.
(76) matches the sum over Poisson modes evaluated both with the original
contour (Eq. (91)) and the shifted contour which makes the $e^{-\mu L}$
scaling manifest (Eq. (94)). In Fig. 10 we show the results for the Gaussian
kernel plotted as a function of $E$ for various values of fixed $\epsilon$ and
$L$. In Fig. 11 we plot $|\Delta\rho^{(2)}_{L,\epsilon}(E)|$ at fixed $E$ and
$\epsilon$ versus $L$ for both the Gaussian and Cauchy kernels. We see that
the volume effects are oscillatory functions of $L$ with an envelope decaying
according to the predicted exponential. The derivations of this appendix
therefore suggest that finite-volume effects are small for the numerical
results in Sec. V.1.

##  Appendix D Simulation algorithm

We employ the single-cluster algorithm and associated cluster estimators from
Ref. Luscher:1990ck  which are briefly summarized here. A cluster update
proceeds as follows. First, a random vector $r\in{\mathbb{R}}^{3},\,|r|=1$ is
drawn uniformly from the unit sphere. Then a ‘seed’ site is chosen uniformly
as the first member of the cluster. For each new site $x$ added to the
cluster, consider all non-cluster sites among the four nearest neighbors of
$x$. A neighbor $y$ is added to the cluster with probability

$\displaystyle p_{\rm
add}=1-\exp\left[\min\\{-2\beta\sigma_{r}(x)\sigma_{r}(y),0\\}\right]$ (96)

where $\sigma_{r}(x)=\sigma(x)\cdot r$. After all cluster neighbors have been
considered for addition, all cluster sites are updated according to

$\displaystyle\sigma^{a}(x)\rightarrow\sigma^{a}(x)-2\sigma_{r}(x)r^{a}.$ (97)

In order to employ a cluster estimator for $C(a\tau)$, which contains four
fundamential fields, a second orthogonal cluster update is required. This
proceeds by choosing a second random vector $u$ from the unit sphere with the
constraint that $r\cdot u=0$ and then performing a second cluster update in
the same manner. The combination of a single $r$-update followed by a
$u$-update is henceforth referred to as a ‘cluster update’, the number of
which are tabulated in Tab. 2 for each ensemble.

The cluster estimator for $C(a\tau)$ is built from

$\displaystyle\tilde{\sigma}_{r,u}(x)=\sqrt{\frac{|\Lambda|}{N_{r,u}}}\theta_{r,u}(x)\sigma_{r,u}(x)$
(98)

where $N_{r,u}/|\Lambda|$ is the cluster fraction and $\theta_{r,u}(x)$ the
characteristic function. Expectation values $\langle\dots\rangle_{\rm 1C}$
denote integration over all possible pairs of orthogonal single-cluster
configurations as well as the usual integration over the field variables.
Straightforward application of Eq. (36) gives

$\displaystyle\langle
j^{a}_{\mu}(x)j^{a}_{\mu}(y)\rangle=2\beta^{2}\,P(ab|cd)\,\langle\sigma_{a}(x)\sigma_{b}(x+a\hat{\mu})\sigma_{c}(y)\sigma_{d}(y+a\hat{\mu})\rangle=12\beta^{2}\,\langle\tilde{\sigma}_{[r}(x)\tilde{\sigma}_{u]}(x+a\hat{\mu})\,\tilde{\sigma}_{[r}(y)\sigma_{u]}(y+a\hat{\mu})\rangle_{\rm
1C}$ (99)

where $P(ab|cd)=\frac{1}{2}(\delta_{ac}\delta_{bd}-\delta_{ad}\delta_{bc})$
and
$\tilde{\sigma}_{[r}(x)\tilde{\sigma}_{u]}(y)=\frac{1}{2}\\{\tilde{\sigma}_{r}(x)\tilde{\sigma}_{u}(y)-\tilde{\sigma}_{u}(x)\tilde{\sigma}_{r}(y)\\}$.
The correlator is then given by

$\displaystyle
C(a\tau)=\frac{12\beta^{2}}{L/a}\langle\Phi_{ru}(a\tau)\,\Phi_{ru}(0)\rangle_{\rm
1C},\qquad\Phi_{ru}(a\tau)=\sum_{\bm{x}}\tilde{\sigma}_{[r}(a\tau,\bm{x})\tilde{\sigma}_{u]}(a\tau,\bm{x}+a\hat{1})\,,$
(100)

where (although not denoted explicitly) time translation invariance is
employed to average over all equivalent time separations on the periodic
torus.

The estimator for $C(a\tau)$ in Eq. (100) is not positive definite, in
contrast to the cluster estimator for the two-point function of the
fundamental fields, which requires only a single cluster. The computational
cost of $\Phi_{ru}(a\tau)$ for a single update scales only weakly with the
lattice size at fixed $\beta$, since only pairs of sites which are neighbors
in the spatial direction and belong to different clusters contribute. However
this overlap becomes increasingly unlikely, resulting in an increasing
statistical error for a fixed number of cluster updates. Despite this ‘cluster
cutoff’, the variance of the estimator in Eq. (100) decays exponentially with
increasing $\tau$, while the one for standard estimator (given by the middle
expression in Eq. (99)) approaches a constant. The cluster estimator thus
results in a significant improvement in the signal-to-noise ratio, which
empirically decays with a rate roughly similar to $m$.

## References

  * (1) M. Hansen, A. Lupo, and N. Tantalo, Extraction of spectral densities from lattice correlators, Phys. Rev. D99 (2019), no. 9 094508, [arXiv:1903.06476].
  * (2) M. Lüscher, Two particle states on a torus and their relation to the scattering matrix, Nucl. Phys. B354 (1991) 531–578.
  * (3) L. Lellouch and M. Lüscher, Weak transition matrix elements from finite volume correlation functions, Commun. Math. Phys. 219 (2001) 31–44, [hep-lat/0003023].
  * (4) M. T. Hansen and S. R. Sharpe, Lattice QCD and Three-particle Decays of Resonances, Ann. Rev. Nucl. Part. Sci. 69 (2019) 65–107, [arXiv:1901.00483].
  * (5) A. Rusetsky, Three particles on the lattice, PoS LATTICE2019 (2019) 281, [arXiv:1911.01253].
  * (6) J. C. A. Barata and K. Fredenhagen, Particle scattering in Euclidean lattice field theories, Commun. Math. Phys. 138 (1991) 507–520.
  * (7) K.-F. Liu, Parton Distribution Function from the Hadronic Tensor on the Lattice, PoS LATTICE2015 (2016) 115, [arXiv:1603.07352].
  * (8) M. T. Hansen, H. B. Meyer, and D. Robaina, From deep inelastic scattering to heavy-flavor semileptonic decays: Total rates into multihadron final states from lattice QCD, Phys. Rev. D96 (2017), no. 9 094513, [arXiv:1704.08993].
  * (9) P. Gambino and S. Hashimoto, Inclusive Semileptonic Decays from Lattice QCD, Phys. Rev. Lett. 125 (2020), no. 3 032001, [arXiv:2005.13730].
  * (10) J. Bulava and M. T. Hansen, Scattering amplitudes from finite-volume spectral functions, Phys. Rev. D100 (2019), no. 3 034521, [arXiv:1903.11735].
  * (11) M. Bruno and M. T. Hansen, Variations on the Maiani-Testa approach and the inverse problem, arXiv:2012.11488.
  * (12) H. B. Meyer, Transport Properties of the Quark-Gluon Plasma: A Lattice QCD Perspective, Eur. Phys. J. A 47 (2011) 86, [arXiv:1104.3708].
  * (13) G. Backus and F. Gilbert, The resolving power of gross earth data, Geophysical Journal International 16 (1968), no. 2 169–205.
  * (14) G. Backus and F. Gilbert, Uniqueness in the inversion of inaccurate gross earth data, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 266 (1970), no. 1173 123–192, [http://rsta.royalsocietypublishing.org/content/266/1173/123.full.pdf].
  * (15) F. Pijpers and M. Thompson, Faster formulations of the optimally localized averages method for helioseismic inversions, Astronomy and Astrophysics 262 (08, 1992) L33–L36.
  * (16) G. Bailas, S. Hashimoto, and T. Ishikawa, Reconstruction of smeared spectral function from Euclidean correlation functions, PTEP 2020 (2020), no. 4 043B07, [arXiv:2001.11779].
  * (17) P. Sabatier, ed., Inverse Methods in Action.  Springer Verlag, 1990.
  * (18) R.-A. Tripolt, P. Gubler, M. Ulybyshev, and L. Von Smekal, Numerical analytic continuation of Euclidean data, Comput. Phys. Commun. 237 (2019) 129–142, [arXiv:1801.10348].
  * (19) M. Luscher and U. Wolff, How to Calculate the Elastic Scattering Matrix in Two-dimensional Quantum Field Theories by Numerical Simulation, Nucl. Phys. B339 (1990) 222–252.
  * (20) H. B. Meyer, Lattice QCD and the Timelike Pion Form Factor, Phys. Rev. Lett. 107 (2011) 072002, [arXiv:1105.1892].
  * (21) A. Francis, G. von Hippel, H. B. Meyer, and F. Jegerlehner, Vector correlator and scale determination in lattice QCD, PoS LATTICE2013 (2013) 320, [arXiv:1312.0035].
  * (22) M. Hasenbusch, P. Hasenfratz, F. Niedermayer, B. Seefeld, and U. Wolff, Nonstandard cutoff effects in the nonlinear sigma model, Nucl. Phys. B Proc. Suppl. 106 (2002) 911–913, [hep-lat/0110202].
  * (23) J. Balog, F. Niedermayer, and P. Weisz, Logarithmic corrections to O(a**2) lattice artifacts, Phys. Lett. B 676 (2009) 188–192, [arXiv:0901.4033].
  * (24) J. Balog, F. Niedermayer, and P. Weisz, The Puzzle of apparent linear lattice artifacts in the 2d non-linear sigma-model and Symanzik’s solution, Nucl. Phys. B 824 (2010) 563–615, [arXiv:0905.1730].
  * (25) E. C. Poggio, H. R. Quinn, and S. Weinberg, Smearing the Quark Model, Phys. Rev. D 13 (1976) 1958.
  * (26) A. Francis, P. Fritzsch, M. Lüscher, and A. Rago, Master-field simulations of O($a$)-improved lattice QCD: Algorithms, stability and exactness, Comput. Phys. Commun. 255 (2020) 107355, [arXiv:1911.04533].
  * (27) L. Giusti and M. Lüscher, Topological susceptibility at $T>T_{\rm c}$ from master-field simulations of the SU(3) gauge theory, Eur. Phys. J. C 79 (2019), no. 3 207, [arXiv:1812.02062].
  * (28) M. Cè, M. Bruno, J. Bulava, A. Francis, P. Fritzsch, J. R. Green, M. T. Hansen, and A. Rago, Approaching the master-field: Hadronic observables in large volumes, in 38th International Symposium on Lattice Field Theory, 10, 2021.  arXiv:2110.15375.
  * (29) M. Karowski and P. Weisz, Exact Form-Factors in (1+1)-Dimensional Field Theoretic Models with Soliton Behavior, Nucl. Phys. B 139 (1978) 455–476.
  * (30) A. B. Zamolodchikov and A. B. Zamolodchikov, Relativistic Factorized S Matrix in Two-Dimensions Having O(N) Isotopic Symmetry, Nucl. Phys. B133 (1978) 525. [JETP Lett.26,457(1977)].
  * (31) J. Balog and M. Niedermaier, Off-shell dynamics of the O(3) NLS model beyond Monte Carlo and perturbation theory, Nucl. Phys. B 500 (1997) 421–461, [hep-th/9612039].
  * (32) J. Balog, M. Niedermaier, and T. Hauer, Perturbative versus nonperturbative QFT: Lessons from the O(3) NLS model, Phys. Lett. B 386 (1996) 224–232, [hep-th/9604161].
  * (33) W. H. Press and G. R. Farrar, Recursive stratified sampling for multidimensional monte carlo integration, Computers in Physics 4 (1990), no. 2 190–195, [https://aip.scitation.org/doi/pdf/10.1063/1.4822899].
  * (34) F. Johansson, Arb: efficient arbitrary-precision midpoint-radius interval arithmetic, IEEE Transactions on Computers 66 (2017) 1281–1292.
  * (35) L. Fousse, G. Hanrot, V. Lefèvre, P. Pélissier, and P. Zimmermann, Mpfr: A multiple-precision binary floating-point library with correct rounding, ACM Trans. Math. Softw. 33 (June, 2007) 13¿es.
  * (36) M. Mori and M. Sugihara, The double-exponential transformation in numerical analysis, Journal of Computational and Applied Mathematics 127 (2001), no. 1 287–296. Numerical Analysis 2000. Vol. V: Quadrature and Orthogonal Polynomials.
  * (37) R. A. Briceño, J. V. Guerrero, M. T. Hansen, and A. M. Sturzu, Role of boundary conditions in quantum computations of scattering observables, Phys. Rev. D 103 (2021), no. 1 014506, [arXiv:2007.01155].