for the Hadron Spectrum Collaboration

# Timelike meson form-factors beyond the elastic region from lattice QCD

Felipe G. Ortega-Gama<EMAIL_ADDRESS>Department of Physics, William &
Mary, Williamsburg, VA 23187, USA  Jozef J. Dudek<EMAIL_ADDRESS>Department
of Physics, William & Mary, Williamsburg, VA 23187, USA  Thomas Jefferson
National Accelerator Facility, 12000 Jefferson Avenue, Newport News, VA 23606,
USA  Robert G. Edwards  Thomas Jefferson National Accelerator Facility, 12000
Jefferson Avenue, Newport News, VA 23606, USA

###### Abstract

We present a calculation of the vector-isovector timelike form-factors of the
pion and the kaon using lattice quantum chromodynamics. We calculate two-point
correlation functions with $m_{\pi}\sim 280$ MeV, extracting both the finite-
volume spectrum and matrix elements for these states created from the vacuum
by a vector current. After determining the coupled-channel
$\pi\pi,K\overline{K}$ scattering amplitudes, we perform the necessary
correction for the significant finite-volume effects present in the current
matrix elements, leading to the timelike form-factors. We find these to be
dominated by the presence of the $\rho$ resonance, and we extract its decay
constant by an analytic continuation of the amplitudes to the resonance pole.
In addition, the spacelike pion form-factor is determined on the same lattice
configurations, and a dispersive parameterization is used to simultaneously
describe the spacelike and elastic timelike regions.

††preprint: JLAB-THY-24-4132

##  I Introduction

The lightest octet of pseudoscalar mesons plays a unique role in QCD, being
bound states of quarks and gluons, while also serving as the pseudo-Goldstone
bosons of broken chiral symmetry. How this is manifested can be explored via a
comprehensive description of their internal structure, and this can advance
our understanding of the non-perturbative regime of QCD. Quantitative measures
of internal structure can come from transitions mediated by external probes
such as electroweak currents; and one of the simplest transitions, which is
widely studied both theoretically and experimentally, is the vector form-
factor of the pion, $f_{\pi}(P^{2})$, corresponding to the interaction vertex
between a photon of virtuality $P^{2}$ and two charged pions. Depending upon
the external kinematics, this can describe two kinematically disjoint
processes: for _spacelike_ $P^{2}<0$ it corresponds to the elastic
electromagnetic response of a $\pi^{\pm}$ meson, while for _timelike_
$P^{2}>4m_{\pi}^{2}$ it describes the production amplitude of a
$\pi^{+}\pi^{-}$ pair from the QCD vacuum. Experimentally, the spacelike form-
factor is extracted from elastic $e^{-}\pi^{\pm}$ scattering, either from a
pion beam on atomic electrons Amendolia  _et al._ (1986), or electroproduction
of charged pions from a nucleus Huber  _et al._ (2008), while the timelike
form-factor can be retrieved from the $e^{+}e^{-}\to\pi^{+}\pi^{-}$ cross-
section Fang  _et al._ (2021), or the tau-lepton decay rate,
$\tau^{-}\to\pi^{-}\pi^{0}\nu_{\tau}$ Fujikawa  _et al._ (2008). Ref.
Colangelo  _et al._ (2019) presents a recent comprehensive review. Below the
$K\overline{K}$ threshold, the largest contribution to the timelike pion form-
factor comes from the $\rho$ resonance which generates a broad peak in
$f_{\pi}(P^{2})$. Once past the $K\overline{K}$ threshold, $P^{2}>4m_{K}^{2}$,
the production process becomes inelastic, and the amplitude describing
$\gamma\to K\overline{K}$ can be expressed in terms of the kaon timelike form-
factor.

Theoretical calculation of the pion form-factor from first principles in QCD
is possible within lattice QCD, a systematically improvable numerical
implementation of the QCD path integral considered on a discrete grid of
space-time points of finite extent. Owing to the Euclidean signature used in
this approach, real-time dynamics like scattering are not directly accessible,
but the discrete spectrum of finite-volume states and local matrix elements
can be extracted from correlation functions. The form-factor in the spacelike
region can be obtained from matrix elements directly accessible from three-
point correlation functions having only exponentially suppressed finite-volume
corrections, with many previous calculations being reported (a summary with
references can be found in Ref. Aoki  _et al._ (2022)). In the timelike
region, on the other hand, the vector current produces $\pi\pi$ pairs which in
infinite volume strongly rescatter, and in finite volume these pions can go
on-shell and sample the boundary of the lattice, leading to significant
multiplicative finite-volume corrections to the matrix elements extracted from
two-point correlation functions. Evaluation of these correction factors
requires knowledge of the meson-meson scattering amplitude with vector quantum
numbers, but this can be determined applying the Lüscher approach Luscher
(1986a, b) to the discrete spectrum of states extracted from a matrix of two-
point correlation functions, a technique that is now well established even in
the case of coupled-channels Briceno  _et al._ (2018a).

The finite-volume correction factor required in production amplitudes follows
from the pioneering work of Lellouch and Lüscher Lellouch and Luscher (2001)
with the case of the timelike form-factor of the pion presented in Ref. Meyer
(2011). The formalism has been applied in explicit lattice QCD calculations,
restricted to the elastic region of $\pi\pi$ production, in Refs. Feng  _et
al._ (2015); Andersen  _et al._ (2019); Erben  _et al._ (2020). Many more
lattice calculations of $\pi\pi$ scattering have been performed, which is a
prerequisite for the form-factor volume corrections, see Ref. Mai  _et al._
(2023) for a recent review of these.

In this paper, we will perform a lattice QCD calculation extended into the
inelastic energy region where production of $K\overline{K}$ as well as
$\pi\pi$ is possible, showing that the appropriate finite-volume corrections
can be applied in this coupled-channel situation. In addition, we will compute
the spacelike pion form-factor on the same lattice configurations, and present
amplitudes with appropriate analytic properties which describe the form-factor
across both spacelike and elastic timelike regions.

##  II Meson-meson production

The amplitude describing the production of a pair of pseudoscalars from the
vacuum by a conserved vector current, $\mathcal{J}^{\mu}$, can be decomposed
into a kinematic factor multiplying a single Lorentz-invariant form factor,

$\displaystyle\mathcal{H}_{a}^{\mu}(p_{1},p_{2})$
$\displaystyle\equiv\matrixelement{p_{1},p_{2};a}{\mathcal{J}^{\mu}(x=0)}{0}$
$\displaystyle=(p_{1}-p_{2})^{\mu}\,f_{a}(s)\,,$ (1)

where $a$ indicates the pair ($\\{\pi\pi,K\overline{K},\ldots\\}$) and
$p_{i=1,2}$ the four-momentum of each pseudoscalar meson. The form factors,
$f_{a}(s)$, are expressed in terms of the square of the center-of-momentum
(COM) frame energy, ${s=(p_{1}+p_{2})^{2}}$, which corresponds to the
_timelike_ virtuality, $s>0$, of the current in the $s$-channel production.
The crossed-channel process in which a pion (or kaon) absorbs an off-shell
current is described by the same form factor for _spacelike_ virtualities,
$s<0$, as indicated in Fig. 1.

(a) Spacelike $\pi$ form factor. (b) Timelike $\pi$ form factor. (c) Timelike
$K$ form factor.

Figure 1: Diagrams representing the processes considered in this paper. The
spacelike and timelike form factors of the pion are aspects of a single
function of the photon virtuality, $f_{\pi}(P^{2})$.

While this amplitude decomposition by construction describes only a single
$J^{P}=1^{-}$ partial-wave, explicit partial-wave projection into definite
spin-projection, $m$, brings it into a form more natural for application of
the finite-volume formalism,

$\mathcal{H}^{\mu}_{a,m}(P)=\int\\!\frac{\differential\Omega^{\star}}{\sqrt{4\pi}}\,Y^{*}_{1m}(\Omega^{\star})\,\mathcal{H}^{\mu}_{a}(p_{1},p_{2})\,,$

in a frame with four-momentum $P^{\mu}=p_{1}^{\mu}+p_{2}^{\mu}$. Evaluating
the integral yields

$\mathcal{H}_{a,m}^{\mu}(P)=K^{\mu}_{m}\,k_{a}^{\star}(s)\,f_{a}(s)\,,$ (2)

with $K^{\mu}_{m}=\sqrt{\tfrac{4}{3}}\,\epsilon^{\mu*}(P,m)$, which features
the polarization vector, and, characteristic of a $P$-wave process, one power
of $k^{\star}_{a}(s)$, the magnitude of the momentum of a pseudoscalar in the
COM-frame.

As is also the case for the scattering $S$-matrix describing the hadronic
processes $\pi\pi\to\pi\pi,\pi\pi\to K\overline{K},\ldots$, the form factors
are constrained by the principles of unitarity, causality, and crossing-
symmetry. Of these, unitarity has the simplest implementation, and the most
direct consequence – it can be expressed as a constraint on the imaginary part
of the amplitude leading to

$k_{a}^{\star}(s)\,\text{Im}f_{a}(s)=\sum\nolimits_{b}\,\mathcal{M}_{ab}^{*}(s)\,\rho_{b}(s)\,k_{b}^{\star}(s)\,f_{b}(s)\,,$
(3)

where the sum runs over all kinematically open pair channels at COM-frame
energy $\sqrt{s}$. The presence of the $P$-wave hadronic scattering matrix,
$\mathcal{M}$, in this expression can be interpreted as reflecting the strong
hadron rescattering that must occur after a meson-meson pair is produced by
the current. The corresponding unitarity constraint on $\mathcal{M}$ reads,

$\text{Im}\,\mathcal{M}_{ac}=\sum\nolimits_{b}\mathcal{M}_{ab}^{*}\,\rho_{b}\,\mathcal{M}_{bc}\,,$
(4)

where the presence of the phase-space for each channel,
$\rho_{a}(s)=\tfrac{k_{a}^{\star}(s)}{8\pi\sqrt{s}}$, indicates the unitarity-
enforced analytic structure with a branch cut opening at each new threshold.
Moving down from real energy values through any one of these cuts leads to a
neighbouring Riemann sheet, on which pole singularities may appear, reflecting
the presence of _resonances_. For the case of $\pi\pi$ scattering, this is how
the $\rho$ appears mathematically.

The similarity in structure between Eqs. 3 and 4 indicates that the imaginary
parts of the functions $f_{a}(s)$ and $\mathcal{M}_{ab}(s)$ are related, and
in the case of _elastic_ scattering, this is the origin of Watson’s theorem,
which states that the _phase_ of $f(s)$ must be equal to that of
$\mathcal{M}$, i.e. the elastic scattering phase-shift. In the more general
coupled-channel case, a parameterization which trivially satisfies Eq. 3 is,

$f_{a}=\sum\nolimits_{b}\tfrac{1}{k_{a}^{\star}}\mathcal{M}_{ab}\tfrac{1}{k_{b}^{\star}}\,\mathcal{F}_{b}\,,$
(5)

where functions, $\mathcal{F}_{a}(s)$, have been introduced which do not have
the branch cut structure present in $\mathcal{M}(s)$, but rather are smooth,
real functions in the kinematically accessible scattering region. The presence
of explicit factors of $k_{a}^{\star}(s)$ reflects the $P$-wave nature of both
production and scattering processes.

Analyticity of the production and scattering amplitudes, which is related to
both unitarity and causality, is less straightforward to implement as a
constraint on explicit parameterizations. It can be introduced by means of
Cauchy’s integral formula applied paying attention to the allowed
singularities of the amplitude, such as the branch-cut discontinuities imposed
by unitarity, and any bound-state poles. The resulting integrals, known as
dispersion relations, typically relate the real part of amplitudes at any
energy to integrals over semi-infinite energy regions of their imaginary part,
modulo some known kernel function.

In the case that scattering is treated as being elastic at all energies, a
solution of the dispersion relations known as the Omnès–form exists Omnes
(1958), which suggests the following decomposition for the elastic form
factor,

$f(s)=\Omega(s)\,\mathcal{F}_{\Omega}(s)\,,$ (6)

where the Omnès factor, $\Omega(s)$, replaces the factor $\mathcal{M}/k^{\star
2}$ of the elastic version of Eq. 5. The Omnès factor has the same phase as
the scattering amplitude, $\mathcal{M}$, as dictated by unitarity, but its
magnitude is determined by a dispersion relation – an explicit form will be
presented later in this manuscript. The function $\mathcal{F}_{\Omega}(s)$ in
Eq. 6 is a smooth real function particular to the pion form-factor.

The analytic properties of the Omnès factor, which by construction features no
singularities for $s<0$, makes Eq. 6 applicable in _both_ the timelike and
spacelike regions, while to use Eq. 5 in the spacelike region would require
the extrapolation of $\mathcal{M}(s)$ into a region where crossed-channel
singularities appear,111Although in practical parameterizations of
$\mathcal{M}(s)$ we may choose not include them when considering only energy
regions away from their influence. and these would need to be cancelled by
singularities in $\mathcal{F}(s)$.

The Omnès form presented above is restricted to _elastic_ scattering, and the
extension to the coupled-channel case is not so simple, with a set of coupled
integral equations known as the Muskhelishvili–Omnès problem needing to be
solved. For the case presented in this paper, there is limited benefit to such
an undertaking, as we will see later.

##  III Finite volume formalism

Lattice QCD calculations necessarily work in a finite spatial volume, and as
such do not feature a continuum of multiparticle states, but rather a discrete
spectrum sensitive to the volume. For energies below three-particle
thresholds, the finite-volume spectrum is related to two-particle scattering
amplitudes by the Lüscher determinant condition Luscher (1986a, b, 1991a,
1991b); Rummukainen and Gottlieb (1995); He  _et al._ (2005); Kim  _et al._
(2005); Christ  _et al._ (2005); Fu (2012); Guo  _et al._ (2013); Hansen and
Sharpe (2012); Briceno and Davoudi (2013); Gockeler  _et al._ (2012); Leskovec
and Prelovsek (2012); Briceno (2014); Briceno  _et al._ (2018a),

$\det[\mathcal{M}(E)+F^{-1}(E,L)\big{]}=0\,,$ (7)

where $F$ is a matrix of known functions of essentially kinematic origin
dependent on the $L\times L\times L$ volume of the periodic lattice, and
$\mathcal{M}$ is a matrix containing scattering amplitudes, diagonal in
partial-waves, but dense in channel-space when coupled-channels are
kinematically accessible. The finite-volume spectrum,
$\\{E_{\mathfrak{n}}(L)\\}$, corresponds to the discrete set of solutions to
this equation for a given $\mathcal{M}(E)$, and as such, in general, finite-
volume eigenstates, $\ket{\mathfrak{n}}_{L}$, cannot be associated with a
particular channel or partial-wave.

We respect the cubic symmetry of the spatial lattice boundary by computing
spectra in irreducible representations, ‘ _irreps_ ’, of the reduced
rotational symmetry group, and these contain subductions of multiple values of
angular momentum, leading to the partial-wave space in the determinant
condition. More constraint on scattering amplitudes can be obtained by
computing spectra in frames in which the two-particle system moves relative to
the fixed lattice, and in this case irreps of the _little group_ that
preserves rotations of the cube around the direction of motion are used.

Parameterizing the energy dependence of partial-wave scattering amplitudes,
$\mathcal{M}(E)$, we can solve Eq. 7 for volumes and irreps in which we have
computed the lattice QCD spectrum, to find ‘model’ spectra. Free parameters in
the amplitudes can then be adjusted to bring the model spectra into agreement
with the computed spectra, thus providing hadron scattering amplitudes
constrained by QCD dynamics. This approach has been successfully applied to
several systems of coupled-channel scattering, see Refs. Dudek  _et al._
(2014); Wilson  _et al._ (2015a, b); Dudek  _et al._ (2016); Moir  _et al._
(2016); Briceno  _et al._ (2018b); Woss  _et al._ (2019, 2021); Bulava  _et
al._ (2024); Prelovsek  _et al._ (2021); Wilson  _et al._ (2024), with an
efficient approach to solving Eq. 7 in the coupled-channel case presented in
Ref. Woss  _et al._ (2020). The approach is reviewed in Ref. Briceno  _et al._
(2018a).

_Production_ amplitudes where a two-particle state is generated by the action
of a current on the vacuum (or on a single-particle state), as introduced in
the previous section, can be accessed by computing current matrix elements
featuring in the initial or final state the discrete finite-volume states
discussed above. The relationship of these finite-volume matrix elements to
the infinite volume amplitudes at the same energy is given by,

$\big{|}\matrixelement{\mathfrak{n}}{\mathcal{J}^{\mu}(x=0)}{0}_{L}\big{|}^{2}=\\\
\tfrac{1}{2E_{\mathfrak{n}}L^{3}}\,\sum\nolimits_{a,b}\mathcal{H}^{\mu}_{a}(E_{\mathfrak{n}})\,\widetilde{\mathcal{R}}_{a,b}(E_{\mathfrak{n}},L)\,\mathcal{H}^{\mu}_{b}(E_{\mathfrak{n}})\,,$
(8)

where the matrix $\widetilde{\mathcal{R}}$, the coupled-channel generalization
of the Lellouch-Lüscher factor Lellouch and Luscher (2001), is related to the
residue of the finite-volume two-particle propagator at the energy,
$E_{\mathfrak{n}}(L)$, where the propagator has a pole Briceño  _et al._
(2015); Briceño and Hansen (2015). The explicit form of the matrix
$\widetilde{\mathcal{R}}$ is,

$\widetilde{\mathcal{R}}(E_{\mathfrak{n}},L)=2E_{\mathfrak{n}}\cdot\lim_{E\to
E_{\mathfrak{n}}}\frac{E-E_{\mathfrak{n}}}{\mathcal{M}(E)+F^{-1}(E,L)}\,,$

where the matrix in the denominator can only have a single vanishing
eigenvalue at $E_{\mathfrak{n}}$ in order to generate _simple poles_ in
correlation functions, consistent with causality. This means that
$\widetilde{\mathcal{R}}$ is a rank-one matrix, which can be decomposed in
terms of an outer product featuring the eigenvector associated with the
vanishing eigenvalue, although in practice it is more convenient Briceño  _et
al._ (2021) to use a matrix $\mathcal{M}\widetilde{\mathcal{R}}\mathcal{M}$,

$\displaystyle\mathcal{M}\,\widetilde{\mathcal{R}}\,\mathcal{M}$
$\displaystyle=-2E_{\mathfrak{n}}\cdot\lim_{E\to
E_{\mathfrak{n}}}\frac{E-E_{\mathfrak{n}}}{\mathcal{M}^{-1}(E)+F(E,L)}\,,$
$\displaystyle=-\frac{2E^{\star}_{\mathfrak{n}}}{\mu^{\prime\star}_{0}}\,\mathbf{w}^{\phantom{\intercal}}_{0}\,\mathbf{w}_{0}^{\intercal}\,,$
(9)

where $\mu_{0}^{\prime\star}$ is the slope in energy of the vanishing
eigenvalue of $\mathcal{M}^{-1}+F$ at $E_{\mathfrak{n}}$, which is seen to
control the normalization, while the eigenvector, $\mathbf{w}_{0}$,
distributes strength across the channel (and/or partial-wave) space.

Once $\widetilde{\mathcal{R}}$ is factorized in this way, Eq. 8 becomes a
linear relationship between the finite-volume and infinite volume amplitudes,
and the introduction of explicit factors of $\mathcal{M}$ in Eq. III removes
the strong rescattering factor in Eq. 5 leaving,

$\displaystyle\mathcal{F}_{\mathfrak{n}}^{(L)}$
$\displaystyle=\sum\nolimits_{a}\sqrt{\tfrac{2E^{\star}_{\mathfrak{n}}}{-\mu^{\prime\star}_{0}}}\,(\mathbf{w}_{0})_{a}\tfrac{1}{k_{a}^{\star}}\,\mathcal{F}_{a}\quantity(s=E^{\star
2}_{\mathfrak{n}})\,,$
$\displaystyle=\sum\nolimits_{a}\tilde{r}_{\mathfrak{n},a}(L)\,\mathcal{F}_{a}\quantity(s=E^{\star
2}_{\mathfrak{n}})\,,$ (10)

where the ‘finite volume form factor’, $\mathcal{F}_{\mathfrak{n}}^{(L)}$, is
defined in terms of the matrix-element on the left-hand-side of Eq. 8 as,

$\matrixelement{\mathfrak{n}}{\mathcal{J}^{\mu}(x=0)}{0}_{L}=\tfrac{1}{\sqrt{2E_{\mathfrak{n}}L^{3}}}\,K^{\mu}\,\mathcal{F}_{\mathfrak{n}}^{(L)}\,,$
(11)

where $K^{\mu}$ is given by the straightforward subduction of $K^{\mu}_{m}$ in
Eq. 2 into the irrep of the finite-volume state $\ket{\mathfrak{n}}_{L}$, as
described in Appendix D.

Equation 10 provides the relationship between the computed finite-volume form
factors, of which there is one per finite volume energy level,
$\mathcal{F}_{\mathfrak{n}}^{(L)}$, and the infinite-volume form factors,
$\mathcal{F}_{a}$, for each kinematically open channel $a$, evaluated at
$s=E^{\star 2}_{\mathfrak{n}}$. The vector of correction factors,
$\tilde{r}_{\mathfrak{n},a}(L)$, is available to us provided we have
determined the relevant scattering amplitudes, $\mathcal{M}(s)$, from the
finite-volume spectrum.

Clearly, solving Eq. 10 energy-by-energy for the $\mathcal{F}_{a}$ in a
coupled-channel situation is not possible, but if the form factors have their
energy dependence parameterized, then by considering multiple energy-levels, a
system of equations can be set up and a best description of multiple
$\mathcal{F}_{\mathfrak{n}}^{(L)}$ values found, yielding the energy
dependence of the form factors.

##  IV Lattice technology and correlator computation

The calculation to be presented in this paper was performed on an
$(L/a_{s})^{3}\times(T/a_{t})=24^{3}\times 256$ lattice with two degenerate
light $u,d$ quarks and a heavier $s$ quark. The light quark mass is set to be
heavier than the physical value, yielding pions of mass 284 MeV and kaons of
mass 519 MeV. Anisotropy introduced in the lattice action gives rise to a
temporal spacing $\xi=3.455(6)$ times finer than the spatial spacing.
Comparing to the physical $\Omega$ baryon mass yields $a_{t}^{-1}\approx
5.988(17)~{}\mathrm{GeV}$. In the temporal lattice units used throughout this
paper, the pion and kaon masses are $a_{t}\,m=0.04735(22),\,0.08659(14)$
respectively. This lattice has been used in previous computations of $\pi\pi$
scattering Rodas  _et al._ (2023a, b), $\pi K$ scattering Wilson  _et al._
(2019) and the electromagnetic transition $\gamma K\to K\pi$ Radhakrishnan
_et al._ (2022).

The finite-volume spectrum needed to constrain scattering amplitudes is
extracted from matrices of two-point correlation functions computed using a
basis of operators described below which make use of _smeared_ quark fields,
and we opt to smear with the _distillation_ Peardon  _et al._ (2009) approach,
using the lowest 162 eigenvectors of the gauge covariant spatial Laplacian. In
this way two-point functions factorize into pieces dependent on the particular
operator constructions, and pieces which encode the propagation of the quarks,
known as perambulators.

In order to access _production_ amplitudes for $\gamma\to\pi\pi$ and
$\gamma\to K\overline{K}$, two-point correlation functions featuring the
_local_ vector current operator $\mathcal{J}^{\mu}(x)$ on the sink timeslice
are required, and the quark fields appearing in the vector current should not
be smeared. This requires the computation of _generalized perambulators_ which
feature an insertion of the unsmeared vector current. Those used in the
present calculation are a re-used subset of those computed for use in Ref.
Radhakrishnan  _et al._ (2022). Only light-quark generalized perambulators are
required, since only the isovector component of the electomagnetic current can
feature in this case. The correlation functions describing production were
computed on an ensemble of 348 configurations.222These 348 configurations are
a subset of an ensemble of 400 configurations used to compute the spectrum.

An established procedure to obtain the discrete spectrum of finite-volume
eigenstates in a lattice calculation proceeds by diagonalizing the matrix
whose elements are ${C_{ij}(t)=\langle
0|\mathcal{O}^{\phantom{{\dagger}}}_{i}(t)\,\mathcal{O}^{\dagger}_{j}(0)|0\rangle}$,
with $\\{\mathcal{O}_{i}\\}_{i=1\ldots N}$ being a basis of operators having
the quantum numbers of the states of interest. In practice the diagonalization
is best achieved through solution of a generalized eigenvalue problem,

$\mathbf{C}(t)\,v_{\mathfrak{n}}=\lambda_{\mathfrak{n}}(t,t_{0})\,\mathbf{C}(t_{0})\,v_{\mathfrak{n}}\,,$
(12)

where the time-dependence of the eigenvalues can be fitted to give the
energies $E_{\mathfrak{n}}$, and where the particular linear combination of
basis operators given by the eigenvector weights,
$\Omega^{\dagger}_{\mathfrak{n}}=\sqrt{2E_{\mathfrak{n}}}e^{-E_{\mathfrak{n}}t_{0}/2}\sum\nolimits_{i}(v_{\mathfrak{n}})_{i}\,\mathcal{O}^{\dagger}_{i}$,
optimally interpolates the discrete state $|\mathfrak{n}\rangle$, and is
approximately orthogonal to other states in the spectrum in the sense that
$\matrixelement{\mathfrak{n}}{\Omega^{\dagger}_{\mathfrak{m}}(0)}{0}=\sqrt{2E_{\mathfrak{n}}}\,\delta_{\mathfrak{n},\mathfrak{m}}$,
up to small corrections.

The specific operator basis used in this calculation includes fermion
bilinears which resemble $q\bar{q}$ constructions and meson-meson-like
operators resembling $\pi\pi$ and $K\overline{K}$ pairs with relative
momentum. The set of allowed relative momenta is guided by the energy such
pairs would have given no dynamical interactions between them, with inclusion
of all operators producing a state below the $\phi\pi$ threshold at
$a_{t}E^{\star}\sim 0.22$. The opening of the $\omega\pi$ and $\phi\pi$
channels in the energy region of interest led us to also add operators
resembling these pairs, but it was found that upon their addition there was
negligible change in the energy of states with large overlap onto
$\pi\pi$-like and $K\overline{K}$-like operators, and as such we will present
results excluding $\omega\pi,\phi\pi$ operators from the basis.333See Appendix
A for a presentation of this, and Refs. Woss  _et al._ (2018, 2019) for some
discussion of the construction of pseudoscalar-vector meson pair operators.

The optimized operators found by solving the generalized eigenvalue problem
described above are used in the subsequent computation of correlation
functions featuring the vector current, $\langle
0|\mathcal{J}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle$, and these
correlation functions are expected to be dominated by the contribution of the
specific finite-volume state $|\mathfrak{n}\rangle$, with other states,
including those of lower energy, being suppressed Shultz  _et al._ (2015). The
extent to which this is true in practical calculation will be explored in
Section VI. In order to compute these correlation functions featuring the
unsmeared vector current, diagrams of the type shown in Figure 2 must be
computed, featuring the generalized perambulators with current insertion at
time $t$.

Figure 2: Wick contraction diagrams required to compute $\langle
0|\mathcal{J}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle$, where single
red and green lines represent light-quark and strange quark perambulators
respectively, and where the double line represents the light-quark generalized
perambulator housing the local current insertion.

##  V Energy spectrum and partial wave scattering amplitudes

Figure 3: Finite volume spectra obtained for all irreps considered in this
analysis. Dashed lines show meson-meson thresholds, while solid lines show
non-interacting meson-meson energies as a function of $L$, color-coded
according to the meson pair (dot-dashed indicating non-interacting energies
where the associated operator was not included in the variational basis). Grey
dotted lines indicate thresholds for higher-multiplicity scattering
($\pi\pi\pi\pi$, $\pi\pi\eta$, $K\overline{K}\pi$).

Matrices of correlation functions using the basis of operators described in
the previous section were computed for several irreps at rest and in moving
frames. Variational analysis of these using Eq. 12 leads to spectra that
densely span the elastic scattering region and the lower part of the coupled-
channel region.444See also Ref. Rodas  _et al._ (2023b) for a similar
extraction on the same lattices, focusing on the elastic region. Together they
provide a high degree of constraint on the $P$-wave partial-wave scattering
amplitudes over the region of interest, and the corresponding optimized
operators will be used to determine the production amplitudes at each of the
dense set of discrete energies.

Figure 3 presents the extracted spectra, with the levels below $K\overline{K}$
threshold showing clear departures from non-interacting $\pi\pi$ energies, as
well as a counting indicating an ‘additional’ level. An avoided-level-crossing
is apparent whenever a non-interacting level lies near ${a_{t}E^{\star}\sim
0.13}$. These observations strongly suggest the presence of a narrow resonance
as confirmed below through determination of scattering amplitudes. In the
coupled-channel region above $K\overline{K}$ threshold, the extracted energy
levels are seen to lie very close to non-interacting values, suggesting that
here the interactions are weak and that the scattering system is free of
additional resonances below ${a_{t}E^{\star}\sim 0.23}$.

###  V.1 Elastic $\pi\pi$ scattering

In the $\pi\pi$ elastic region below $K\overline{K}$ threshold, if the effect
of $J=3$ and higher spin amplitudes can be neglected, the Lüscher quantization
condition of Eq. (7) becomes an algebraic relation for the elastic $\pi\pi$
scattering amplitude, with each discrete energy level yielding a value of the
$J^{P}=1^{-}$ phase-shift at that energy. Alternatively, all energy levels can
be described together by parameterizing the energy-dependence of the elastic
scattering amplitude, and minimizing a $\chi^{2}$ as explained in Ref. Wilson
_et al._ (2015a). These two approaches are presented in Figure 4, where four
parameterization choices are shown, each of which prove capable of describing
17 energy levels below $a_{t}E^{\star}=0.173$ with a
$\chi^{2}/N_{\text{dof}}\approx 1.2$ (see Appendix B). As can be seen the
result is not sensitive to the details of the parameterization, and in the
remainder of this manuscript we will focus on a _reference elastic_ amplitude
constructed from a $K$-matrix built from a pole and a constant, and a
dispersively improved ‘Chew-Mandelstam’ phase-space555See e.g. Appendix B of
Ref. LABEL:Wilson:2014cna subtracted at $s=m^{2}$, the $K$-matrix pole
location,

$\displaystyle\mathcal{M}(s)$
$\displaystyle=\frac{16\pi}{\tfrac{1}{(2k^{\star})^{2}}\mathcal{K}^{-1}(s)-I_{\text{CM}}(s)}\,,$
$\displaystyle\mathcal{K}(s)$ $\displaystyle=\frac{g^{2}}{m^{2}-s}+\gamma\,,$
(13)

where the best fit to the lattice energy levels is with parameter values,

$m=$ | $0.1335\,(5)\cdot a_{t}^{-1}$ | $\begin{bmatrix}[r]1&-0.3&-0.3\\\ &1&0.7\\\ &&1\end{bmatrix}$
---|---|---
$g=$ | $0.445\,(10)$
$\gamma=$ | $(3.4\pm 2.2)\cdot a_{t}^{2}$
$\chi^{2}/N_{\text{dof}}=\frac{17.0}{17-3}=1.21$ .
(14)

The behavior observed in Fig. 4 is clearly that of a narrow resonance, and
indeed the _reference elastic_ amplitude is found to feature a pole on the
unphysical Riemann sheet at

$a_{t}\sqrt{s_{R}}=0.1328(5)-\tfrac{i}{2}0.0096(5)\,,$

with a coupling to $\pi\pi$ defined at the pole, ${\mathcal{M}(s\sim
s_{R})\sim 16\pi\,\frac{c_{\pi\pi}^{2}}{s_{R}-s}}$, of value,

$a_{t}\,c_{\pi\pi}=0.0426(11)\,e^{-i\pi\cdot 0.047(3)}\,.$

Figure 4: Elastic $\pi\pi$ $P$-wave phase-shift constrained by energy levels
below $K\overline{K}$ threshold in Fig. 3. Bands shows parameterizations
including the _reference elastic_ parameterization of Eq.14.

###  V.2 Coupled $\pi\pi$/$K\overline{K}$ scattering

Above $K\overline{K}$ threshold, each energy level is in principle sensitive
to all elements of the coupled-channel scattering matrix, and as such we must
proceed using parameterizations of this matrix. The structure of Eq. 7 is such
that only parameterizations consistent with unitarity will generate solutions.
Use of $K$-matrix forms ensures unitarity, and parameterizations in which the
$K$-matrix is constructed from poles plus polynomials in $s$ have been
successful in past attempts to describe lattice spectra, and in particular in
Ref. Wilson  _et al._ (2015b) where the $\pi\pi,K\overline{K}$ $P$-wave system
was considered at a lighter pion mass.

A simple form, which we will refer to as the _reference coupled-channel_
amplitude, capable of describing the entire spectrum, is given by,

$\displaystyle\big{[}\mathcal{M}^{-1}\big{]}_{ab}(s)$
$\displaystyle=\tfrac{1}{16\pi}\left(\tfrac{1}{2k_{a}^{\star}}\big{[}\mathcal{K}^{-1}(s)\big{]}_{ab}\tfrac{1}{2k_{b}^{\star}}-\delta_{ab}\,I_{\text{CM},a}(s)\right)\,,$
$\displaystyle\mathcal{K}_{ab}(s)$
$\displaystyle=\frac{g_{a}\,g_{b}}{m^{2}-s}+\gamma_{ab}\,,$ (15)

(where a Chew-Mandelstam phase-space is used) with parameter values,

$m=$ | $0.1338\,(5)\cdot a_{t}^{-1}$ | $\begin{bmatrix}[r]1&-0.2&0.0&-0.2&0.2&-0.1\\\ &1&-0.4&0.6&-0.4&-0.4\\\ &&1&-0.2&0.8&0.9\\\ &&&1&-0.2&-0.1\\\ &&&&1&0.8\\\ &&&&&1\\\ \end{bmatrix}$
---|---|---
$g_{\pi\pi}=$ | $0.441\,(9)$
$g_{K\overline{K}}=$ | $0.17\,(30)$
$\gamma_{\pi\pi,\pi\pi}=$ | $(2.9\pm 0.9)\cdot a_{t}^{2}$
$\gamma_{\pi\pi,K\overline{K}}=$ | $-(2.4\pm 5.0)\cdot a_{t}^{2}$
$\gamma_{K\overline{K},K\overline{K}}=$ | $-(2.2\pm 4.0)\cdot a_{t}^{2}$
$\chi^{2}/N_{\text{dof}}=\frac{28.7}{32-6}=1.10$ .
(16)

This amplitude is presented in Fig. 5 in terms of the $\pi\pi$ and
$K\overline{K}$ phase-shifts and the inelasticity, whose departure from value
1 indicates channel coupling666The definition of these quantities in terms of
the scattering amplitude is given by Eq. (8) of Ref. Wilson  _et al._
(2015b).. As was suggested by the spectra shown in Figure 3, the amplitudes
indicate no strong scattering above the $K\overline{K}$ threshold, and only
mild evidence for $\pi\pi,K\overline{K}$ channel coupling. The $\rho$
resonance is clearly present in the $\pi\pi$ phase-shift, and because we use a
coupled-channel parameterization at all energies, in principle the $\rho$–pole
will have a coupling to the $K\overline{K}$ channel. Such a resonance coupling
can be rigorously defined by factorizing the amplitude at the location of the
complex resonance pole on the unphysical Riemann sheet,

$\mathcal{M}_{ab}(s\sim s_{R})\sim 16\pi\,\frac{c_{a}\,c_{b}}{s_{R}-s}\,,$
(17)

however, because the $\rho$ lies well below the $K\overline{K}$ threshold, the
coupling to $K\overline{K}$ will not be very well constrained as the impact of
the $K\overline{K}$ components of $\mathcal{M}$ is exponentially suppressed in
Eq. 7 below $K\overline{K}$ threshold. The amplitude presented above has a
pole777On the sheet closest to physical scattering below $K\overline{K}$
threshold where the imaginary parts of the $\pi\pi$ and $K\overline{K}$
momenta are negative and positive, respectively – sheet II in the usual
nomenclature. at

$a_{t}\sqrt{s_{R}}=0.1331(4)-\tfrac{i}{2}0.0095(4)\,,$ (18)

with channel couplings,

$\displaystyle a_{t}\,c_{\pi\pi}$ $\displaystyle=0.0424(8)\,e^{-i\pi\cdot
0.047(2)}\,,$ $\displaystyle a_{t}\,c_{K\overline{K}}$
$\displaystyle=0.019(33)\,e^{i\pi\cdot 0.47(5)}\,,$ (19)

which are observed to be close to real valued for the kinematically open
$\pi\pi$ channel, and close to imaginary valued for the kinematically closed
$K\overline{K}$ channel.

The pole location in this coupled-channel analysis is seen to be compatible
with the one found in the previous elastic analysis, as expected for a
resonance lying well below the $K\overline{K}$ threshold. Further, also as
expected, the $K\overline{K}$ coupling is not precisely determined. Owing to
the $P$-wave nature of the resonance, the (dimensionless) values of the
couplings with the barrier factor divided out are also relevant,

$\displaystyle\sqrt{16\pi}\,\absolutevalue{\hat{c}_{\pi\pi}}$
$\displaystyle\equiv\sqrt{16\pi}\left|\frac{c_{\pi\pi}}{k^{\star}_{\pi\pi}(s_{R})}\right|=6.41\pm
0.13\,,$ $\displaystyle\sqrt{16\pi}\,\absolutevalue{\hat{c}_{K\overline{K}}}$
$\displaystyle\equiv\sqrt{16\pi}\left|\frac{c_{K\overline{K}}}{k^{\star}_{K\overline{K}}(s_{R})}\right|=2.4\pm
4.0\,.$

As well as this _reference coupled-channel_ amplitude, parameterization
variations listed in Appendix B also prove capable of describing the spectra,
generating plots which closely resemble Fig. 5 and which have a $\rho$ pole
and couplings broadly compatible with those given above.

Figure 5:  Phase-shifts, $\delta_{1}^{a}$, and inelasticity, $\eta$, for
coupled-channel $\pi\pi,K\overline{K}$ scattering corresponding to the
amplitude given in Eq. 16. The constraining energy level locations are shown
between the panels, with levels lying close to $K\overline{K}$ non-interacting
energies colored in red. White circles show the thresholds for $\pi\pi\pi\pi$,
$\omega\pi$ and $\phi\pi$ production, in that order.

##  VI Lattice matrix elements

To obtain _production_ amplitudes, the required finite-volume matrix elements
are extracted from two-point correlation functions, $\langle
0|\mathcal{J}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle$, which feature
the electromagnetic current. This current can be decomposed into components of
definite isospin,

$\mathcal{J}=Z_{V}^{l}\tfrac{1}{\sqrt{2}}\left(\mathcal{J}_{\rho,{\text{lat}}}+\tfrac{1}{3}\mathcal{J}_{\omega_{l},\text{lat}}\right)+Z_{V}^{s}\left(-\tfrac{1}{3}\mathcal{J}_{\omega_{s},\mathrm{lat}}\right)\,,$
(20)

where the isospin-basis currents are

$\mathcal{J}_{\rho}\equiv\tfrac{1}{\sqrt{2}}\big{(}\bar{u}\Gamma
u-\bar{d}\Gamma
d\big{)},\mathcal{J}_{\omega_{l}}\equiv\tfrac{1}{\sqrt{2}}\big{(}\bar{u}\Gamma
u+\bar{d}\Gamma d\big{)},\mathcal{J}_{\omega_{s}}\equiv\bar{s}\Gamma s\,,$

and the spatially-directed vector current whose improvement at
$\mathcal{O}(a)$ is consistent with the anisotropic Clover quark action is
Shultz  _et al._ (2015),

$\bar{q}\Gamma
q=\bar{q}\gamma^{k}q+\tfrac{1}{4}(1-\xi)\,a_{t}\partial_{4}\big{(}\bar{q}\sigma^{4k}q\big{)}\,.$
(21)

In the case where production of $\pi\pi$ is considered, only the isovector
component $\mathcal{J}_{\rho}$ appears, and the multiplicative renormalization
factor, $Z_{V}^{l}$, is determined non-perturbatively using the pion form
factor at zero virtuality extracted from three-point correlation functions,
$\langle 0|\Omega_{\pi}(\Delta
t)\,\mathcal{J}(t)\,\Omega^{\dagger}_{\pi}(0)|0\rangle$, as described in Ref.
Radhakrishnan  _et al._ (2022), and Appendix H.

We compute two-point current correlation functions, $\langle
0|\mathcal{J}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle$, for each
energy level in Figure 3 using the corresponding optimized operator determined
by solving the relevant generalized eigenvalue problem. The current operator
is subduced into the irrep $\Lambda$ and projected to spatial momentum
$\mathbf{P}$ of the state $\mathfrak{n}$, and the correlation function is
computed for all $0\leq t/a_{t}<32$. The use of an optimized operator for
state $\mathfrak{n}$ should ensure a leading time-dependence of

$\langle
0|\mathcal{J}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle=e^{-E_{\mathfrak{n}}t}\,\sqrt{2E_{\mathfrak{n}}}\,L^{3}\matrixelement{0}{\mathcal{J}(x=0)}{\mathfrak{n}}_{L}+\dots\,,$

where the subleading time-dependence is expected to arise from overlap with
high-lying states. Hence, fitting $e^{E_{\mathfrak{n}}t}\cdot\langle
0|\mathcal{J}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle$ to a constant
(at late times) or a constant plus an exponential modelling excited states
(when including earlier times) would yield the desired matrix element in the
fitted constant.

In practice, timeslice-to-timeslice data correlation is reduced if the
following ratio is instead considered:

$\frac{\langle
0|\mathcal{J}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle}{\langle
0|\Omega_{\mathfrak{n}}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle}=\frac{L^{3}}{\sqrt{2E_{\mathfrak{n}}}}\matrixelement{0}{\mathcal{J}(x=0)}{\mathfrak{n}}_{L}+\dots\,,$

which can be fitted to a constant or a constant plus an exponential. For
convenience, we define the dimensionless ratio,

$R_{\mathfrak{n}}(t)\equiv\frac{\sqrt{2E_{\mathfrak{n}}}}{K(\Lambda)}\frac{\langle
0|\mathcal{J}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle}{\langle
0|\Omega_{\mathfrak{n}}(t)\,\Omega^{\dagger}_{\mathfrak{n}}(0)|0\rangle}=\sqrt{\tfrac{L^{3}}{2E_{\mathfrak{n}}}}\mathcal{F}_{\mathfrak{n}}^{(L)}+\dots\,,$
(22)

where we have also divided out the kinematic factor $K(\Lambda)$, which
corresponds to the factor $K^{\mu}$ of Eq. (11) once the current operator has
been subduced into irrep $\Lambda$ (see Tab. 4 in App. D).

Figure 6: Time dependence of $R_{\mathfrak{n}}(t)$ defined in Eq. 22 for the
three lowest discrete energy levels in the ${[111]\,A_{1}}$ irrep.
Squares/triangles show data for the unimproved/improvement term currents
(first and second terms in Eq. (21)). Curves show the timeslice description
having largest AIC weight, with the constant fit value and the
$\chi^{2}/N_{\text{dof}}$ of the fit also shown. Variations in the constant
fit value for the unimproved current over different timeslice fit windows are
shown in the right column, together with the AIC weight (in orange), and the
“model average” (in red).

The timeslice behavior of $R_{\mathfrak{n}}(t)$ for the three lowest lying
states of the [111] $A_{1}$ irrep is shown in Fig. 6. The observed relaxation
to a constant value for ${t\gtrsim t_{0}=10\,a_{t}}$, even for states above
the ground state, is a result of making use of optimized operators. We show
separately the contributions of the two terms in Eq. (21), where the
$\mathcal{O}(a)$ improvement term is observed to impact at the level of
$5-10\%$ of the leading term. For the anisotropic lattice considered in this
work, the improvement enters proportional to the energy difference between the
initial and final states of the matrix element, which is larger in magnitude
in the current case of production than in the case of three-point functions
used to extract the spacelike pion form factor, see App. H.

The $R_{\mathfrak{n}}(t)$ timeslice data can be described using constant or
constant plus exponential forms fitted over various windows. We show the
result of using a version Jay and Neil (2021) of the Akaike Information
Criterion (AIC) which assigns a weight to each fit depending on its correlated
$\chi^{2}$ and the number of degrees of freedom, allowing us to average over a
range of fits. In practice we observe that the ‘model average’ is dominated by
only four or five fits which have large weight and the averaged error is
dominated by the statistical error on these. The model average in each case is
used for all subsequent analysis in this paper.

Figure 7: Finite volume form factor in the $\pi\pi$ elastic region and above
the $K\overline{K}$ threshold. The fractional statistical uncertainties for
the form factor values are shown in the bottom panel.

Levels found lying near $K\overline{K}$ non-interacting energies (for example
the fourth state in $[111]\,A_{1}$) demand a slightly different timeslice
analysis. The heavily suppressed value of $\mathcal{F}^{(L)}_{\mathfrak{n}}$
for these states relative to the large $\mathcal{F}^{(L)}_{\mathfrak{m}}$ for
lower-lying states near the resonance energy, is such that even though the
optimized operator $\Omega_{\mathfrak{n}}$ has only a tiny overlap onto lower-
lying states $|\mathfrak{m}\rangle$, the combination
$\mathcal{F}_{\mathfrak{m}}^{(L)}\cdot\matrixelement{\mathfrak{m}}{\Omega_{\mathfrak{n}}^{\dagger}(0)}{0}$
can be significant. This gives rise to a ‘subleading’ contribution to
$R_{\mathfrak{n}}(t)$ with exponentially _growing_ time-dependence,

$R_{\mathfrak{n}}(t)=\dots+\sum_{\mathfrak{m}<\mathfrak{n}}\varepsilon_{\mathfrak{m},\mathfrak{n}}\,e^{(E_{\mathfrak{n}}-E_{\mathfrak{m}})t}+\dots\,,$

where
$\varepsilon_{\mathfrak{m},\mathfrak{n}}\propto\mathcal{F}_{\mathfrak{m}}^{(L)}\cdot\matrixelement{\mathfrak{m}}{\Omega_{\mathfrak{n}}^{\dagger}(0)}{0}$.
Approaches to handle these cases are presented in Appendix F.

The resulting values of $\mathcal{F}_{\mathfrak{n}}^{(L)}$ from timeslice fits
for energy levels across all irreps considered are summarized in Fig. 7. The
data is visibly enhanced near the resonance energy ${a_{t}E^{\star}\sim
0.135}$, while for higher energies it does not appear to have any simple
energy dependence, with a strong dependence on the irrep of state
$\mathfrak{n}$ being observed even for very similar values of
$a_{t}E^{\star}_{\mathfrak{n}}$. This is as expected given the need for
finite-volume corrections, which will be addressed in the next section.

##  VII Finite volume correction factor

The finite-volume correction factors, $\tilde{r}_{\mathfrak{n},a}(L)$, defined
in Eq. 10, are obtained making use of parameterized scattering matrices,
$\mathcal{M}(s)$, that describe the finite-volume spectra as described in
Section V. The matrix $\mathcal{M}^{-1}+F$ is eigen-decomposed at the energy
where one of its eigenvalues crosses zero, which approximates within the
scattering ‘model’ the lattice QCD computed energy, $E_{\mathfrak{n}}(L)$. At
this point the eigenvector, $\mathbf{w}_{0}$, corresponding to the zero
eigenvalue is obtained, and this is used in the computation of the required
eigenvalue slope in a finite-difference,

$\displaystyle\mu^{\prime\star}_{0}=\frac{1}{2\,\Delta
E}\,\mathbf{w}_{0}^{\intercal}\cdot\bigg{(}\big{[}$
$\displaystyle\mathcal{M}^{-1}+F\big{]}_{E_{\mathfrak{n}}+\Delta E}$
$\displaystyle-\big{[}\mathcal{M}^{-1}+F\big{]}_{E_{\mathfrak{n}}-\Delta
E}\bigg{)}\cdot\mathbf{w}_{0}\,,$ (23)

which enters into

$\tilde{r}_{\mathfrak{n},a}(L)=\sqrt{\frac{2E^{\star}_{\mathfrak{n}}}{-\mu^{\prime\star}_{0}}}\frac{\mathbf{w}_{0,a}}{k_{a}^{\star}}\,.$
(24)

Uncertainties are propagated into this quantity using jackknife via the
ensemble of scattering matrix parameter values.

As shown in Appendix B of Ref. Radhakrishnan  _et al._ (2022), in the case
that $\mathcal{M}(s)$ houses a _narrow_ resonance, for energies near to the
resonance mass, the finite-volume correction factor becomes volume-independent
and has elements proportional to the coupling of the resonance to each
channel, $a$,

$\tilde{r}_{\mathfrak{n},a}(L)\approx\sqrt{16\pi}\;\frac{c_{a}}{k^{\star}_{a}}+\mathcal{O}\left(\tfrac{\Gamma_{R}}{m_{R}}\right)\,.$

This reflects the dominance in the scattering at these energies of a spatially
localized state whose wavefunction does not sample the boundary of the finite-
volume.

At energies where $\mathcal{M}$ describes _weak_ scattering and the solutions
to Eq. 7 lie close to non-interacting meson-pair energies, the value of
$\tilde{r}_{\mathfrak{n},a}(L)$ is set largely by properties of the geometric
matrix $F(E,L)$ related to how the meson-pair relative momentum directions
subduce into the irrep under consideration. Between these two extremes the
finite-volume correction factor is sensitive to both $\mathcal{M}$ and
$F(E,L)$.

###  VII.1 Elastic finite volume correction

Below $K\overline{K}$ threshold, where only elastic $\pi\pi$ scattering is
relevant, $\tilde{r}_{\mathfrak{n}}$ can be considered to be just a single
number correcting the normalization of the matrix element for each finite-
volume energy level. In Fig. 8 we present values of $\tilde{r}_{\mathfrak{n}}$
for the 17 energy levels below $K\overline{K}$ threshold, computed using the
_reference elastic_ amplitude.

The very similar values observed for four energy levels in different irreps
close to $a_{t}E^{\star}=0.135$ are explained by the presence of the narrow
$\rho$ resonance, with the values of $\tilde{r}_{\mathfrak{n}}$ being seen to
be almost equal to the value of $\sqrt{16\pi}|\hat{c}_{\pi\pi}|$ for this
amplitude.

Differing values of $\tilde{r}_{\mathfrak{n}}$ are observed for the two
highest levels presented, each with $\mathbf{P}=[111]$, lying close to the
$\pi_{110}\pi_{100}$ non-interacting energy at $a_{t}E^{\star}_{\text{non-
int}}\approx 0.16$. These values reflect the difference between a nearly non-
interacting helicity 0 $\pi\pi$ state subduced into the $A_{1}$ irrep, and
helicity $\pm 1$ states subduced into $E_{2}$, and the same pattern of
magnitudes can be observed in the lattice QCD computed matrix elements for the
corresponding states in Fig. 7.

Figure 8: Elastic $\pi\pi$ finite-volume factors, $\tilde{r}_{\mathfrak{n}}$,
computed using the _reference elastic_ scattering amplitude parameterization.

###  VII.2 Coupled channel finite volume correction

In the coupled-channel case, the object in parentheses in Eq. VII is a matrix,
with the projection into the zero-crossing eigenvalue achieved with the vector
$\mathbf{w}_{0}$. Above the $K\overline{K}$ threshold, the matrix is real and
symmetric, and the components of $\mathbf{w}_{0}$ are real. Above $\pi\pi$
threshold, but below $K\overline{K}$ threshold this is not the case, and the
presence of $\ell=1$ angular momentum barrier factors in $\mathcal{M}$ and $F$
causes $(\mathbf{w}_{0})_{a}\propto k_{a}^{\star}$ such that kinematically
closed channels become imaginary components of the vector. In the computation
of $\tilde{r}_{\mathfrak{n}}$, this phase is cancelled by the explicit factor
of $1/k_{a}^{\star}$ in Eq. 24, yielding a real valued correction vector (see
also the discussion in Appendix A of Ref. Briceño  _et al._ (2021)).

Explicit values of $\tilde{r}_{\mathfrak{n},a=\pi\pi,K\overline{K}}$ are shown
in Fig. 9, computed using the _reference coupled-channel_ amplitude
parameterization which successfully described the finite-volume spectra
presented in Section V. Levels below $K\overline{K}$ threshold have values of
$\tilde{r}_{\mathfrak{n},\pi\pi}$ in close agreement with the values computed
using the _reference elastic_ amplitude, as presented in Fig 8. The
corresponding $K\overline{K}$ components are observed to be quite uncertain,
as expected given the lack of constraint on $K\overline{K}$ well below its
kinematic threshold.888See Appendix A1 of Ref. Briceño  _et al._ (2021) for a
discussion of closed channels within this formalism. Above $K\overline{K}$
threshold, the $\pi\pi,K\overline{K}$ components are of similar magnitude –
that the $K\overline{K}$ components do not obviously dominate for levels lying
close to $K\overline{K}$ non-interacting energies is a result of having
already extracted a factor of the scattering matrix, $\mathcal{M}$, from the
quantity being finite-volume corrected, see Eq. III.

Figure 9:  Histograms for each energy level show the magnitudes of finite-
volume factors, $\big{|}\tilde{r}_{\mathfrak{n},a}\big{|}$, for $a=\pi\pi$
(blue bar) and $a=K\overline{K}$ (red bar) computed using the _reference
coupled-channel_ scattering amplitude. Non-interacting $\pi\pi$ (blue) and
$K\overline{K}$(red) energies are indicated by the horizontal dotted lines.

These values of $\tilde{r}_{\mathfrak{n}}$, combined with the lattice QCD
computed finite-volume form factors in Fig 7, will yield our infinite-volume
production form factors.

##  VIII Timelike pion form factor

The finite-volume form factor values extracted from lattice QCD computed
correlation functions were presented in Figure 7. Focussing initially on those
at energies lying below the $K\overline{K}$ threshold, we can correct with the
factor $\tilde{r}_{\mathfrak{n}}$ presented in Figure 8 to yield discrete
values of

$\mathcal{F}\quantity(s=E^{\star
2}_{\mathfrak{n}})=\frac{1}{\tilde{r}_{\mathfrak{n}}(L)}\mathcal{F}^{(L)}_{\mathfrak{n}}\,,$
(25)

as shown in the top panel of Fig. 10. This data is clearly consistent with
being describable by a smooth function of scattering energy, as demanded by
unitarity. When weighted with the _reference elastic_ scattering amplitude,
this data gives the pion form factor in the elastic timelike region (c.f. Eq.
5),

$f_{\pi}(s)=\frac{1}{k^{\star\,2}_{\pi\pi}}\mathcal{M}(s)\,\mathcal{F}(s)\,,$
(26)

which we plot in the middle panel of Figure 10. The form factor energy
dependence is clearly dominated by the presence of the narrow $\rho$ resonance
in $\mathcal{M}(s)$.

In order to accurately determine the uncertainty on the function
$\mathcal{F}(s)$, we must propagate the correlated errors on the finite-volume
correction factors, $\tilde{r}_{\mathfrak{n}}(L)$, and this poses an
implementation challenge. As seen in Figure 7, the finite-volume form factor
data has tiny fractional errors, such that the dominant uncertainty in Eq. 10
comes from $\tilde{r}_{\mathfrak{n}}(L)$ (see Figures 8, 9), and in practice,
the determined values of $\tilde{r}_{\mathfrak{n}}(L)$ for different energy
levels, $\mathfrak{n}$, are highly correlated with each other, owing to being
generated from scattering amplitude parameterizations featuring a relatively
small number of free parameters. Naively propagating such a high degree of
correlation into the $\chi^{2}$ minimization which will determine the
parameters in a parameterization of $\mathcal{F}(s)$ leads to results which do
not reflect the actual level of precision on the original energy level and
matrix-element data.

Our approach is to resample the ensembles of $\tilde{r}_{\mathfrak{n}}(L)$,
maintaining their computed mean and variance (as shown in Section VII), while
matching the distribution over lattice configurations of the energy level
values, $E_{\mathfrak{n}}$, from which they were obtained.999The only
exception to this are those levels lying ‘on-resonance’, i.e. the four levels
in the inset of Fig. 8. The data correlations for these levels are handled
slightly differently, reflecting the dominant sensitivity of
$\tilde{r}_{\mathfrak{n}}(L)$ for these levels to the parameter describing the
resonance coupling to $\pi\pi$. See Appendix C for details. This ensures that
the data correlation for the $\tilde{r}_{\mathfrak{n}}(L)$ inherits (up to
signs) the relatively mild data correlation of the corresponding
$E_{\mathfrak{n}}$, in accordance with a linearised approximation to error
propagation. Details of the implementation are presented in Appendix C.

Figure 10:  Elastic infinite-volume pion form-factor in the timelike region
determined from Eqns. 25, 26. Bottom panel shows the phase-shift for the
_reference elastic_ scattering amplitude, which by Watson’s theorem is also
the phase of the production amplitude.

Over the energy region of elastic scattering between the $\pi\pi$ threshold
and the $K\overline{K}$ threshold, $\mathcal{F}(s)$ can be parameterized using
a low-order polynomial in $s$,

$\mathcal{F}(s)/m_{\pi}^{2}=\sum_{n=0}c_{n}\cdot\left(\frac{s-s_{0}}{s_{0}}\right)^{n}\,,$
(27)

where we choose $\sqrt{s_{0}}=0.135\,a_{t}^{-1}=2.85\,m_{\pi}$, centering on
the resonance peak, for convenience.

As seen in Figure 11, fits linear or quadratic in $s$ can capture the observed
energy dependence, with the quadratic fit having parameters,

$c_{0}=$ | $0.2359\,(28)$ | $\begin{bmatrix}[r]1&-0.3&-0.6\\\ &1&0.7\\\ &&1\\\ \end{bmatrix}$
---|---|---
$c_{1}=$ | $0.265\,(26)$
$c_{2}=$ | $0.20\,(7)$
$\chi^{2}/N_{\text{dof}}=\frac{25.4}{17-3}=1.82$ ,
(28)

where parameter correlations are seen to not be excessive. We note here that
this fit has not imposed any additional constraints on the amplitude, such as
the fact that the form factor must take value 1 at $s=0$ to reflect the charge
of the pion. This is owing to the expectation that the simple _reference
elastic_ amplitude parameterization will not be suitable for extrapolation far
below the $\pi\pi$ threshold.

Figure 11: Polynomial descriptions of $\mathcal{F}(s)$ – linear in $s$ (el. I)
and quadratic in $s$ (el. II).

While the _timelike_ form factor, and its required finite-volume correction,
has been our focus thus far, the _spacelike_ form factor can also be extracted
on the same lattice. A consistent description of both energy regions can be
obtained by parameterizations with the appropriate analytic structure, e.g.
the Omnès function, as we describe in the next section.

##  IX Timelike & spacelike pion form factor

The pion electromagnetic form factor in the _spacelike_ region can be
extracted from calculations of three-point functions, $\langle
0|\Omega_{\pi}(\Delta t)\,\mathcal{J}(t)\,\Omega^{\dagger}_{\pi}(0)|0\rangle$,
where the so-determined matrix-elements do not require significant finite-
volume correction, as there is no two-hadron system that can go on-shell and
sample the spatial lattice boundary. A description of our procedures for
three-point function computation and analysis can be found in Ref.
Radhakrishnan  _et al._ (2022) with some details presented in Appendix H. In
this way, we determined the spacelike form factor over a range of $s<0$ (or
equivalently $Q^{2}=-s>0$ using the usual notation for virtualities) finding
it to undergo monotonic decrease from value $1$ at $s=0$ (imposition of which
serves to set the vector current renormalization factor as previously
mentioned).

As discussed in Section II, a _dispersive_ description of elastic $\pi\pi$
scattering can be constructed which is appropriate for analytic continuation
from the timelike region into the spacelike region. Use of this ‘Omnès’ form
allows us to reliably connect our determinations of the spacelike and timelike
form factors through a single amplitude. In brief, the relevant factor,
$\Omega(s)$, in Eq. 6 can be constructed using the elastic unitarity condition
on the production amplitude,
$f(s)-f^{*}(s)=2i\,\mathcal{M}^{*}(s)\rho(s)f(s)$, and the Schwarz reflection
principle to determine the discontinuity across the unitarity cut,
${f(s+i\epsilon)=e^{2i\delta(s+i\epsilon)}\,f(s-i\epsilon)}$. A dispersive
integral can be constructed which has the appropriate discontinuity,

$\log\Omega(s+i\epsilon)=\frac{s}{\pi}\int_{s_{\mathrm{thr}}}^{\infty}\\!ds^{\prime}\,\frac{\delta(s^{\prime})}{s^{\prime}\,(s^{\prime}-s-i\epsilon)}\,,$

where a subtraction has been introduced that normalizes $\Omega(0)=1$ while
also suppressing the dependence on the elastic phase-shift at high energies.

Our _reference elastic_ scattering amplitude provides an elastic phase-shift
from threshold up to $K\overline{K}$ threshold, and comparison with the
_reference coupled-channel_ amplitude shows that the $\pi\pi$ phase-shift
continues as given by the parameterization to somewhat higher energies. In
practice we will use this form in the integral above up to
$\sqrt{s_{a}}=1.2\cdot(2m_{K})$, and at energies higher than this a simple
parameterization,

$\delta(s>s_{a})=\pi-\big{(}\pi-\delta_{\mathrm{ref}}(s)\big{)}\frac{2}{1+\left(s/s_{a}\right)^{3/4}}\,,$

similar to the one proposed in Ref. Colangelo  _et al._ (2019).101010This form
enforces an asymptotic return of the phase-shift to $\pi$, correcting for the
threshold behavior built into the parameterization which ceases to be
appropriate at high energies. In fact Ref. Colangelo  _et al._ (2019) used
$\delta_{\mathrm{ref}}(s_{a})$ in place of $\delta_{\mathrm{ref}}(s)$ in this
expression but since the phase-shift of the _elastic reference_
parameterization has almost reached its asymptotic value at
$\sqrt{s_{a}}=1.2\cdot(2m_{K})$, there is no practical difference between the
two choices, and any small difference generated in $\Omega(s)$ can be absorbed
into the smooth function $\mathcal{F}_{\Omega}(s)$. In practice the Omnès
factor at low energies, where we require it, is not particularly sensitive to
the details of this continuation, owing to the subtraction suppressing the
contribution of high energies in the integral.

The finite-volume corrected form factor data from Figure 10 is presented in
Figure 12(b) with the Omnès factor, $\Omega(s)$, superimposed, where we
clearly see that the bulk of the energy dependence is captured, with need only
for a mild energy dependence in $\mathcal{F}_{\Omega}(s)$. Similarly Figure
12(a) shows our spacelike form factor data which also lies close to the
(analytically continued) $\Omega(s)$ for $s<0$.

(a) Spacelike region. (b) Timelike region.

Figure 12:  Pion form factor in the (a) spacelike and (b) timelike (as Fig 10)
regions. Superimposed are the Omnès function calculated from the _reference
elastic_ amplitude, and the Gounaris-Sakurai (GS) form determined by
describing the finite-volume spectra.

Also shown in these figures is a commonly-used parameterization known as the
Gounaris-Sakurai form Gounaris and Sakurai (1968), which effectively
corresponds to an elastic $K$-matrix pole describing the $\rho$ resonance,
along with a pole-subtracted Chew-Mandelstam phase-space. The parameters of
this amplitude were constrained by the lattice spectra in the elastic energy
region. The use of a dispersively improved phase-space removes the spurious
singularity otherwise present at $s=0$ and thus makes an extrapolation into
the space-like region somewhat plausible.

The function $\mathcal{F}_{\Omega}(s)$ obtained by dividing the Omnès factor
out from the timelike and spacelike form factor data can be parameterized and
a description of the lattice data obtained. Describing all data simultaneously
requires spanning a large energy region, and it proves helpful to make use of
a _conformal mapping_ of $s$ into a variable $z(s)$, defined as,

$z(s)=\frac{\sqrt{s_{c}-s_{0}}-\sqrt{s_{c}-s}}{\sqrt{s_{c}-s_{0}}+\sqrt{s_{c}-s}}\,,$
(29)

where this variable maps the entire complex plane of $s$, excluding the real
$s$ axis above $s_{c}$, into a unit disk in $z$ centered so that $z(s_{0})=0$.
Suitable choice of $s_{c}$ can reflect the fact that we expect there to be
cuts due to inelastic channels at higher energies, and $s_{0}$ can be selected
to conveniently distribute the data to be fitted around $z=0$. A polynomial in
$z(s)$ will, as desired, have no singularities in the timelike elastic
scattering region, or the spacelike region.

Figure 13:  Ratio of the form factor data, $|f_{\pi}(s)|$ (as presented in
Fig. 10), to the Omnès function, $\Omega(s)$, as a function of the conformal
map variable of Eq. 29. The curve shows the polynomial fit of Eq. 31.

Writing a parameterization of $\mathcal{F}_{\Omega}(s)$ as a low-order
polynomial in $z(s)$,

$\mathcal{F}_{\Omega}(s)=Q_{\pi^{+}}+\sum_{n=1}d_{n}\cdot\big{(}z(s)^{n}-z(0)^{n}\big{)}\,,$
(30)

where the constraint from the fixed electric charge of the pion,
$f_{\pi}(s\\!=\\!0)=\mathcal{F}_{\Omega}(s\\!=\\!0)=Q_{\pi^{+}}=1$, is simply
imposed, we can describe the lattice data across both spacelike and timelike
regions, as shown in Figure 13. With $\sqrt{s_{0}}=0.135\,a_{t}^{-1}$ and
$\sqrt{s_{c}}=0.22\,a_{t}^{-1}=1.27\cdot(2m_{K})$, the parameter values for a
quadratic description are,

$d_{1}=$ | $1.85\,(12)$ | $\begin{bmatrix}[r]1&0.8\\\ &1\\\ \end{bmatrix}$
---|---|---
$d_{2}=$ | $7.0\,(5)$
$\chi^{2}/N_{\text{dof}}=\frac{85.4}{37-2}=2.44$ .
(31)

The corresponding pion form factor for this description of
$\mathcal{F}_{\Omega}(s)$ is presented in Figure 14. The somewhat large
$\chi^{2}/N_{\text{dof}}$ is dominated by points in the spacelike region which
have high statistical precision. Lattice discretization effects which are
small in absolute terms, but of comparable size to the small statistical
errors could explain this observation.

(a) Energy dependence across spacelike and timelike regions. (b) Zoom over the
spacelike region.

Figure 14: Pion form factor across spacelike and timelike regions, and energy-
dependent description by Omnès modulated by the quadratic $F_{\Omega}(s)$ of
Eq. 31.

Figure 15 shows that over the timelike elastic scattering region, the Omnès
approach and the previous approach using Eq. 26 yield compatible descriptions
of the form factor. It is clear that in the region of elastic scattering, the
production amplitude can be determined with little systematic error. We now
turn to the problem of determining production amplitude above the lowest
inelastic threshold where we can produce real $K\overline{K}$ pairs.

Figure 15:  Timelike pion form-factor described by _reference elastic_
amplitude via Eq. 26 (green) or via Omnès form (orange). These descriptions
share a common phase given by the scattering amplitude phase-shift shown in
the bottom panel.

##  X Timelike coupled-channel production amplitudes

The additional constraint provided by the finite-volume matrix elements for
discrete energy levels above $K\overline{K}$ threshold allows us to access the
$\gamma\to K\overline{K}$ amplitude in addition to the $\gamma\to\pi\pi$
process at higher energies. The extension of the finite-volume production
formalism to multiple channels leads to a description of each finite-volume
form-factor, $\mathcal{F}_{\mathfrak{n}}^{(L)}$, as a linear combination of
the infinite-volume kaon and pion smooth production functions at $s=E^{\star
2}_{\mathfrak{n}}$, as indicated in Eq. 10. Since there is no longer a one-to-
one mapping from finite to infinite-volume, it is necessary to parameterize
the energy dependence of the production functions, and globally describe
multiple finite-volume form-factor values via a $\chi^{2}$ minimization, as
proposed in Ref. Briceño  _et al._ (2021).

We implement the effect of the correlated uncertainty on the finite volume
correction factors, $\tilde{r}_{\mathfrak{n},a}$, in terms of a ‘systematic’
contribution to the data covariance in the fit $\chi^{2}$,

$\displaystyle\chi^{2}=\sum_{\mathfrak{n},\mathfrak{m}}\Big{[}$
$\displaystyle\mathcal{F}^{(L)}_{\mathfrak{n}}-\big{(}\tilde{r}_{\mathfrak{n},\pi\pi}\mathcal{F}_{\pi\pi}(s)+\tilde{r}_{\mathfrak{n},K\overline{K}}\mathcal{F}_{K\overline{K}}(s)\big{)}\Big{]}$
$\displaystyle\cdot\big{(}C^{\mathrm{stat.}}+C^{\mathrm{syst.}}\big{)}^{-1}_{\mathfrak{n},\mathfrak{m}}$
$\displaystyle\cdot\Big{[}$
$\displaystyle\mathcal{F}^{(L)}_{\mathfrak{m}}-\big{(}\tilde{r}_{\mathfrak{m},\pi\pi}\mathcal{F}_{\pi\pi}(s)+\tilde{r}_{\mathfrak{m},K\overline{K}}\mathcal{F}_{K\overline{K}}(s)\big{)}\Big{]}\,,$

which augments the ‘statistical’ covariance of the
$\mathcal{F}^{(L)}_{\mathfrak{n}}$. This ‘systematic’ covariance is computed
using the resampled $\tilde{r}_{\mathfrak{n},a}$ introduced in Section VIII
whose data covariance inherits that of the energy levels. For each energy
level, only the dominant component, $a=\pi\pi$ or $K\overline{K}$ is used to
compute the covariance. The detailed implementation is presented in Appendix
C.

Figure 16:  Finite volume form factor (open black circles) with errorbars
modified to account for the fractional uncertainty on the finite-volume
correction factors, coming from the diagonal elements of the ‘systematic’
covariance, as described in the text. Global description using the _reference
coupled-channel_ scattering amplitude, and order-1 polynomial forms for
$\mathcal{F}_{\pi\pi}(s),\mathcal{F}_{K\overline{K}}(s)$ (green squares).

To describe the infinite volume form factors we use the parameterization of
Eq. 5, with the _reference coupled-channel_ scattering amplitude determined in
Sec. V.2. The smooth functions $\mathcal{F}_{a}$ are parameterized with low
order polynomials in $s$, analogous to Eq. 27,

$\mathcal{F}_{a}(s)/m_{\pi}^{2}=\sum_{n=0}^{N_{a}}\,h_{a,n}\cdot\quantity(\tfrac{s-s_{0}}{s_{0}})^{n}\,.$
(32)

We emphasize here that the function $\mathcal{F}_{\pi\pi}(s)$ obtained in the
coupled-channel case does not have to resemble the similarly named function in
the _elastic_ case, presented in Sec. VIII. This can be seen by considering
$f_{\pi}(s)$, which should be similar between the two cases, and which has a
representation,

$f_{\pi}=\frac{1}{k^{\star\,2}_{\pi\pi}}\mathcal{M}_{\pi\pi,{\pi\pi}}\quantity[\mathcal{F}_{\pi\pi}+\frac{k_{\pi\pi}^{\star}}{k_{K\overline{K}}^{\star}}\frac{\mathcal{M}_{\pi\pi,K\overline{K}}}{\mathcal{M}_{\pi\pi,{\pi\pi}}}\mathcal{F}_{K\overline{K}}]\,,$
(33)

in the coupled-channel case, where we have factored out the elastic $\pi\pi$
scattering in analogy to Eq. 26. Below the $K\overline{K}$ threshold, the
expression within square brackets is real and smooth (because the ratio of
scattering amplitude components does not posses the $\pi\pi$ branch cut), and
serves as the effective $\mathcal{F}$ in the prior elastic case.

Describing all the matrix element values previously presented in Figure 7
using finite-volume factors $\tilde{r}_{\mathfrak{n},a}(L)$ computed from the
_reference coupled-channel_ amplitude and smooth functions
$\mathcal{F}_{\pi\pi},\mathcal{F}_{K\overline{K}}$ each described by
polynomials of order 1, leads to the result summarized in Fig. 16. The
description given by the green points is observed to be in good qualitative
agreement with the constraining lattice data, with just a few isolated points
having significant deviations which produce a somewhat large
$\chi^{2}/N_{\text{dof}}$.

Examining Figure 15 we see that the points which are most discrepant in the
coupled-channel description were also somewhat discrepant in the elastic
description, and that the larger effect in the current case may reflect the
different approach for propagating the finite-volume factor uncertainty.

Table 1 presents variations of both the polynomial orders used in
$\mathcal{F}_{\pi\pi}(s)$, $\mathcal{F}_{K\overline{K}}(s)$, and the form used
for the coupled-channel scattering amplitude, where we consider a second
functional form of the $K$-matrix, given in Eq. 38, which includes an
additional linear term in $s$ with a set of $\gamma^{(1)}$ coefficients. The
reduction in the $\chi^{2}$ values for the varied scattering amplitude can be
entirely associated to the decreased precision of the finite-volume
corrections factors, which is implemented as a larger systematic error on the
data.

The final column of Tab. 1 shows an alternative choice for the determination
of the ‘systematic’ covariance, where the magnitude of the diagonal elements
are determined iteratively from the mean value of the model
$\sum_{a}\tilde{r}_{\mathfrak{n},a}\,\mathcal{F}_{a}$. This prescription was
inspired by the fitting procedure described in Ref. McRae and Blunden (2024),
more details can be found in that reference and in Appendix C. We found
consistent results from both methods to determine the ‘systematic’ covariance,
including the determination of the decay constant of the $\rho$-meson
described in the following section, and shown in Fig. 19.

Model | $N_{\pi\pi}$ | $N_{K\overline{K}}$ | $N_{\text{dof}}$ | $\chi^{2}/N_{\text{dof}}$
---|---|---|---|---
reference | $\gamma^{(1)}$ | ref. iter.
a | 1 | 0 | $32-3$ | 5.51 | 2.62 | 5.26
b | 0 | 1 | $32-3$ | 4.39 | 2.41 | 4.90
c | 1 | 1 | $32-4$ | 4.53 | 2.48 | 4.92
d | 2 | 0 | $32-4$ | 4.87 | 2.55 | 5.28
Table 1: Global description of full set of finite-volume matrix elements.
Variations in polynomial order in Eq. 32 used with _reference coupled-channel_
scattering amplitude (fifth column), an alternative coupled channel amplitude
with more parameter freedom, Eq. 38 (sixth column), and _reference coupled-
channel_ amplitude with a modified iterative fitting strategy described in
Appendix C (seventh column).

The timelike pion form factor for these descriptions is shown in Figure 17,
where we see only modest spread over the different parameterizations, and
rather close agreement with the result of Sec. IX described in terms of the
Omnès function over the elastic region. The bottom panel of Fig. 17 shows the
difference between the form factor phase and the $\pi\pi$ $P$-wave phase-
shift, which above the $K\overline{K}$ threshold need not
coincide111111Watson’s theorem equating the production phase to the scattering
phase-shift applies rigorously only in the elastic region. However, due to the
inelasticity in the scattering amplitude being close to 1, indicating little
$\pi\pi,K\overline{K}$ channel coupling, we expect a small phase difference,
and indeed we observe differences compatible with zero.

In Figure 18 we show the timelike pion and kaon (isovector) form-factors for
energies above the $K\overline{K}$ threshold, where we observe a much greater
sensitivity to the choice of scattering amplitude parameterization in the kaon
case. In particular the kaon form-factor using the “$\gamma^{(1)}$”
parameterization inherits a larger coupling to the $\rho$ from the scattering
amplitude, and more of a resemblance to the tail of the resonance.
Nevertheless, the kaon amplitudes broadly agree within the fairly large
uncertainties. They differ most at and below $K\overline{K}$ threshold, but
this is precisely where the finite-volume approach becomes unconstraining, as
the contribution of the production amplitude to the finite volume matrix
elements is exponentially suppressed below its threshold (see Appendix A of
Ref. Briceño  _et al._ (2021)). In the current case, this can be understood in
terms of the poorly constrained value of the coupling of the $K\overline{K}$
channel with the $\rho$-resonance. For a narrow resonance the value of $f_{a}$
at the resonance mass $m_{R}$ is approximated by,

$f_{a}=-i\sqrt{16\pi}\frac{\hat{c}_{a}}{m_{R}\Gamma_{R}}\mathcal{F}^{(L)}_{\mathfrak{n}}+\mathcal{O}\left(\tfrac{\Gamma_{R}}{m_{R}}\right)\,,$
(34)

and the value of $\hat{c}_{K\overline{K}}$ varies significantly between the
_reference coupled-channel_ and the “$\gamma^{(1)}$” scattering
parameterizations, with a greater value in the latter case. We report the
resonance parameters of the “$\gamma^{(1)}$” scattering amplitude in Appendix
B.

Figure 17: Timelike pion form factor using coupled-channel amplitude
parameterizations showing variations detailed in Table 1(blue curves). Also
shown the discrete elastic determinations of Section VIII and the
corresponding elastic Omnès description. The bottom panel shows the form
factor phase difference with respect to the $\pi\pi$ phase-shift. The points
in between panels show the energies of the levels used to constrain the form
factor energy dependence.

(a) Pion timelike form factor. (b) Kaon isovector timelike form factor.

Figure 18: Coupled-channel parameterization variations of Table 1 plotted over
the inelastic energy region.

##  XI Pion and rho parameters from form factors

The spacelike form factor is typically characterized by its slope at zero
momentum transfer, what is commonly referred to as the ‘charge radius’,
$\expectationvalue{r_{\pi}^{2}}^{1/2}$, where121212See Ref. Jaffe (2021), and
references therein, for a discussion about the proper interpretation of this
quantity when it is comparable to the Compton wavelength of the system.

$\expectationvalue{r_{\pi}^{2}}\equiv
6\evaluated{\frac{d}{ds}f_{\pi}(s)}_{s=0}\,.$ (35)

We determine this quantity from the parameterization in Eq. 31, obtaining
$\expectationvalue{r_{\pi}^{2}}^{1/2}=0.614(7)$ fm. This value sits in between
the experimental pion charge radius,
$\expectationvalue{r_{\pi}^{2}}^{1/2}=0.659(4)$ fm, and kaon charge radius,
$\expectationvalue{r_{K}^{2}}^{1/2}=0.560(31)$ fm, reported by the PDG Workman
and Others (2022), as we might expect given our lattice pion mass lying
between the physical pion and kaon masses. This result is also consistent with
previous analyses of timelike form factors determined at similar pion masses
in lattice QCD Andersen  _et al._ (2019); Erben  _et al._ (2020).

The timelike pion form-factor is dominated by the contribution of the $\rho$
resonance, and at the resonance pole, the $\pi\pi$ production amplitude
factorizes into the coupling to the $\pi\pi$ system and to the electromagnetic
current,

$\mathcal{H}^{\mu}_{\pi\pi,m}=\frac{\sqrt{16\pi}\,c_{\pi\pi}}{s_{R}-s}\,\big{\langle}\rho(\vec{P},m)\big{|}\mathcal{J}^{\mu}\big{|}0\big{\rangle}\,.$

The current-vector matrix element can be expressed in terms of a Lorentz-
invariant $\rho$-photon coupling, which we choose to parameterize in a
dimensionless form,131313Other normalizations have been used in the
literature, see for example Ref. Gasser and Leutwyler (1982) that uses
$f_{\rho}=\sqrt{2}\,m_{R}\,f_{V}$.

$\big{\langle}\rho(\vec{P},m)\big{|}\mathcal{J}^{\mu}\big{|}0\big{\rangle}\
=\epsilon^{\mu*}(P,m)\,f_{V}\,m_{R}^{2}\,,$ (36)

where $m_{R}^{2}=\text{Re}(s_{R})$.

One definition of the $\rho$ meson decay ‘partial width’ into an electron-
positron pair,

$\Gamma_{\rho\to
e^{+}e^{-}}=\frac{4\pi\alpha^{2}}{3}m_{R}\absolutevalue{f_{V}}^{2}\,,$

uses this coupling, and use of the PDG average for $\Gamma_{\rho\to
e^{+}e^{-}}$ can thus provide an estimate for $|f_{V}|$. A more consistent
approach, which was taken in Ref. Hoferichter  _et al._ (2017), would be to
describe the $e^{+}e^{-}\to\pi\pi$ cross-section energy dependence using the
infinite-volume amplitude parameterizations presented earlier, and
analytically continue them to the pole, yielding a complex-valued $f_{V}$.

Practical extraction of the $f_{V}$ coupling from our lattice-constrained
amplitudes varies slightly depending upon the form of the amplitude
construction. In the elastic case using $K$-matrix parameterizations, we may
simply use the $\pi\pi$ coupling from the scattering amplitude and the
singularity-free function, $\mathcal{F}(s)$, evaluated at the pole location
Hoferichter  _et al._ (2017),

$f_{V}=-\sqrt{\frac{4}{3}}\frac{1}{m_{R}^{2}}\sqrt{16\pi}\,\hat{c}_{\pi\pi}\,\mathcal{F}(s_{R})\,.$

In the case of the Omnès parameterization of Eq. 6 one needs to analytically
continue the Omnès function $\Omega$ into the unphysical Riemann sheet. This
can be achieved by multiplying $\Omega$ by the $S$-matrix evaluated in the
unphysical Riemann sheet where the scattering amplitude houses the resonance
pole.

In the coupled-channel case, the extracted coupling value is independent of
whether it is pulled from the pion or kaon form factor, in both cases being

$f_{V}=-\sqrt{\frac{4}{3}}\frac{1}{m_{R}^{2}}\,\sqrt{16\pi}\,\Big{(}\hat{c}_{\pi\pi}\mathcal{F}_{\pi\pi}(s_{R})-\hat{c}_{KK}\mathcal{F}_{KK}(s_{R})\Big{)}\,.$
(37)

Figure 19: Vector decay constant, $f_{V}$, for the $\rho$ meson. The first two
values correspond to the experimental extraction reported by the PDGWorkman
and Others (2022) and the dispersive analysis of Hoferichter et. al. [HKZ]
Hoferichter  _et al._ (2017). Subsequent values are the result of this
analysis, for only the elastic region (circles) and also including the coupled
channel region (squares).

In Figure 19 we summarize our determinations of $f_{V}$ using extractions from
elastic and coupled-channel parameterizations of various forms. The values are
seen to be in good agreement, and as expected for a calculation at an
unphysical quark mass, they differ somewhat from the experimental value
following from the PDG’s $\rho\to e^{+}e^{-}$ partial width. The errorbars
include the uncertainties on the resonance parameters, $c_{a}$ and $s_{R}$,
extracted from the scattering amplitude, as well as those of the smooth
functions $\mathcal{F}_{a}$, which are estimated from the variance of the
$\mathcal{F}^{(L)}_{\mathfrak{n}}$ data and the finite-volume correction
factors. The orange and purple points show two different strategies for
handling the error propagation from the finite-volume correction factors, as
discussed in Appendix C. Green points differ from the orange ones by the
functional form of the scattering amplitude used to describe the spectra and
to calculate the finite-volume factors, see Eq. 38 in Appendix B.

The value coming from the Omnès parameterization serves as a conservative
estimate at $m_{\pi}=284$ MeV,

$\displaystyle f_{V}$ $\displaystyle=0.224(6)\,e^{-i\pi\cdot
0.0418(27)}\quad[-0.09]\,,$ $\displaystyle{}\Gamma_{\rho\to e^{+}e^{-}}$
$\displaystyle=8.9(5)\text{ keV}\,,$

where the coupling is observed to be close to being real valued, with only
weak correlation between magnitude and phase.

A mild pion mass dependence of $f_{V}$ was observed in Ref. Jansen  _et al._
(2009) where the $\rho$ was studied using only two light quark flavors and
without explicit finite-volume correction (as it was either stable or very
narrow at the pion masses considered), and in Refs. Sun  _et al._ (2018);
DeGrand and Neil (2020) similarly for three flavors. However, a rigorous study
of the pion mass dependence of $f_{V}$ will require the implementation, as
done in this work, of the finite-volume correction to the lattice matrix
elements, and the analytic continuation to the resonance pole. These effects
become more significant as the pion mass decreases towards the physical point,
because the decay phase-space to pions, hence the width of the resonance,
increases, moving the resonance pole away from the real energy axis.

##  XII Summary and outlook

In this work we have presented a QCD calculation of the timelike form factor
of the pion in a kinematic region where the final state interactions are
sensitive to the rescattering of coupled $\pi\pi$ and $K\overline{K}$ systems.
The formalism correcting for the finite spatial volume of the lattice
necessarily gave us access to both the timelike pion and kaon form factors.

We also described _simultaneously_ the spacelike and the elastic timelike pion
form factor with a dispersive representation, the Omnès function, satisfying
the analyticity requirement of the amplitude in addition to the elastic
unitarity requirement. This provides an improvement over the Gounaris-Sakurai
form, which only satisfies the latter.

We found that considering a variety of functional forms for the energy
dependence of the form factors, when determined from lattice matrix elements
in the elastic timelike region, the spacelike plus elastic timelike region,
and the timelike region including above the $K\overline{K}$ threshold, all
gave consistent results for the decay constant of the $\rho$, the resonance
which dominates the timelike pion form factor at low energies. This
observation provides a consistency check on the coupled-channel transition
technology, motivating its application to more challenging cases of resonances
decaying into multiple two-meson channels, e.g. the $\pi_{1}$ studied in Ref.
Woss  _et al._ (2021).

The current calculation was restricted to the _isovector_ component of the
electromagnetic current, which is the only non-zero contribution in the case
of the pion form factor. However, in the case of the electromagnetic _kaon_
timelike form factor, an isoscalar component is allowed, and
phenomenologically this is known to be dominated at low energies by the
isoscalar vector poles of the $\omega$ and $\phi$ resonances Stamen  _et al._
(2022). A calculation of the isoscalar component, in a lattice with pion
masses $\lesssim 350$ MeV, would require a different technique to remove the
power-law finite-volume corrections, given that the isoscalar $K\overline{K}$
channel couples to a system of _three_ pions. Progress on the development of
this technology can be found in Refs. Briceño  _et al._ (2017); Hansen  _et
al._ (2021).

The techniques employed here can also be applied to study electroweak
processes involving multiple two-hadron final states. For example, to
calculate $CP$-violation in $D$-meson decays to $\pi\pi$ and $K\overline{K}$,
which could be then compared to the recent experimental determination of Ref.
Aaij  _et al._ (2019).

###### Acknowledgements.

We thank our colleagues within the Hadron Spectrum Collaboration for their
continued assistance. The authors acknowledge support from the U.S. Department
of Energy contract DE-SC0018416 at William & Mary, and contract DE-
AC05-06OR23177, under which Jefferson Science Associates, LLC, manages and
operates Jefferson Lab. This work contributes to the goals of the U.S.
Department of Energy _ExoHad_ Topical Collaboration, Contract No. DE-
SC0023598. The software codes Chroma Edwards and Joo (2005), QUDA Clark  _et
al._ (2010); Babich  _et al._ (2010a), QUDA-MG Clark  _et al._ (2016), QPhiX
Joó  _et al._ (2013), MG_PROTO Joó , and QOPQDP Osborn  _et al._ (2010);
Babich  _et al._ (2010b) were used for the computation of the quark
propagators. The authors acknowledge support from the U.S. Department of
Energy, Office of Science, Office of Advanced Scientific Computing Research
and Office of Nuclear Physics, Scientific Discovery through Advanced Computing
(SciDAC) program. Also acknowledged is support from the Exascale Computing
Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy
Office of Science and the National Nuclear Security Administration. This work
was also performed on clusters at Jefferson Lab under the USQCD Collaboration
and the LQCD ARRA Project. This research was supported in part under an ALCC
award, and used resources of the Oak Ridge Leadership Computing Facility at
the Oak Ridge National Laboratory, which is supported by the Office of Science
of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This
research used resources of the National Energy Research Scientific Computing
Center (NERSC), a DOE Office of Science User Facility supported by the Office
of Science of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231. The authors acknowledge the Texas Advanced Computing Center
(TACC) at The University of Texas at Austin for providing HPC resources. Gauge
configurations were generated using resources awarded from the U.S. Department
of Energy INCITE program at the Oak Ridge Leadership Computing Facility, the
NERSC, the NSF Teragrid at the TACC and the Pittsburgh Supercomputer Center,
as well as at the Cambridge Service for Data Driven Discovery (CSD3) and
Jefferson Lab.

## References

  * Amendolia  _et al._ (1986) S. R. Amendolia _et al._ (NA7), Nucl. Phys. B 277, 168 (1986).
  * Huber  _et al._ (2008) G. M. Huber _et al._ (Jefferson Lab), Phys. Rev. C 78, 045203 (2008).
  * Fang  _et al._ (2021) S.-s. Fang, B. Kubis, and A. Kupsc, Prog. Part. Nucl. Phys. 120, 103884 (2021).
  * Fujikawa  _et al._ (2008) M. Fujikawa _et al._ (Belle), Phys. Rev. D 78, 072006 (2008).
  * Colangelo  _et al._ (2019) G. Colangelo, M. Hoferichter, and P. Stoffer, JHEP 02, 006 (2019).
  * Aoki  _et al._ (2022) Y. Aoki _et al._ (Flavour Lattice Averaging Group (FLAG)), Eur. Phys. J. C 82, 869 (2022).
  * Luscher (1986a) M. Luscher, Commun. Math. Phys. 104, 177 (1986a).
  * Luscher (1986b) M. Luscher, Commun. Math. Phys. 105, 153 (1986b).
  * Briceno  _et al._ (2018a) R. A. Briceno, J. J. Dudek, and R. D. Young, Rev. Mod. Phys. 90, 025001 (2018a).
  * Lellouch and Luscher (2001) L. Lellouch and M. Luscher, Commun. Math. Phys. 219, 31 (2001).
  * Meyer (2011) H. B. Meyer, Phys. Rev. Lett. 107, 072002 (2011).
  * Feng  _et al._ (2015) X. Feng, S. Aoki, S. Hashimoto, and T. Kaneko, Phys. Rev. D 91, 054504 (2015).
  * Andersen  _et al._ (2019) C. Andersen, J. Bulava, B. Hörz, and C. Morningstar, Nucl. Phys. B 939, 145 (2019).
  * Erben  _et al._ (2020) F. Erben, J. R. Green, D. Mohler, and H. Wittig, Phys. Rev. D 101, 054504 (2020).
  * Mai  _et al._ (2023) M. Mai, U.-G. Meißner, and C. Urbach, Phys. Rept. 1001, 1 (2023).
  * Omnes (1958) R. Omnes, Nuovo Cim. 8, 316 (1958).
  * Luscher (1991a) M. Luscher, Nucl. Phys. B 354, 531 (1991a).
  * Luscher (1991b) M. Luscher, Nucl. Phys. B 364, 237 (1991b).
  * Rummukainen and Gottlieb (1995) K. Rummukainen and S. A. Gottlieb, Nucl. Phys. B 450, 397 (1995).
  * He  _et al._ (2005) S. He, X. Feng, and C. Liu, JHEP 07, 011 (2005).
  * Kim  _et al._ (2005) C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl. Phys. B 727, 218 (2005).
  * Christ  _et al._ (2005) N. H. Christ, C. Kim, and T. Yamazaki, Phys. Rev. D 72, 114506 (2005).
  * Fu (2012) Z. Fu, Phys. Rev. D 85, 014506 (2012).
  * Guo  _et al._ (2013) P. Guo, J. Dudek, R. Edwards, and A. P. Szczepaniak, Phys. Rev. D 88, 014501 (2013).
  * Hansen and Sharpe (2012) M. T. Hansen and S. R. Sharpe, Phys. Rev. D 86, 016007 (2012).
  * Briceno and Davoudi (2013) R. A. Briceno and Z. Davoudi, Phys. Rev. D 88, 094507 (2013).
  * Gockeler  _et al._ (2012) M. Gockeler, R. Horsley, M. Lage, U. G. Meissner, P. E. L. Rakow, A. Rusetsky, G. Schierholz, and J. M. Zanotti, Phys. Rev. D 86, 094513 (2012).
  * Leskovec and Prelovsek (2012) L. Leskovec and S. Prelovsek, Phys. Rev. D 85, 114507 (2012).
  * Briceno (2014) R. A. Briceno, Phys. Rev. D 89, 074507 (2014).
  * Dudek  _et al._ (2014) J. J. Dudek, R. G. Edwards, C. E. Thomas, and D. J. Wilson (Hadron Spectrum), Phys. Rev. Lett. 113, 182001 (2014).
  * Wilson  _et al._ (2015a) D. J. Wilson, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. D 91, 054008 (2015a).
  * Wilson  _et al._ (2015b) D. J. Wilson, R. A. Briceno, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. D 92, 094502 (2015b).
  * Dudek  _et al._ (2016) J. J. Dudek, R. G. Edwards, and D. J. Wilson (Hadron Spectrum), Phys. Rev. D 93, 094506 (2016).
  * Moir  _et al._ (2016) G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, and D. J. Wilson, JHEP 10, 011 (2016).
  * Briceno  _et al._ (2018b) R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. D 97, 054513 (2018b).
  * Woss  _et al._ (2019) A. J. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. D 100, 054506 (2019).
  * Woss  _et al._ (2021) A. J. Woss, J. J. Dudek, R. G. Edwards, C. E. Thomas, and D. J. Wilson (Hadron Spectrum), Phys. Rev. D 103, 054502 (2021).
  * Bulava  _et al._ (2024) J. Bulava _et al._ (Baryon Scattering (BaSc)), Phys. Rev. D 109, 014511 (2024).
  * Prelovsek  _et al._ (2021) S. Prelovsek, S. Collins, D. Mohler, M. Padmanath, and S. Piemonte, JHEP 06, 035 (2021).
  * Wilson  _et al._ (2024) D. J. Wilson, C. E. Thomas, J. J. Dudek, and R. G. Edwards (Hadron Spectrum), Phys. Rev. Lett. 132, 241901 (2024).
  * Woss  _et al._ (2020) A. J. Woss, D. J. Wilson, and J. J. Dudek (Hadron Spectrum), Phys. Rev. D 101, 114505 (2020).
  * Briceño  _et al._ (2015) R. A. Briceño, M. T. Hansen, and A. Walker-Loud, Phys. Rev. D 91, 034501 (2015).
  * Briceño and Hansen (2015) R. A. Briceño and M. T. Hansen, Phys. Rev. D 92, 074509 (2015).
  * Briceño  _et al._ (2021) R. A. Briceño, J. J. Dudek, and L. Leskovec, Phys. Rev. D 104, 054509 (2021).
  * Rodas  _et al._ (2023a) A. Rodas, J. J. Dudek, and R. G. Edwards, (2023a).
  * Rodas  _et al._ (2023b) A. Rodas, J. J. Dudek, and R. G. Edwards (Hadron Spectrum), Phys. Rev. D 108, 034513 (2023b).
  * Wilson  _et al._ (2019) D. J. Wilson, R. A. Briceno, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. Lett. 123, 042002 (2019).
  * Radhakrishnan  _et al._ (2022) A. Radhakrishnan, J. J. Dudek, and R. G. Edwards (Hadron Spectrum), Phys. Rev. D 106, 114513 (2022).
  * Peardon  _et al._ (2009) M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R. G. Edwards, B. Joo, H.-W. Lin, D. G. Richards, and K. J. Juge (Hadron Spectrum), Phys. Rev. D 80, 054506 (2009).
  * Woss  _et al._ (2018) A. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards, and D. J. Wilson, JHEP 07, 043 (2018).
  * Shultz  _et al._ (2015) C. J. Shultz, J. J. Dudek, and R. G. Edwards, Phys. Rev. D 91, 114501 (2015).
  * Jay and Neil (2021) W. I. Jay and E. T. Neil, Phys. Rev. D 103, 114502 (2021).
  * Gounaris and Sakurai (1968) G. J. Gounaris and J. J. Sakurai, Phys. Rev. Lett. 21, 244 (1968).
  * McRae and Blunden (2024) C. McRae and P. G. Blunden, Phys. Rev. C 109, 015503 (2024).
  * Jaffe (2021) R. L. Jaffe, Phys. Rev. D 103, 016017 (2021).
  * Workman and Others (2022) R. L. Workman and Others (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Gasser and Leutwyler (1982) J. Gasser and H. Leutwyler, Phys. Rept. 87, 77 (1982).
  * Hoferichter  _et al._ (2017) M. Hoferichter, B. Kubis, and M. Zanke, Phys. Rev. D 96, 114016 (2017).
  * Jansen  _et al._ (2009) K. Jansen, C. McNeile, C. Michael, and C. Urbach (ETM), Phys. Rev. D 80, 054510 (2009).
  * Sun  _et al._ (2018) W. Sun, A. Alexandru, Y. Chen, T. Draper, Z. Liu, and Y.-B. Yang ($\chi$QCD), Chin. Phys. C 42, 063102 (2018).
  * DeGrand and Neil (2020) T. DeGrand and E. T. Neil, Phys. Rev. D 101, 034504 (2020).
  * Stamen  _et al._ (2022) D. Stamen, D. Hariharan, M. Hoferichter, B. Kubis, and P. Stoffer, Eur. Phys. J. C 82, 432 (2022).
  * Briceño  _et al._ (2017) R. A. Briceño, M. T. Hansen, and S. R. Sharpe, Phys. Rev. D 95, 074510 (2017).
  * Hansen  _et al._ (2021) M. T. Hansen, F. Romero-López, and S. R. Sharpe, JHEP 04, 113 (2021).
  * Aaij  _et al._ (2019) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 122, 211803 (2019).
  * Edwards and Joo (2005) R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD),  _Lattice field theory. Proceedings, 22nd International Symposium, Lattice 2004, Batavia, USA, June 21-26, 2004_ , Nucl. Phys. Proc. Suppl. 140, 832 (2005).
  * Clark  _et al._ (2010) M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010).
  * Babich  _et al._ (2010a) R. Babich, M. A. Clark, and B. Joo, in  _SC 10 (Supercomputing 2010) New Orleans, Louisiana, November 13-19, 2010_ (2010).
  * Clark  _et al._ (2016) K. Clark, B. Joo, A. Strelchenko, M. Cheng, A. Gambhir, and R. Brower, in  _Proceedings of SC 16 (Supercomputing 2016) Salt Lake City, Utah, November 2016_ (2016).
  * Joó  _et al._ (2013) B. Joó, D. Kalamkar, K. Vaidyanathan, M. Smelyanskiy, K. Pamnany, V. Lee, P. Dubey, and W. Watson, in  _Supercomputing_, Lecture Notes in Computer Science, Vol. 7905, edited by J. Kunkel, T. Ludwig, and H. Meuer (Springer Berlin Heidelberg, 2013) pp. 40–54.
  * (71) B. Joó, “MG_PROTO: A Multigrid Library for QCD,” https://github.com/JeffersonLab/mg_proto/.
  * Osborn  _et al._ (2010) J. C. Osborn, R. Babich, J. Brannick, R. C. Brower, M. A. Clark, S. D. Cohen, and C. Rebbi,  _Proceedings, 28th International Symposium on Lattice field theory (Lattice 2010): Villasimius, Italy, June 14-19, 2010_ , PoS LATTICE2010, 037 (2010).
  * Babich  _et al._ (2010b) R. Babich, J. Brannick, R. C. Brower, M. A. Clark, T. A. Manteuffel, S. F. McCormick, J. C. Osborn, and C. Rebbi, Phys. Rev. Lett. 105, 201602 (2010b).
  * Aubert  _et al._ (2005) B. Aubert _et al._ (BaBar), Phys. Rev. D 71, 052001 (2005).
  * Aubert  _et al._ (2007) B. Aubert _et al._ (BaBar), Phys. Rev. D 76, 092005 (2007), [Erratum: Phys.Rev.D 77, 119902 (2008)].
  * Thomas  _et al._ (2012) C. E. Thomas, R. G. Edwards, and J. J. Dudek, Phys. Rev. D 85, 014507 (2012).
  * Segner  _et al._ (2023) A. Segner, A. Hanlon, A. Risch, and H. Wittig, in  _Proceedings of The 39th International Symposium on Lattice Field Theory — PoS(LATTICE2022)_, Vol. 430 (2023) p. 084.

##  Appendix A Spectrum sensitivity to use of $\pi\omega$ operators

The analysis of the $\pi\pi,K\overline{K}$ coupled channel system performed on
this lattice in principle is only justified up to the next inelastic
threshold, and the finite-volume approach we use is only rigorously correct up
to the lowest $n>2$ particle threshold. However, in practice, if the
$\pi\pi,K\overline{K}$ sector is decoupled from other channels as they open,
the analysis can be carried out to higher energies.

The $\pi\omega$ and $\pi\pi\pi\pi$ thresholds are statistically consistent
with each other, with the $\omega$ being effectively stable on this lattice.
Experimentally the $\pi\pi\pi\pi$ channel is found to be weakly interacting up
to 300 MeV above its threshold Aubert  _et al._ (2005), and we expect the lack
of $\pi\pi\pi\pi$-like operators to have a minimal impact in our analysis. A
similar situation occurs for the $\pi\pi\eta$ channel with a small cross-
section at energies below the effective threshold for $\rho\eta$, so that the
amplitude can be described by an isobar model Aubert  _et al._ (2007), and the
(negligible) impact of this channel in finite-volume was explored in a study
similar to the present one in Ref. Woss  _et al._ (2019).

Figure 20: Extracted finite-volume spectrum in the $[111]\,E_{2}$ irrep from
variational analysis of (left) correlation matrix including $\pi\omega$-like
and $\pi\phi$-like operators, and (right) excluding such operators. Dashed
lines indicate meson-meson thresholds, and solid lines the non-interacting
energy levels in this volume, with color coding as indicated below. The
histograms show the relative sizes of overlap factors for each state for each
operator in the basis, color coded as: subduced single-meson with
$J^{P}=1^{-}$ (dark blue), $J^{P}=1^{+}$ (brown), and $J^{P}=2^{-}$ (purple),
$\pi\pi$ (teal), $K\overline{K}$ (red), $\pi\omega$ (orange) and $\pi\phi$
(green). The spectrum from the left panel is reproduced on the right panel as
the grey points to aid comparison.

Since in-flight irreps subduced from non-zero helicity contain contributions
from $J^{P}=1^{+}$ as well as our desired $J^{P}=1^{-}$, the $\pi\omega$ can
appear in $S$-wave, and as such we might anticipate it having a significant
scattering amplitude141414indeed we expect it to resonate as the $b_{1}$ in
some energy region slightly higher than we have explored., and thus an impact
on the finite-volume spectrum.

We explicitly computed matrices of two-point correlations including operators
resembling both $\pi\omega$ and $\pi\phi$ pairs, and carried out a variational
energy-level extraction. An example of this is shown in Fig. 20, showing the
extracted discrete spectrum in the $[111]E_{2}$ irrep including such operators
(left), and excluding them (right).

As expected, the low-lying spectrum, well below the $\pi\omega$ threshold, is
unaffected by the inclusion of the extra operators whose non-interacting
energies lie much higher, and the overlaps,
$Z_{\mathfrak{n}}^{i}=\matrixelement{\mathfrak{n}}{\mathcal{O}_{i}^{\dagger}}{0}$,
for these states with those operators are observed to be negligible.

Above the $\pi\omega$ threshold we can identify by the operator-state overlaps
that some states are generated predominantly by the $\pi\omega$-like or
$\pi\phi$-like operators.

Because the variational analysis approach is able to separate contributions to
the matrix of correlation functions from multiple levels, even if they are
almost degenerate, we can observe that the spectrum on the left contains a
very precise energy level right on top of the $K_{110}\overline{K}_{100}$ non-
interacting energy with dominant overlap onto a $K\overline{K}$-like operator,
and which has energy statistically consistent with another three levels
overlapping dominantly with $\pi\omega$-like and $\pi\phi$-like operators.
This $K\overline{K}$-like level remains essentially unchanged in an analysis
that does not include the $\pi\omega$-like and $\pi\phi$-like operators,
justifying their removal for the spectrum determination, and their further
exclusion on the timelike form factor calculation. Comparable observations can
be made in the other irreps considered in this paper.

##  Appendix B Parameterization variations of the scattering amplitude

In order to determine the sensitivity of the extracted production amplitudes
to the detailed description of scattering, we used a variety of
parameterizations of the scattering amplitude to describe the finite-volume
spectrum extracted from lattice QCD computed correlation functions.

A general feature of the amplitudes that are able to successfully describe the
energy levels is a rapid phase-shift increase around
$E^{\star}=0.135\,a_{t}^{-1}$, which is efficiently parameterized by including
a pole term in the $K$-matrix (Eq. 13). Our variations concern whether this
pole is allowed to couple to the $K\overline{K}$ channel, what degree of
polynomial is added to the pole, and whether a simple, or dispersively
improved phase-space is used. Table 2 summarizes three variations used for
elastic analysis, and Table 3 summarizes 23 variations in the coupled-channel
case.

In the elastic case we found little variation in the scattering energy
dependence from different functional forms, as illustrated in Fig. 4. In the
coupled-channel analysis we find the parameterizations listed in Tab. 3 in
overall agreement, with slight variations seen on the energy dependence of the
inelasticity, but statistically consistent values of the phase-shifts
$\delta^{\pi\pi}_{1}$ and $\delta^{K\overline{K}}_{1}$. The variations for the
resonance pole position $s_{R}$ and coupling to the $\pi\pi$ channel
$c_{\pi\pi}$ fall within the statistical uncertainty of the data. However, as
mentioned in the main text, the coupling to the $K\overline{K}$ channel is
poorly determined in each parameterization, and some models differ from each
other at the level of up to two standard deviations.

To illustrate the results with largest variation we select two of these
coupled-channel amplitudes for propagation into the production amplitude
analysis. The first is the reference amplitude of Eq. 15, the first entry of
Tab. 3, while the second, which we label “$\gamma^{(1)}$”, has functional
form,

$\mathcal{K}_{ab}(s)=\frac{g_{a}\,g_{b}}{m^{2}-s}+\gamma_{ab}^{(0)}+\gamma_{ab}^{(1)}s\,,$
(38)

corresponding to the tenth entry of Table 3. The resonance pole of the
“$\gamma^{(1)}$” parameterization is located at

$a_{t}\sqrt{s_{R}}=0.1327(5)-\tfrac{i}{2}0.0096(4)\,,$ (39)

consistent with Eq. 18, while the channel couplings are,

$\displaystyle a_{t}\,c_{\pi\pi}$ $\displaystyle=0.0424(9)\,e^{-i\pi\cdot
0.055(4)}\,,$ $\displaystyle a_{t}\,c_{K\overline{K}}$
$\displaystyle=0.097(24)\,e^{i\pi\cdot 0.484(10)}\,,$ (40)

where we observe a noticeable variation on $c_{K\overline{K}}$ as compared to
the reference amplitude case in Eq. 19. Even with this variation of the
coupling, the value of $f_{V}$ determined with Eq. 37 shows little dependence
on the chosen scattering amplitude, see Fig. 19.

Parameterization | $N_{\text{pars}}$ | $\chi^{2}/N_{\text{dof}}$
---|---|---
Relativistic Breit-Wigner | 2 | 1.23
$\mathcal{K}=\frac{g^{2}}{m^{2}-s}$ (Gounaris-Sakurai) | 2 | 1.29
$\mathcal{K}=\frac{g^{2}}{m^{2}-s}+\gamma$ (With $-i\rho$ phase space) | 3 | 1.23
$\mathcal{K}=\frac{g^{2}}{m^{2}-s}+\gamma$ | 3 | 1.21
Table 2: $P$-wave elastic amplitude parameterization variations. Second and fourth entries use Chew-Mandelstam phase-space Parameterization | Restrictions | $N_{\text{pars}}$ | $\chi^{2}/N_{\text{dof}}$
---|---|---|---
$\mathcal{K}_{ab}=\frac{g_{a}g_{b}}{m^{2}-s}+\gamma_{ab}$ | – | 6 | 1.10
$g_{K\overline{K}}=0$ | 5 | 1.07
$\gamma_{K\overline{K},K\overline{K}}=0$ | 5 | 1.07
$\gamma_{\pi\pi,K\overline{K}}=0$ | 5 | 1.08
$\gamma_{\pi\pi,\pi\pi}=0$ | 5 | 1.47
$\mathcal{K}_{ab}=\frac{g_{a}g_{b}}{m^{2}-s}+\gamma^{(1)}_{ab}s$ | – | 6 | 1.05
$g_{K\overline{K}}=0$ | 5 | 1.03
$\gamma^{(1)}_{K\overline{K},K\overline{K}}=0$ | 5 | 1.02
$\gamma^{(1)}_{\pi\pi,K\overline{K}}=0$ | 5 | 1.02
$\mathcal{K}_{ab}=\frac{g_{a}g_{b}}{m^{2}-s}+\gamma_{ab}(s)$ $\gamma_{ab}(s)=\gamma^{(0)}_{ab}+\gamma^{(1)}_{ab}s$ | – | 9 | 1.05
$g_{K\overline{K}}=0$ | 8 | 1.11
$\gamma^{(1)}_{K\overline{K},K\overline{K}}=0$ | 8 | 1.07
$\gamma^{(1)}_{\pi\pi,K\overline{K}}=0$ | 8 | 1.06
$\gamma^{(1)}_{\pi\pi,\pi\pi}=0$ | 8 | 1.08
$\mathcal{K}_{ab}=\frac{g_{a}(s)g_{b}(s)}{m^{2}-s}+\gamma_{ab}$ $g_{a}(s)=g_{a}^{(0)}+g_{a}^{(1)}s$ | – | 8 | 1.15
$g_{\pi\pi}^{(1)}=0$ | 7 | 1.11
$\gamma_{\pi\pi,\pi\pi}=0$, $\gamma_{\pi\pi,K\overline{K}}=0$ | 6 | 1.10
$\gamma_{\pi\pi,\pi\pi}=0$, $\gamma_{K\overline{K},K\overline{K}}=0$ | 6 | 1.11
|  $\gamma_{\pi\pi,K\overline{K}}=0$, $\gamma_{K\overline{K},K\overline{K}}=0$ | 6 | 1.10
$\mathcal{K}_{ab}=\frac{g_{a}g_{b}}{m^{2}-s}+\gamma_{ab}$ $I_{ab}=-i\,\delta_{ab}\,\rho_{a}$ | – | 6 | 1.12
$g_{K\overline{K}}=0$ | 5 | 1.11
$\gamma_{K\overline{K},K\overline{K}}=0$ | 5 | 1.12
$\gamma_{\pi\pi,K\overline{K}}=0$ | 5 | 1.11
Table 3: $P$-wave coupled-channel amplitude parameterization variations.

##  Appendix C Data correlation between finite-volume correction factors

The $\tilde{r}_{\mathfrak{n}}(L)$ factors are calculated from scattering
amplitudes described by a small number of parameters determined as ensembles
over the set of lattice configurations. This produces a high degree of
correlation among them, which poses an implementation challenge to properly
propagate their uncertainty into the determination of production amplitudes.
This high degree of correlation does not reflect the relatively modest amount
of correlation between the energy levels in the finite-volume spectrum which
constrain the scattering amplitude, as shown in Fig. 21.

Figure 21:  Data correlation matrix for the energy-level values,
$E^{\star}_{\mathfrak{n}}$, with rows and columns ordered by increasing energy
from left to right, and from top to bottom, respectively.

In place of the data correlation for $\tilde{r}_{\mathfrak{n}}(L)$ computed
using the amplitude parameterization, we adopt one which inherits the energy-
level correlation, with the following motivation: In a linearised
approximation to error propagation, if the Jacobian were known, we could
calculate the covariance of the $\tilde{r}_{\mathfrak{n}}(L)$ using

$\sigma(\tilde{r}_{\mathfrak{n}},\tilde{r}_{\mathfrak{m}})=\sum\nolimits_{\mathfrak{n}^{\prime}\mathfrak{m}^{\prime}}\partialderivative{\tilde{r}_{\mathfrak{n}}}{E_{\mathfrak{n}^{\prime}}}\,\sigma(E_{\mathfrak{n}^{\prime}},E_{\mathfrak{m}^{\prime}})\,\partialderivative{\tilde{r}_{\mathfrak{m}}}{E_{\mathfrak{m}^{\prime}}}\,.$
(41)

However, to a good approximation, we expect the finite-volume correlation
factor computed with Eqns. VII and 24 to be influenced predominantly by the
_local_ behavior of $\mathcal{M}$ and $F$ around $E_{\mathfrak{n}}$. This
would imply that the Jacobian between $\tilde{r}_{\mathfrak{n}}(L)$ and the
energies, $\\{E_{\mathfrak{n}^{\prime}}\\}$, is largely dominated by the
diagonal elements, $\mathfrak{n}=\mathfrak{n}^{\prime}$.

For instance, in an ideal scenario, where the lattice spectrum was determined
for a large number of closely spaced volumes, $L$, and lattice irreps,
$\Lambda$, we would have access to a high density of states per energy unit.
From this spectrum we could determine the scattering amplitude
$\mathcal{M}(E)$ algebraically at a similarly high density, and in that case,
the calculation of $\tilde{r}_{\mathfrak{n}}(L)$ would not require an explicit
parameterization of $\mathcal{M}$, and the uncertainty of each finite-volume
correction factor would only be correlated to the uncertainty of the nearby
energy levels.

Based on this observation, we approximate the Jacobian of Eq. 41 by a diagonal
matrix
$\partialderivative{\tilde{r}_{\mathfrak{n}}}{E_{\mathfrak{n}^{\prime}}}\propto\delta_{\mathfrak{n},\mathfrak{n}^{\prime}}$.
In practice, we implement this procedure by shifting and reweighing the
jackknife ensemble of each energy level to match the mean and error of the
respective finite-volume correction factor according to the values computed
using the parameterization of $\mathcal{M}$ and Eqns. VII, 24. With this
prescription it is not necessary to compute the magnitude of the Jacobian
elements, however we still require knowledge of their sign. These are
extracted from the correlation between the ‘model’ energy, i.e. the solution
of Eq. 7, and $\tilde{r}_{\mathfrak{n}}(L)$, both extracted from a
parameterization of $\mathcal{M}$ 151515Support for this prescription is found
empirically in the fact that the absolute value of this correlation value is
consistent with 1 for most energy levels., these are shown on the top row of
Figure 22.

An exception to this behavior are those levels with energies very close to the
resonance mass, i.e. the four levels in the inset of Fig. 8. As mentioned
before, the value of $\tilde{r}_{\mathfrak{n}}(L)$ for these levels is equal
to the resonance coupling, modulo rather small finite-width corrections, and
we observe that for these levels $\tilde{r}_{\mathfrak{n}}(L)$ is mostly
uncorrelated with the respective energy value. This is illustrated in Figure
22, where we see that their strongest correlation (for the _reference elastic_
parameterization) is with the $g$ parameter. From the spectra of Fig. 3 we
note that these levels are the farthest from any non-interacting energy level,
and as such they are the least sensitive to finite-volume effects.

Figure 22: Data correlation between $\tilde{r}_{\mathfrak{n}}$ and the
respective Lüscher energy solution $E^{\star}_{\mathfrak{n}}$, as well as the
scattering parameters in the _reference elastic_ parameterization.

To account for this observation, for these levels we adopt a modified
procedure, where instead of only using the energy level data, we add a linear
combination of the fluctuations on the scattering parameters, with the
parameter $g$ having the highest contribution. For
$\tilde{r}_{\mathfrak{n}}(L)$ for ‘on-resonance’ levels, we form a linear
combination of the ensembles of the scattering parameters and the
corresponding energy level,

$\vec{\tilde{r}}_{\mathfrak{n}}(L)=x_{E}\,\vec{E}_{\mathfrak{n}}+x_{g}\,\vec{g}+x_{m}\,\vec{m}+x_{\gamma}\,\vec{\gamma}=\sum\nolimits_{i}x_{i}\,\vec{p}_{i}.$
(42)

where we use the vector notation $\vec{p}_{i}$ to emphasize that we are
manipulating the quantity $a$ in the ensemble configuration space. The
quantities $\vec{a}$ have been previously shifted and reweighed, so that they
have zero mean and unit covariance, i.e. $\expectationvalue{\vec{p}_{i}}=0$
and $\expectationvalue{\vec{p}_{i}\,^{2}}=1$.

To obtain the coefficients $x_{i}$ we demand that the correlation between the
generated $\tilde{r}_{\mathfrak{n}}$ ensemble and each of the ${p}_{i}$
quantities reproduces the correlation
$\text{corr}(p_{i},\tilde{r}_{\mathfrak{n}}(L))$ shown in Fig. 22.
Mathematically, this translates into the system of equations

$\displaystyle\sum\nolimits_{j}\expectationvalue{\vec{p}_{i}\,\vec{p}_{j}}x_{j}$
$\displaystyle=\text{corr}(p_{i},\tilde{r}_{\mathfrak{n}}(L))\,,$ (43)

where the matrix $\expectationvalue{\vec{p}_{i}\,\vec{p}_{j}}$ is calculated
using the corresponding lattice energy level and the scattering parameters.
Equation 43 represents a linear system of equations from which the
coefficients $x_{i}$ can be easily obtained. We find that the value of
$\expectationvalue{\vec{\tilde{r}}_{\mathfrak{n}}(L)^{2}}$ extracted in this
way is equal to unity, providing a self-consistency check of this procedure.
Finally we reweight and shift $\vec{\tilde{r}}_{\mathfrak{n}}(L)$ to obtain a
jackknife ensemble with the mean and error of $\tilde{r}_{\mathfrak{n}}(L)$.

It is not immediately obvious how to extend the procedure we just described
from the elastic analysis to the coupled channel case, given that multiple
components $\tilde{r}_{\mathfrak{n},a}(L)$ are associated with each level.
Furthermore, we need the correlation between the components of
$\tilde{r}_{\mathfrak{n}}$ to reproduce the correlation between different
elements of the scattering amplitude. For instance, to stay consistent with
unitarity the following relation (see Appendix A1 of Ref. Briceño  _et al._
(2021)) needs to be satisfied in the elastic energy region:

$\frac{\tilde{r}_{\mathfrak{n},K\overline{K}}}{\tilde{r}_{\mathfrak{n},\pi\pi}}=\frac{k_{\pi\pi}^{\star}\,\mathcal{M}_{\pi\pi,K\overline{K}}}{k_{K\overline{K}}^{\star}\,\mathcal{M}_{\pi\pi,\pi\pi}}\,.$
(44)

We instead choose to estimate the impact of the uncertainty of the finite-
volume correction by means of a “systematic” covariance, $C^{\mathrm{syst.}}$,
applied to the finite-volume matrix elements
$\mathcal{F}_{\mathfrak{n}}^{(L)}$. We add this systematic covariance matrix
to the statistical covariance matrix of $\mathcal{F}_{\mathfrak{n}}^{(L)}$ to
obtain a total covariance, which in turn is used for the $\chi^{2}$ function
used to determine the smooth functions $\mathcal{F}_{a}(s)$.

To calculate the magnitude of the diagonal elements of this covariance we
first note that for levels below the $K\overline{K}$ threshold the component
$\tilde{r}_{\mathfrak{n},K\overline{K}}$ cannot be precisely determined (by
the nature of the finite-volume formalism). Therefore, we assign the relative
uncertainty of the “relevant” $\tilde{r}_{\mathfrak{n},\pi\pi}$ as the
systematic relative uncertainty of the corresponding
$\mathcal{F}_{\mathfrak{n}}^{(L)}$. We note that following this procedure for
the elastic analysis we obtain results consistent with the alternative
procedure of replacing the correlation of the finite-volume correction factors
by that of the energies and propagating the error over jackknife.

Above the $K\overline{K}$ threshold we have a weakly interacting system, where
each state in the spectrum is close to a non-interacting energy. This allows
us to unambiguously assign a dominant component
$\tilde{r}_{\mathfrak{n},a}(L)$ to each level. Furthermore, following the
expectations outlined about the correlation between the finite-volume
correction factor and the energy value, we find that out of the two, this
component is always the most correlated with the energy solution of the
quantization condition. Once we picked the dominant component $a$, we again
assign the relative uncertainty of $\tilde{r}_{a,\mathfrak{n}}$ as the
“systematic” relative uncertainty of the corresponding
$\mathcal{F}_{\mathfrak{n}}^{(L)}$.

Finally, to obtain off-diagonal elements of $C^{\mathrm{syst.}}$, showing the
correlation of the finite-volume correction across the spectrum, we mimic the
prescription used in the elastic analysis. This means that we will once again
use the correlation matrix of the energy levels, or the appropriately
constructed linear combinations for the levels close to the resonance. This
correlation matrix is multiplied by the magnitude of the diagonal elements
determined above to obtain $C^{\mathrm{syst.}}$.

We also explored an alternative procedure to determine the diagonal elements
of $C^{\mathrm{syst.}}$ to corroborate the consistency of our prescription.
This second option was inspired by the iterative fit method described in Ref.
McRae and Blunden (2024). In this case we use the value of the “model”,
$\sum_{a}\tilde{r}_{a,\mathfrak{n}}(L)\mathcal{F}_{a}(E_{\mathfrak{n}}^{\star
2})$, to multiply the relative uncertainty of the finite-volume correction
factor and obtain an absolute systematic error for each diagonal element of
$C^{\mathrm{syst.}}$. In the first iteration we pick the solution from our
previous prescription to determine $C^{\mathrm{syst.}}$ and minimize the
$\chi^{2}$ with this systematic covariance. This yields new values for
$\mathcal{F}_{a}(s)$, which are then used to recompute $C^{\mathrm{syst.}}$
and repeat the minimization. We find that this process converges after a few
iterations. The results of this iterative procedure show consistency with the
previous choice of using $\mathcal{F}_{\mathfrak{n}}^{(L)}$ to determine the
magnitude of the systematic errors. This can be seen in the very similar
values of the $\chi^{2}$ for each option, shown in Table 1, and the
determination of the decay constant, $f_{V}$, shown in Figure 19.

##  Appendix D Subduction of kinematic factor $K(\Lambda)$

To obtain the kinematic factor $K(\Lambda)$ we need to project the components
of the Lorentz vector $K^{\mu}_{m}$ of Eq. 2 into a basis of definite
helicity, to then be able to apply the subduction coefficients
$\mathcal{S}^{\lambda}_{\Lambda}$ of Table II of Ref. Thomas  _et al._ (2012).
To illustrate this procedure, we will focus on the polarization vector
$\epsilon^{\mu}(P,m)$, and include all other scalar factors at the end.

We begin by calculating the components of the polarization vector in the frame
of momentum $P^{\mu}$, with a final hadron pair in a state of definite
helicity $\lambda$. This is achieved by boosting the rest-frame polarization
vector, with components $(0,\epsilon^{i}_{\lambda})^{\mu}$, to the frame of
$P^{\mu}$ in a two-step process: first a boost is applied in the $z$-direction
with velocity $\beta=\absolutevalue{\mathbf{P}}/P^{0}$, followed by a rotation
$R_{\hat{\mathbf{P}}}$ from the $z$-axis to the direction of the spatial part
of $P^{\mu}$. Applying these transformations to the rest-frame polarization
vector we obtain

$\epsilon^{\mu}(P,\lambda)=\begin{pmatrix}\gamma\beta\,\delta_{\lambda,0}\\\
\epsilon^{i}_{\lambda}(\hat{\mathbf{P}})\big{(}\gamma\,\delta_{\lambda,0}+\delta_{\lambda,\pm}\big{)}\end{pmatrix}\,,$

where
$\epsilon^{i}_{\lambda}(\hat{\mathbf{P}})\equiv[R_{\hat{\mathbf{P}}}]^{i}{}_{j}\epsilon^{j}_{\lambda}$,
and $\gamma$ is the relativistic factor associated to $\beta$.

Once we know the cartesian components of the polarization vector, we follow
Ref. Thomas  _et al._ (2012) to project the different $\mu$ components of
$K^{\mu}_{\lambda}$, according to their properties under spatial rotations,
into a helicity basis with spin-components $\hat{\lambda}$. This process
converts $K^{\mu}_{\lambda}$ into a diagonal matrix in the
$\\{\hat{\lambda},\lambda\\}$ space, with the temporal component only spanning
elements with zero helicity, while the spatial components span helicities $0$
and $\pm 1$,

$\displaystyle K^{0}_{\lambda}$ $\displaystyle\to
K^{[0]}_{\hat{\lambda}\lambda}=\sqrt{\tfrac{4}{3}}\,\gamma\beta\,\delta_{\hat{\lambda},0}\,\delta_{\lambda,0}\,\,,$
$\displaystyle K^{i}_{\lambda}$ $\displaystyle\to
K^{[i]}_{\hat{\lambda}\lambda}=\sqrt{\tfrac{4}{3}}\,\quantity(\gamma\delta_{\hat{\lambda},0}\,\delta_{\lambda,0}+\delta_{\hat{\lambda},\pm}\,\delta_{\lambda,\pm})\,.$

Finally, the kinematic factor $K(\Lambda)$ is calculated by applying the
appropriate subduction coefficients

$\sum\nolimits_{\hat{\lambda},\lambda}\mathcal{S}_{\Lambda}^{\hat{\lambda}}\,K^{[\mu]}_{\hat{\lambda}\lambda}\,\mathcal{S}_{\Lambda^{\prime}}^{\lambda}=\delta_{\Lambda\Lambda^{\prime}}K(\Lambda)\,,$

the result of which listed in Table 4. This last step is actually trivial
because the subduction matrices are unitary, and the matrix
$K^{[\mu]}_{\hat{\lambda}\lambda}$ is proportional to the identity in the
subspace where the subduction matrices are non-zero.

$\Lambda(\mathcal{J}^{\mu})$ | $T_{1}^{-}(\mathcal{J}^{i})$ | $A_{1}(\mathcal{J}^{0})$ | $A_{1}(\mathcal{J}^{i})$ | $E_{2},B_{1},B_{2}(\mathcal{J}^{i})$
---|---|---|---|---
$K(\Lambda)$ | $\sqrt{\tfrac{4}{3}}$ | $\sqrt{\tfrac{4}{3}}\gamma\beta$ | $\sqrt{\tfrac{4}{3}}\gamma$ | $\sqrt{\tfrac{4}{3}}$
Table 4: Kinematic factor of the subduced finite volume matrix elements into
irrep $\Lambda$ from the temporal, $\mathcal{J}^{0}$, or spatial,
$\mathcal{J}^{i}$, components of a vector current.

##  Appendix E Timeslice fit algorithm of the finite-volume matrix elements

Our approach to extracting a single value of the finite-volume matrix element
from timeslice correlation functions, as presented in Fig. 6, follows the
‘model averaging’ prescription suggested by Jay and Neil Jay and Neil (2021).

Insisting that a minimal number of timeslices $N_{t}$ must be used in any fit,
we perform fits to a constant $C$ with varying values $t_{\text{min}}$ but a
fixed value of $t_{\text{max}}$. For each such fit we compute a version of the
Akaike Information Criterion (AIC), by combining the correlated $\chi^{2}$ and
the number of degrees of freedom,
$N_{\text{dof}}=N_{\mathrm{data}}-N_{\mathrm{params}}$, according to
${w=\exp[-(\chi^{2}/2-N_{\text{dof}})]}$.

We label $t^{\star}$ the value of $t_{\text{min}}$ having the maximum value of
$w$, which according to the criterion identifies this range as the dataset for
a given $t_{\text{max}}$ best described by a constant. Subsequently fits are
attempted for all $t_{\text{min}}<t^{\star}$ using a constant plus an
exponential.

We repeat this procedure for all values of $t_{\text{max}}$ allowed by
$N_{t}$, yielding a large number of fits to a constant or a constant plus an
exponential for many timeslice fit windows. Each fit, $\alpha$, yields a value
and a statistical error for the constant, $C_{\alpha}\pm\sigma_{C,\alpha}$,
which can be ranked according to their respective value of AIC weight
$w_{\alpha}$, and a _model average_ computed from a weighted average of a
representative set of the fits with the largest weight,

$\displaystyle C$
$\displaystyle=\frac{\sum_{\alpha}w_{\alpha}C_{\alpha}}{\sum_{\beta}w_{\beta}}\,,$
(45) $\displaystyle\sigma_{C}^{2}$
$\displaystyle=\frac{\sum_{\alpha}w_{\alpha}\sigma_{C,\alpha}^{2}}{\sum_{\beta}w_{\beta}}+\frac{\sum_{\alpha\beta}w_{\alpha}w_{\beta}(C_{\alpha}-C_{\beta})^{2}}{2\quantity(\sum_{\beta}w_{\beta})^{2}}\,.$
(46)

To make this process compatible with the jackknife resampling technique used
for error propagation, we employ the following prescription, based on the
method presented in Ref. Segner  _et al._ (2023). For each fit we can obtain
an ensemble $\\{C_{\alpha,i}\\}$, where the index $i$ indicates an entry in
the jackknife ensemble of model $\alpha$. Then we can define the model
averaged ensemble,

$C_{i}=\frac{\sum_{\alpha}w_{\alpha}\,C_{\alpha,i}}{\sum_{\beta}w_{\beta}}\,.$

The average of this ensemble is equal to the value given in Eq. 45, but the
variance is bounded from above by the value of Eq. 46 (typically we find it to
be a few percent lower). We can fix this mismatch by adding Gaussian noise to
$C_{i}$, as long as the random variable $\eta_{i}$ has zero mean and a
variance such that the variance of

$\overline{C}_{i}=C_{i}+\eta_{i}\,$ (47)

matches Eq. 46. The Gaussian noise variables, $\eta_{i}$, are drawn from a
multidimensional uncorrelated Gaussian distribution, such that this
prescription does not impact the covariance among different matrix elements.
Reference Jay and Neil (2021) does not address the model average of the
covariance among different variables, but we assume that it is best to not
modify it. We employ the ensembles of Eq. 47 to calculate the form-factors
presented in this work.

##  Appendix F Late time pollution to matrix elements

The use of optimized operators at the source in correlation functions having
the vector current at the sink allows us access to finite-volume matrix
elements for states above the ground state owing to their overlap with other
states in the spectrum being highly suppressed,
$\matrixelement{\mathfrak{m}}{\Omega_{\mathfrak{n}}^{\dagger}(0)}{0}\ll 1$,
for $\mathfrak{m}\neq\mathfrak{n}$. This is a powerful technique Shultz  _et
al._ (2015), but we do encounter cases where a hierarchy of current matrix
elements for different states can compensate for this suppression, i.e. where
${\mathcal{F}_{\mathfrak{m}}^{(L)}\cdot\matrixelement{\mathfrak{m}}{\Omega_{\mathfrak{n}}^{\dagger}(0)}{0}\sim\mathcal{F}_{\mathfrak{n}}^{(L)}}$.

For those cases the ratio of correlation functions of Eq. 22 has non-
negligible additive contributions from states other than $\mathfrak{n}$ of the
form,

$\varepsilon_{\mathfrak{m},\mathfrak{n}}(t)=\varepsilon_{\mathfrak{m},\mathfrak{n}}\,e^{(E_{\mathfrak{n}}-E_{\mathfrak{m}})t}\,,$
(48)

such that states lighter than $\mathfrak{n}$ will cause a rising time-
dependence, visible at late times, that is not accounted for in our default
timeslice fitting form.

Use of a GEVP solution at an appropriately large value of $t_{0}$ to form the
optimized operators places some constraints on the scale of these late-time
pollutions. The optimized correlation function matrix,
$\matrixelement{0}{\Omega_{\mathfrak{m}}(t)\Omega_{\mathfrak{n}}^{\dagger}(0)}{0}$,
is diagonal at $t=t_{0}$, and close to diagonal for timeslices close to
$t_{0}$.161616In practice this matrix is not _exactly_ diagonal because the
optimized operators are constructed with the ensemble averaged vectors
$v_{\mathfrak{n}}$. The two-point current correlation functions,
$\matrixelement{0}{\mathcal{J}(t)\Omega^{\dagger}_{\mathfrak{n}}(0)}{0}$, are
calculated on the same ensemble, and with the same time sources as the matrix
of correlators used in the GEVP, leading to a significant correlation between
the contributions $\varepsilon_{\mathfrak{m},\mathfrak{n}}(t)$ and the off-
diagonal elements of the matrix of optimized correlators, suggesting that
$\varepsilon_{\mathfrak{m},\mathfrak{n}}(t)$ will be also suppressed for
timeslices around $t_{0}$ 171717Contributions from terms like Eq. 48 to
$R_{\mathfrak{n}}(t)$ are not relevant for $t<t_{0}$, where pollution from
states _higher_ in energy come to dominate correlation functions..

In the analysis reported in this paper, the only cases found where an excited
state energy level has a current matrix element that is significantly smaller
in magnitude than the current matrix element for a lighter state, and hence
where this ‘late-time’ pollution needs to be considered, are those energy
levels lying close to non-interacting $K\overline{K}$ energies and having
dominant overlap onto $K\overline{K}$-like operators. For these levels, based
upon the logic above, we restrict the time-windows for constant and constant-
plus-exponential fit forms to more modest values of $t_{\text{max}}\sim
1.5\,t_{0}$, thus excluding late-times where the pollution from lighter states
begins to become significant. In practice, we only need to impose this
fitting-window restriction explicitly for three states, the second excited
energy level of each of the irreps $[100]\,A_{1}$, $[110]\,A_{1}$, and
$[111]\,E_{2}$. Four other states lying close to non-interacting
$K\overline{K}$ energies in other irreps show a modest late-time enhancement
of the form suggested by Eq. 48, but the AIC-value quantifying the quality of
fits favored those with lower values of $t_{\text{max}}$ anyway.

We present below an example case supporting the arguments presented above, in
which we will reconstruct the time-dependence of correlation functions in
terms of contributions of the form of Eq. 48. For convenience of presentation
we will suppress kinematic factors not relevant to the illustration, by
introducing unit-normalized optimized operators
$\matrixelement{\mathfrak{n}}{\Omega^{\prime{\dagger}}_{\mathfrak{m}}(0)}{0}=\delta_{\mathfrak{n},\mathfrak{m}}$,
the ratio,

$R^{\prime}_{\mathfrak{n}}(t)\equiv\frac{\langle
0|\mathcal{J}(t)\,\Omega^{\prime{\dagger}}_{\mathfrak{n}}(0)|0\rangle}{\langle
0|\Omega^{\prime}_{\mathfrak{n}}(t)\,\Omega^{\prime{\dagger}}_{\mathfrak{n}}(0)|0\rangle}\,,$

and constants $\mathfrak{f}^{(L)}_{\mathfrak{n}}$ describing the leading
constant time-dependence of ${R}^{\prime}_{\mathfrak{n}}(t)$.

(a) Level 0. (b) Level 1. (c) Level 2. (d) Level 3.

Figure 23: Ratio of two-point correlators $R^{\prime}_{\mathfrak{n}}(t)$
(unimproved current) for discrete states in the $[100]\,A_{1}$ irrep compared
with the extracted constants $\mathfrak{f}^{(L)}_{\mathfrak{n}}$, the expected
contributions from lower energy levels,
$\widehat{\varepsilon}^{\prime}_{\mathfrak{m},\mathfrak{n}}(t)$, and the sum
of all contributions, $\widehat{R}^{\prime}_{\mathfrak{n}}(t)$.

Objects,

$\widehat{\varepsilon}^{\prime}_{\mathfrak{m},\mathfrak{n}}(t)\equiv\mathfrak{f}^{(L)}_{\mathfrak{m}}\frac{\matrixelement{0}{\Omega^{\prime}_{\mathfrak{m}}(t)\Omega_{\mathfrak{n}}^{\prime{\dagger}}(0)}{0}}{\matrixelement{0}{\Omega^{\prime}_{\mathfrak{n}}(t)\Omega_{\mathfrak{n}}^{\prime{\dagger}}(0)}{0}}\,,$
(49)

are analogous to Eq. 48 for the ratio $R^{\prime}_{\mathfrak{n}}(t)$, whenever
$E_{\mathfrak{n}}>E_{\mathfrak{m}}$. It follows that we should be able to
reconstruct as

$\widehat{R}^{\prime}_{\mathfrak{n}}(t)=\mathfrak{f}^{(L)}_{\mathfrak{n}}+\sum_{\mathfrak{m}<\mathfrak{n}}\widehat{\varepsilon}^{\prime}_{\mathfrak{m},\mathfrak{n}}(t)\,,$

and this will have the same time-dependence as $R^{\prime}_{\mathfrak{n}}(t)$
for times after the small-time excited state contributions have died out.

Comparisons between $R^{\prime}_{\mathfrak{n}}(t)$ and
$\widehat{R}^{\prime}_{\mathfrak{n}}(t)$ for levels in the $[100]\,A_{1}$
irrep are shown in Fig. 23. The values of $\mathfrak{f}^{(L)}_{\mathfrak{m}}$
for $\mathfrak{m}<\mathfrak{n}$ needed in Eq. 49 are obtained serially from
timeslice fits to the lower-lying energy level correlation functions. Level 2
meets our criteria for having non-negligible late-time pollution since its
$\mathfrak{f}^{(L)}$ value is significantly smaller than the corresponding
values for levels 0 and 1, with the value for level 1 being over ten times
larger in absolute value. We see that while the contribution of lower-lying
states to levels 1 and 3 is negligible, for level 2, the contribution of the
nominally suppressed level 1 is observably large. We also see that restricting
timeslice fits to values of $t_{\mathrm{max}}\lesssim 1.5t_{0}$ will reduce
the impact of this pollution considerably.

##  Appendix G Coupled-channel form factor parameterization

To describe timelike form factors in the coupled-channel region we use
polynomial parameterizations for the smooth functions in Eq. 5,

$\mathcal{F}_{a}(s)/m_{\pi}^{2}=\sum_{n=0}^{N_{a}}\,h_{a,n}\cdot\quantity(\frac{s-s_{0}}{s_{0}})^{n}\,,$
(50)

and we characterize each of them by the order of the polynomials, i.e.
$\\{N_{\pi\pi},N_{K\overline{K}}\\}$, as listed in Tab. 1.

As an example, the result of the fit to $\mathcal{F}_{\mathfrak{n}}^{(L)}$
using the _reference coupled-channel_ scattering parameterization, and
$N_{\pi\pi}=N_{K\overline{K}}=1$, has the following parameters: