Isoscalar $\pi\pi, K\overline{K}, \eta\eta$ scattering and the $\sigma, f_0, f_2$ mesons from QCD

We present the first lattice QCD study of coupled isoscalar $\pi\pi,K\overline{K},\eta\eta$ $S$- and $D$-wave scattering extracted from discrete finite-volume spectra computed on lattices which have a value of the quark mass corresponding to $m_\pi\sim391$ MeV. In the $J^P=0^+$ sector we find analogues of the experimental $\sigma$ and $f_0(980)$ states, where the $\sigma$ appears as a stable bound-state below $\pi\pi$ threshold, and, similar to what is seen in experiment, the $f_0(980)$ manifests itself as a dip in the $\pi\pi$ cross section in the vicinity of the $K\overline{K}$ threshold. For $J^P=2^+$ we find two states resembling the $f_2(1270)$ and $f_2'(1525)$, observed as narrow peaks, with the lighter state dominantly decaying to $\pi\pi$ and the heavier state to $K\overline{K}$. The presence of all these states is determined rigorously by finding the pole singularity content of scattering amplitudes, and their couplings to decay channels are established using the residues of the poles.


I. INTRODUCTION
The composition of the lightest hadrons with scalar (J P = 0 + ) quantum numbers remains a mystery. From one viewpoint, this is peculiar -being rather light and thus having only limited possible decay channels, we might expect them to provide a relatively simple system to study, but yet they have proven to be rather challenging. The lightest isoscalar scalar mesons are particularly interesting, as they appear to be two quite different objects -an extremely broad and light f 0 (500) (historically called the σ, a name we will continue to use), and a rather narrow f 0 (980) appearing very close to the KK threshold [1,2]. These states appear in ππ scattering in a way which challenges our naive view of hadron resonances as peak-like enhancements in cross-sections, appearing rather as a very broad bump with a sharp dip at the KK threshold [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].
The σ and f 0 (980) have commonly been placed together with the K 0 (800) (or κ) and a 0 (980) resonances into a nonet of (broken) SU (3) flavor. The κ appears to be a strange analogue of the σ, being a very broad resonance enhancing πK scattering at low energies, while the isovector a 0 (980) closely resembles the f 0 (980). The a 0 (980) appears at the KK threshold, and likely shares the f 0 (980)'s strong coupling to KK. It is observed through its decay to πη, but the lack of direct data on πη elastic scattering limits the precision with which it can be studied.
In the simplest interpretation of the constituent quark model, the scalar nonet would arise from qq( 3 P 0 ) constructions, but such an assignment for the states discussed * briceno@jlab.org † dudek@jlab.org ‡ edwards@jlab.org § djwilson@maths.tcd.ie above appears unnatural given that the related qq( 3 P 1,2 ) states lie at much higher masses. In addition, this picture provides no explanation for the near degeneracy of the the f 0 (980) and a 0 (980), nor why the σ and κ are apparently lighter than them. One suggestion which attempts to remedy this is to have the nonet be of dominantly qqqq construction [18], with the a 0 (980) and f 0 (980) having hidden strangeness, making them heavier than the σ which is proposed to contain only light quarks. Alternatively, given their proximity to the KK threshold, a natural suggestion is that the a 0 (980) and f 0 (980) might be dominated by KK molecular configurations, where the binding is provided by residual interhadron forces [19] (see Ref. [20] for a review of resonances as hadronic molecules). Similarly, dynamical modeling which has a qq 'seed' being dressed by strong coupling to its meson-meson decay modes has proven capable of generating light scalar resonances which resemble those in experiment [21][22][23][24]. Interpretation of the scalar mesons can be further complicated if glueball basis states are assumed to also play a role in the isoscalar case [25]. In contrast to the scalar sector, the lightest isoscalar tensor mesons (J P = 2 + ) present a situation that much more closely aligns with our naïve view of hadron resonances. In the constituent quark picture, two isoscalar states can be constructed from combinations of uu + dd and ss pairs in the 3 P 2 configuration. Experimentally the f 2 (1270) and f 2 (1525) appear as "bump-like" enhancements which are well described by Breit-Wigner parametrizations, with the lighter state dominantly decaying to the ππ final state (84%), and the heavier state to KK (89%) [26]. These decay characteristics are taken as support for the quarkmodel assignment f 2 (1270) ∼ uū + dd, f 2 (1525) ∼ ss using the phenomenology of the 'OZI' rule, where qq pair creation is proposed to dominate existing qq annihilation in hadron decays. The rates of two-photon decays of these states have also been presented as support for these assignments (see Ref. [27] and references therein).
The presence of a resonance of angular momentum, J, in a hadron scattering process has a rigorous signature in the form of a pole singularity in the partial-wave amplitude, t J (s). A pole at a complex value of Mandelstam s = s R = m R − i 1 2 Γ R 2 influences scattering for real values of sin particular for small value of the width, Γ R , the effect for an isolated resonance is typically a narrow peak. In general, the scattering amplitude is a matrix in the space of kinematically open channels, and near a resonance pole the elements of the matrix take the form t ij ∼ ci cj s R −s , where c i is the coupling to the i channel. A rigorous presentation of hadron resonances will take the form of a list of pole positions and the associated couplings.
As described above, dynamical models working at the level of constituent quarks and/or hadrons are informative, and offer frameworks within which experimental observations may be placed, but ultimately all observed hadron phenomena must have an origin within QCD. Lattice QCD provides an explicit numerical approach to studying QCD at energies relevant to resonance physics, where it is fundamentally non-perturbative. After discretizing the theory on a finite hypercubic grid, an ensemble of gluon field configurations can be generated, and using these, correlation functions evaluated. By computing appropriate two-point correlation functions, the discrete spectrum of QCD eigenstates in the finite-volume defined by the lattice can be determined. These can be related to scattering amplitudes through the Lüscher formalism [28][29][30] and its extensions [31][32][33][34][35][36][37][38] which connect the volume dependence of the discrete spectrum to scattering amplitudes in an infinite volume.
An approach which has proven successful [39][40][41][42] proceeds by parametrizing the energy-dependence of coupledchannel amplitudes and fitting a large set of energy levels, from one or more lattice volumes, within a kinematic window [35]. A dense spectrum of energy levels will tightly constrain the possible energy-dependence of the scattering t-matrix, and to acquire as many energy levels as possible, we may consider systems with various total momenta. Currently this approach is limited to energies below threehadron and higher thresholds -to go beyond this an extension of this formalism is required, and progress in this direction is being made [43][44][45][46][47][48].
In this first study of excited isoscalar meson resonances, we will calculate a version of QCD featuring light quarks which are heavier than those found experimentally, leading to a pion mass of 391 MeV. In this world, three-and fourmeson thresholds appear at higher energies, providing a larger energy region in which we can perform rigorous extraction of scattering amplitudes. At this pion mass we have previously studied the other symmetry channels that would make up a flavor nonet: I = 1/2, S = ±1 through coupled πK, ηK scattering [39,40], and I = 1, S = 0 through coupled πη, KK [42].
In Refs. [39,40], the scattering matrix in the πK, ηK sector was found to be almost completely decoupled, with ηK → ηK scattering being rather weak, while the πK S-wave was found to be attractive at threshold, and to feature no rapid energy dependence. This was interpreted as being due to a virtual bound-state singularity over 200 MeV below πK threshold, as well as a much heavier broad scalar resonance pole above 1.3 GeV lying far into the complex plane.
In the isovector channel, in Ref. [42], the scattering matrix for πη, KK was found to feature strong channel coupling in S-wave, with very rapid energy dependence around the KK threshold. This was shown to be due to a resonance pole singularity lying close to the KK threshold. The couplings of this resonance to the πη and KK channels were found to be of comparable size.
The isoscalar sector that we now turn to is notoriously more challenging for lattice QCD, owing to the need to compute completely disconnected diagrams in which all quarks and antiquarks annihilate. Using the distillation framework, we have been able to evaluate all required contributions to correlation functions with good statistical precision. The elastic ππ → ππ part of I = 0, S = 0 scattering (below KK threshold) was presented in Ref. [49], where a bound-state pole was found about 24 MeV below the ππ threshold. The same reference also obtained this amplitude using a lighter value of the quark masses corresponding to m π ∼ 236 MeV, and found that the bound state evolved into a broad resonance, closely resembling the experimental σ.
In calculations with m π ∼ 391 MeV, J P = 2 + resonances have also been determined that are found to be narrow. The lowest-lying K 2 resonance in πK is found to have a mass m = 1576(7) MeV and width Γ = 62(12) MeV decaying into πK only. The isovector a 2 was found at m = 1506(4) MeV , Γ = 20(3) MeV with approximately equal couplings to πη and KK . These extractions were somewhat less rigorous than those of the scalar mesons owing to the neglect of possible three-meson decays.
In this paper we will extend the work presented in Ref. [49] at m π ∼ 391 MeV to consider also the energy region above KK threshold, and study isoscalar coupled ππ, KK, ηη scattering in S-wave and D-wave. We will confirm the σ bound-state previously observed, and furthermore identify an f 0 (980)-like state in S-wave, appearing close to KK threshold, with strong couplings to both ππ and KK . In the D-wave scattering amplitudes, we observe two clear peaks in the scattering amplitudes, which are interpreted as being due to two resonances, which resemble the experimental f 2 (1270) and f 2 (1525), being relatively narrow, and with the lighter state dominantly coupling to ππ and the heavier to KK .
The remainder of the paper is structured as follows: Section II presents the finite-volume spectra determined in explicit lattice QCD calculations. Section III briefly describes the technique for determining scattering amplitudes from finite-volume spectra before presenting results for elastic ππ scattering (III A), coupled S-wave ππ, KK, ηη scattering (III B) and coupled D-wave scattering (III C). Section IV examines the pole singularity content of the determined amplitudes, which is interpreted in terms of resonances in Section V. In Section VI, we consider the complete flavor nonet of scalar mesons determined at m π ∼ 391 MeV, before summarizing the current calculation in Section VII.
For a system at rest, spectra are computed according to irreducible representations (irreps) of the cubic group, which contain subductions of the angular momenta which characterize the infinite-volume spectrum. In systems of identical meson pairs with definite G-parity, only even partial waves contribute for even isospin, which subduce according to Table II of [59]. In particular, we are interested in the irreps A + 1 which contains subductions of J P = 0 + , 4 + . . ., and E + , T + 2 containing J P = 2 + , 4 + . . .. For systems with non-zero total momentum, the symmetry is reduced to the little group of the cubic group, which is defined by the allowed rotations of a cube which leave the total momentum invariant. In this work, we consider systems with total momentum up to [200].
Stable meson masses on these lattices can be found in [42] -particularly relevant here are pseudoscalar meson masses: a t m π = 0.06906(13), a t m K = 0.09698(9) and a t m η = 0.10364 (19). The anisotropy determined from stable meson dispersion relations is ξ = a s /a t = 3.444 (6).
The isoscalar sector is notoriously challenging to study within lattice QCD. In previous studies, various approximations have been made in order to make calculations practical [67], including the omission of disconnected diagrams or the use of only a small basis of operators. These approximations cannot in general be justified. In this paper we will compute matrices of correlation functions in a large basis of operators 2 , including many which resemble a pair of mesons each having definite momentum, with the correlator construction achieved using the distillation framework [68]. All required Wick contractions, including those in which quarks at the source or sink annihilate, are included without further approximation.
Matrices of correlation functions are analyzed variationally by solving a generalized eigenvalue problem [69,70] to yield discrete spectra of states. The operators resembling pairs of mesons with definite momentum are themselves constructed using variationally-optimized π, K and η operators [59]. The KK operators are constructed with definite G-parity, and the ηη operators contain both light and strange components. Our procedures for correlator construction and variational analysis have been described in previous papers [51,56,59,68], and we will not repeat details here, except to point out that we use a 'weightingshifting' correction [59] to reduce the (small) pollution due to the finite temporal extent of the lattice. For restframe irreps, this amounts to computing the difference between two timeslices, which also acts to remove the time-independent vacuum contribution in the [000]A + 1 irrep. Figure 1 shows spectra determined in three lattice volumes for all A 1 irreps with total momentum up to [200]. The error bars include, as well as the statistical error due to the use of a finite sample of gauge-field configurations, also an estimate of systematic error obtained by varying the details of the variational analysis, e.g. reasonable variation of t 0 , the extent of fitting-time windows, and the content of the variational operator basis. Significant shifts are observed with respect to non-interacting meson-meson energies shown by the curves, indicating the presence of strong scattering dynamics. The levels shown in black and blue will be used later to constrain three-channel scattering in S-wave. We utilise 57 energy levels densely filling an energy region from below ππ threshold to some way above ηη threshold at a t E cm ∼ 0.24.
Energy levels below KK threshold computed on these lattices previously appeared in [49], where they were used to determine the ππ elastic scattering amplitude which was found to feature a bound-state pole identified with the σ. The spectra in Figure 1 contain some small differences with respect to those presented in [49], typically at the level of statistical fluctuations. For this paper we have undertaken a thorough consideration of the variations seen in correlator extraction, and worked with somewhat larger operator bases in order to access the coupled-channel energy region. Figure 2 shows the [000]A + 1 irrep, and alongside each energy level we plot a histogram showing the relative contribution to each state from each operator in the variational basis. We highlight 5 distinct types of operator: ππ (red), uΓu +dΓd (grey),sΓs (light green), KK (green), and ηη (blue), and the spectrum is seen to not be diagonal in this basis. A level below ππ is seen on all three volumes dominated by overlap with ππ andūΓu +dΓd operators, which is connected to the σ meson that appears as a shallow bound state on these lattices [49] 3 . On all There are no three-meson or higher thresholds in this energy region. Discrete energies determined in variation analysis of correlator matrices are plotted, and the uncertainties include systematic variation as described in the text. Energy levels colored blue are those which have large overlap with ηη-like operators. Ghosted points are not used in the determination of coupled-channel scattering amplitudes that will be presented in Section III B.
three volumes we see there is a state coincident with ηη threshold that has dominant overlap with an operator resembling η[000]η[000]. The statistically negligible shift with respect to the non-interacting level, and relatively small mixing between the ηη operator and the other operators may be a hint that the ηη channel is not as strongly interacting or as strongly coupled as ππ and KK. While operator overlaps are useful for building intuition, they are not a rigorous tool, so we reserve further comment until after the scattering amplitudes have been extracted and analyzed in Section III B. Figure 3 shows the spectra determined in three lattice volumes for irreps which feature J P = 2 + scattering as the lowest angular momentum. We extract 34 energy levels shown in black and blue that are to be used in the extraction of scattering amplitudes in Section III C. Figure 4 shows the spectrum and associated histograms for the E + and T + 2 irreps, where the pattern of operator overlaps is seen to be quite different to that in A + 1 , with notably less 'mixing' between the light and strange sectors. Statistically significant shifts from the non-interacting energies are observed, and there are more levels than would be expected based on counting the non-interacting energies on each volume in this energy region, hinting that narrow states could be present. The histograms go further, suggesting that there are likely to be two resonances, one dominated by light quarks, and another dominated by strange quarks. This is perhaps clearest in  non-interacting levels, but we still see two levels, one near a t E cm = 0.26 dominantly overlapping withūΓu +dΓd constructions and a second near a t E cm = 0.28 dominantly overlapping withsΓs constructions. Again, we reserve further comment until after the amplitudes have been extracted in Section III C. We now turn to extracting the scattering amplitudes from these spectra, beginning with the S-wave elastic region before proceeding to investigate coupled-channel S-and D-wave amplitudes.

III. DETERMINING COUPLED-CHANNEL SCATTERING AMPLITUDES
It has been shown that the discrete spectrum of states in a periodic finite spatial volume is determined by the energy-dependence of infinite-volume hadron-hadron scattering amplitudes [28][29][30][31][32][33][34][35][36][37][38]. The particular form of the relationship relevant to coupled-channel pseudoscalarpseudoscalar scattering described by a t-matrix, t(E cm ), can be compactly expressed as [35] with ρ(E cm ) a diagonal matrix of phase-space factors, finite-volume matrix M is diagonal in channel space, but it is in general not diagonal in angular momentum.
This finite-volume formalism was recently reviewed in some detail in Ref. [71], and more description of our implementation (using the same notation as we will use in this paper) is given in Ref. [42]. Our approach is to parameterize the energy-dependence of the t-matrix, solving Eq. 1 for the discrete spectrum with a given choice of parameter values. This spectrum is then compared to the lattice spectra shown in e.g. Figure 1 in a correlated χ 2 function. Minimizing the χ 2 by varying parameter values we obtain a best estimate for the t-matrix. Subtleties associated with broken rotational symmetry due to the cubic boundary, both at rest and in-flight, expressed mathematically by the subduction of partial waves into irreps of the relevant symmetry group, are discussed in Refs. [42] and [56].
One approach to ensure that the t-matrix in partial-wave satisfies multichannel unitarity (required in order for Eq. 1 to have solutions) is to express it in terms of a K-matrix, where the factors (2k i ) − provide the required kinematic behavior at threshold for the -wave. The elements K ij (E cm ) form a symmetric matrix that is real for real values of E cm , and the elements I ij form a diagonal matrix whose imaginary part is fixed above threshold by unitarity to be −ρ i . An option for the real part of I which ensures the amplitude behaves sensibly below threshold and for complex values of the energy is to use the Chew-Mandelstam phase-space -a discussion can be found in Ref. [40]. We will make use of a range of parameterizations for the matrix K in the S-and D-waves, considering mostly the coupled three-channel problem ππ, KK, ηη. We will later also have cause to also consider an application of the Jost functions, through which we may explicitly control the singularity content of the amplitudes, at the cost of losing guaranteed unitarity for all parameter values. The possibility of extracting information about scalar mesons from lattice QCD spectra has previously been discussed in the context of various amplitude parameterizations in Refs. [72][73][74][75].
A. Elastic S-wave ππ scattering An analysis of the discrete spectrum of states in the energy region below KK in terms of elastic ππ scattering in S-wave was previously presented in [49]. In this case, under the (verified) assumption that higher partial-waves make a negligible contribution to Eq. 1, the problem reduces to a one-to-one mapping from the computed E cm value to a value of the scattering amplitude at that energy. In this way the elastic scattering amplitude was mapped out, and a behavior compatible with a bound-state that we associated with the σ was observed.
For this paper we have considered larger matrices of correlation functions, including as well as those used in Ref. [49], also a significant number of KK and ηη operators, and the resulting variational analysis leads to spectra that differ slightly from those presented in [49], in particular through better estimation of a systematic error. We use these improved energy levels here to determine the elastic scattering amplitude, expressed via the phaseshift (t ππ,ππ = ρ −1 ππ e iδππ sin δ ππ ), assuming that = 2 and higher partial waves have a negligible impact on the finite-volume spectrum 4 . In Figure 5 we present discrete phase-shift points, as well as effective-range-expansion descriptions of the finite-volume spectrum below KK threshold.
Re-expressing the elastic t-matrix as t = Ecm 2 1 k cot δ−ik makes it clear that if a graph of k cot δ intersects the curve ik at energies below threshold, where ik = −|k|, there will be a bound-state pole singularity. Such a crossing is clearly visible in the lower pane of Figure 5, corresponding to a bound-state σ. We will explicitly address the singularity content of scattering amplitudes and their interpretation in terms of bound-states and resonances in Section IV.
The energy dependence of the elastic amplitude can be described within the S-wave effective range expansion, k cot δ 0 = 1 a + 1 2 rk 2 + n=2 P n k 2n . Describing those energies in the region 0.12 ≤ a t E cm ≤ 0.175, by a scattering length plus effective range form we obtain a/a t = −36.6 ± 3.2, r/a t = −30.9 ± 5.2 with a χ 2 /N dof = 16.3/(17 − 2) = 1.09. A larger energy region, 0.1 ≤ a t E cm ≤ 0.185, can be described if we include also a k 4 term in the expansion, yielding a/a t = −43.1 ± 3.9, r/a t = −18.1 ± 4.4 with a χ 2 /N dof = 18.1/(23 − 2) = 0.91. Both descriptions are displayed in Figure 5. In this elastic analysis we avoid considering levels lying just below the KK threshold, as in finite-volume they are impacted by the t ππ→KK and t KK→KK elements of the t-matrix [35,71].
B. Coupled S-wave ππ, KK, ηη scattering A naïve examination of the energy levels presented in Figure 1, where an 'extra' level is present below ππ threshold, strongly suggests the existence of the bound-state σ discussed in the previous section, but such simple analysis does not obviously indicate any narrow resonances at higher energy -while significant shifts of energy levels away from non-interacting energies are observed, there are not clearly any 'extra' levels that one might associate with a nearly-stable state.
We will proceed by attempting to describe 57 energy levels lying below a t E cm = 0.24 shown in black and blue in Figure 1 using a range of three-channel (ππ, KK, ηη) S-wave amplitude forms in Eq. 1. In this section we will make use of K-matrix parameterizations as described previously.
An illustrative example of a successful description of the spectrum is provided by a K-matrix in which the matrix inverse of K is parameterized as where the row (and column) channel ordering is ππ, KK, ηη. There are eight free parameters, namely the constants a . . . h, and the Chew-Mandelstam phasespace used for I is subtracted for each channel at the relevant threshold. The best description of the spectra with this amplitude 5 , shown in Figure 6, has This is clearly a very successful description of the lattice QCD energy levels -the fact that this reduced χ 2 (and that of many other successful descriptions) is slightly below 1 may indicate that we have been a little too conservative in our estimation of systematic error on the energy levels. We note that the fit to 57 energy levels shown in bold also provides a quite reasonable description of levels at slightly higher energies, and on the L/a s = 16 lattice in-flight (shown ghosted), that we conservatively did not include in the fit.
The energy dependence of the resulting amplitude is shown in Figure 7, where what is plotted is a quantity that is proportional to the cross-section for various scattering processes. Above a broad bump in ππ → ππ at low energies caused by the tail of the σ bound-state, is a sharp dip just below the opening of the KK threshold. Along with a slight cusp in ππ → ππ at KK threshold, we observe a rapid turn on of the KK channel. There is clearly very little activity in ηη, and only tiny kinks in the ππ and KK channels at ηη threshold. The small circles shown alongside the energy axis indicate the 57 energy levels used to constrain the amplitude -we note that they densely span the entire energy region, overconstraining the energy dependence of the t-matrix. The statistical uncertainties on the π, K and η masses prove to have a negligible effect on the amplitude determination.
In passing we note that the rest-frame spectrum in the three volumes considered alone would have provided us with only 15 energy levels. Of these only 9 are in the coupled-channel region and only three of these have large ηη overlap. Such limited information would have provided minimal constraint on our parameterization and consequently would have given us little confidence in our final result.
Another way to present the energy dependence of the scattering amplitudes is by plotting the magnitudes and phases of the diagonal elements of the S-matrix 6 The two-channel unitarity constraint is such that between the KK and ηη thresholds, S ππ,ππ = S KK,KK and this magnitude can be associated with what is usually called the (in)elasticity, η, while the phases ϕ ππ,ππ , ϕ KK,KK are identified with the channel phase-shifts, δ ππ , δ KK . Above the ηη threshold, three-channel unitarity comes into play, and such FIG. 6. A1 irrep spectra as described by the amplitude given in Eq. 3. Black/blue points show the lattice QCD spectrum of Figure 1 and orange curves and points (displaced slightly for clarity) the result of the fitted amplitude. simple identifications can no longer be made, although the relatively weak coupling to the ηη channel suggests that a decoupled "2 + 1" channel (coupled ππ, KK plus decoupled ηη) description might be a reasonable approximation. Figure 8 shows these phases and magnitudes. scattering matrix is provided by the Argand diagram, which is shown for ππ → ππ in Figure 9. Starting at ππ threshold, the amplitude initially moves rapidly clockwise along the unitarity circle under the influence of the σ bound-state, before doubling back and more slowly passing though the point t ππ,ππ = 0 corresponding to the dip in Figure 7. Shortly after this, upon reaching the KK threshold, the amplitude moves inside the unitarity circle as probability is lost from ππ into the KK channel. The opening of the ηη channel is marked by only a small kink in the curve. Examining Figures 7, 8, 9, it is clear that we are dealing with a rather non-trivial scattering system whose resonant content cannot immediately be inferred. Certainly the amplitude looks nothing like the traditional view of a simple resonance appearing as a clearly defined bump on a slowly varying background. However, it is important to note that by fully respecting unitarity, significant constraints have been placed on the possible energy dependence of the amplitude, particularly in the elastic region -this, coupled with the presence of the σ bound-state, could cause the appearance of a resonance to be significantly distorted with respect to our naïve intuition.
A rigorous definition of the resonant content of a scattering system will come from examining the t-matrix at complex values of s where pole singularities will appear if the system contains resonances. We will explore the pole content in a later section, but at this stage it is worth noting that the sharp dip in ππ and rapid turn on of KK intensity suggests there may be a singularity in the vicinity of the KK threshold. We also note in passing that although the σ at this quark mass appears as a bound-state rather than the broad resonance seen in experiment, the energy dependence observed in Figure 7 is not dissimilar to that extracted from pion beam experiments [3][4][5]7] as shown in, for example, Fig. 1 of Ref. [76].

Varying the amplitude parameterization
In the previous section we discussed one particular K-matrix amplitude parameterization that successfully described the lattice QCD spectra of Figure 1. We have explored a wide variety of parameterizations, and in this section we will report on the variation observed in the resulting amplitudes. For clarity of presentation we have selected an illustrative set of 20 parameterizations (including the one presented in the previous section), all of which have χ 2 /N dof below 1.05, and which show no excessive parameter correlation 7 . Many other parameterizations were found to successfully describe the lattice 7 we also reject any amplitude which we find to feature a nearby off-axis pole singularity on the physical sheet -such behavior is acausal and reflects the lack of explicit analyticity constraints on our K-matrix amplitudes. QCD spectra with behavior compatible with one or more of the 20 shown here 8 . Figure 10 shows the variation in central values for 19 parameterizations along with the reference amplitude described in the previous section, and what is clear is that the elastic region behavior, including the position of the dip, shows very little variation, and neither does the behavior of amplitudes in the region between the KK threshold and the ηη threshold. The bulk of variation under change of amplitude form lies at and above the ηη threshold, and is ultimately associated with how strongly the fit allows the ηη channel to couple into the strongly coupled ππ, KK system. The variation is, however, only modest, with the qualitative behavior of the amplitudes being the same for all these successful fits, and in particular the rapid energy dependence around the KK threshold is rather well pinned down.
The in-flight A 1 irreps used to constrain the amplitudes do receive contributions from = 2 amplitudes, although we expect these to be significantly suppressed at these low energies by the centrifugal barrier. We have explicitly estimated their effect by including in our fits the = 2 amplitudes which successfully describe the spectra presented in Figure 3, to be discussed in the next section. When these amplitudes were included in Eq. 1 no significant change in the S-wave amplitudes was observed and we conclude that D-wave amplitudes play a negligible role in determining the spectra presented in Figure 1.
One of the amplitudes included in our illustrative set of 20 includes an Adler zero, implemented by multiplying a parameterization of K(s) by a factor (s − s A ) with the position of the zero set at the value suggested by leading-order chiral perturbation theory, s A = 1 2 m 2 π . The resulting fitted amplitude does not show any noteworthy differences with respect to the others considered. We explored the dependence on the position of the zero, by varying s A in the region from −5m 2 π up to 3 2 m 2 π , allowing the other parameters in the amplitude to float freely. We observed that the χ 2 was largely insensitive to the position of the zero, with a very slight preference for large negative values (where its presence becomes of decreasing importance). The appearance of the amplitudes hardly varies under change in s A , and the pole singularity content of the amplitudes (to be discussed in Section IV) is similarly insensitive. We conclude that an Adler zero is not an important feature of the S-wave t-matrix at m π ∼ 391 MeV.  Figure 3, assuming three-channel scattering, ππ, KK, ηη, and ignoring contributions from ≥ 4. We note that only four levels have significant overlap with ηη-like operator constructions, which suggests that we will have only limited constraint on the ηη sector. However, the proximity of these four levels to the corresponding noninteracting ηη levels suggests that the channel may well be only weakly interacting, and that solutions in which ηη is largely decoupled from ππ, KK may be successful.
In Figure 3 we observe the opening of the ηη channel within the energy region under consideration. We considered the difference in the spectrum extracted including and not including ηη operators and observed no statistically significant change, beyond addition of poorly determined energy levels at higher energies. We will neglect the effect of this channel.
In Figure 3 we also see that a three-hadron channel, ηππ, and a four-hadron channel, ππππ, open in the energy region being considered. We did not include any operators resembling these channels when we computed the correlation matrices. A complete formalism for relating finitevolume spectra to amplitudes including three-hadron and higher multiplicity scattering does not yet exist, but progress in that direction is being made [43][44][45][46][47][48]. We do have reason to believe that these channels will not have a significant impact at the energies we are considering. Experimentally, in virtually all three-hadron and higher multiplicity processes, the final state is dominated by the parts of the phase-space where two (or more) hadrons resonate through an isobar. In this case, the lowest-lying contributing isobar systems for ηππ with J = 2 would be a 2 π and f 2 π, and in [42], the a 2 decaying to πη was found with a mass a t m ∼ 0.26. Later in this paper we will find the lightest f 2 resonance at a similar mass. Thus we might estimate the energies at which we expect ηππ to become a considerable influence to be above a t E cm ∼ 0.33. A possible contribution from the isobar process ρρ in ππππ would be expected above a t E cm ∼ 0.30, and this in part motivates our decision not to consider in this first analysis any levels above a t E cm ∼ 0.30.
As was discussed in Section II, spectra and overlaps suggest there may be two narrow enhancements near a t E cm = 0.26, 0.28, and an efficient way to allow for two narrow resonances within a multichannel K-matrix is to use a parameterization which includes two real poles. When 'dressed' by the phase-space, real poles in the K-matrix with relatively small residues can give rise to an amplitude which resembles a Flatté form 9 . An illustrative example of a parameterization of the type given by Eqn. 2, which proves to be successful in describing the spectra in Figure 3 is provided by where in addition no real part is included in I (conventional phase-space). This amplitude features 9 real parameters: a mass and three channel couplings for each real pole plus a constant to allow for some gentle ηη energy dependence. This form allows the ηη channel to decouple if the fit selects small values for the relevant pole couplings (g (1,2) ηη ). The description of the lattice QCD spectra provided by this amplitude, with a χ 2 /N dof = 28.9/(34 − 9) = 1.15, is shown in Figure 11. The description of the input spectrum is seen to be good, and furthermore the fit makes predictions for states which agree rather well with computed levels that were conservatively not included in the fit (ghosted in the figure). Figure 12 shows the energy dependence of the resulting amplitudes, where we clearly observe significant enhancement in ππ slightly above a t E cm = 0.26 and in KK slightly above a t E cm = 0.28. Essentially no activity is seen in ηη as one would expect given the proximity of the relevant energy levels to the non-interacting ηη energies. The behavior observed in the ππ, KK sector would appear to be that of two narrow resonances, the lighter decaying strongly to ππ and more weakly to KK and the heavier likely decaying mainly to KK. We note in passing the inactivity of the various components of the = 2 amplitude in the kinematic region where the = 0 were predominantly constrained in the previous section, namely energies satisfying a t E cm < 0.24. This explains why the contamination from the = 2 partial wave into the A 1 irreps played a negligible role in the analysis of the spectrum.
It is interesting to examine the Argand diagrams for ππ → ππ and KK → KK for this t-matrix -they are presented in Figure 13. The lower-energy bump has what appears to be a canonical resonance behavior in ππ → ππ, with the amplitude going counterclockwise through the vertical at a t E cm = 0.264 and where the curve lying inside the unitarity circle is indicative of a loss of probability into the KK final state. The lower bump is visible in KK → KK as a strong kink at a t E cm = 0.264, and above this energy the amplitude goes back onto the unitarity circle, suggesting an approximate decoupling from ππ at these energies. Another canonical resonance behavior near a t E cm = 0.284 corresponding to the higher energy bump is present in KK → KK.
Ultimately a rigorous approach to the resonance content and their relative decay strengths to open channels will come through consideration of the pole singularities of the t-matrix -this will be discussed in Section IV B.

Varying the amplitude parameterization
The parameterization of Eq. 4 successfully describes the finite-volume spectra shown in Figure 3 in terms of two narrow bump structures. We find that only amplitudes featuring two such bumps are able to describe the spectrum, and in this section we present the result of considering a broader set of parameterizations, retaining two real poles in K, but adjusting some features: whether we allow the poles to couple to ηη, the form of the quantity added to the poles, and the nature of I (Chew-Mandelstam versus naive phase-space). In total 16 amplitudes are presented 10 , all of which have χ 2 /N dof < 1.2, and which do not feature overly large parameter correlations -the resulting amplitudes are shown in Figure 14. There is clearly very little variation under changes in the parameterization. We do observe a slight variation in the magnitude of the 'shoulder' in ππ → ππ at a t E cm ∼ 0.28 and in the strength of ππ → ππ above the energy region where we have constrained the amplitudes (a t E cm > 0.30). One particular amplitude parameterization finds a somewhat larger ηη → KK component (the lone visible green line in the lower pane), but this behavior begins outside the energy region where we have constrained the amplitudes.
In summary it appears that the ππ, KK sector is rather well determined, in particular the position of the 'bumps'. The ηη channel is less well constrained, but there are very strong hints that it is largely decoupled. In Section IV B we will examine the pole singularity content of these amplitudes and propose a resonance interpretation.

IV. RESONANCE POLES
The partial-wave amplitudes we have determined, like those extracted from experimental measurements, are evaluated at real values of the scattering energy, but an understanding of features of these amplitudes, such as peaks and cusps, comes from considering their singularity structure in the complex s-plane. As well as the branch points required by unitarity at the opening of each new channel, amplitudes can have pole singularities that we may identify with bound-states (if they lie on the real axis below threshold) and resonances (if they lie off the real axis). In the region of a pole singularity at s = s 0 , the elements of the t-matrix behave like and the real and imaginary parts of the pole position are often identified with the mass and width of a resonance as √ s 0 = m R ± i 2 Γ R . The residue at the pole can be factorized into complex-valued couplings, c i , which indicate how strongly the resonance couples to scattering channels.
The presence of branch cuts associated with each new channel opening means that the complex s-plane is multisheeted. In single channel scattering there are two sheets which can be labelled by the sign of the imaginary part of the cm-frame scattering momentum, k. Im k > 0 is called the physical sheet, since it includes the region s = E 2 cm +i lying just above the real axis, where physical scattering occurs. Moving down from the real axis into the complex plane for any energy above threshold, we pass though the cut onto the unphysical sheet where Im k < 0. Resonance poles appear in complex conjugate pairs on the unphysical sheet, but usually only the pole in the lower half-plane is close to physical scattering. Complex poles cannot appear on the physical sheet as their presence would indicate a violation of causality.
In the case of multichannel scattering, the sheet structure becomes more complicated. For n chan channels there is still a physical sheet on which all channel momenta have a positive imaginary part, but there are now 2 n chan − 1 unphysical sheets. We will attempt to avoid confusion in nomenclature by labelling sheets by the sign of the imaginary part of the momentum in the three channels, ordered as (ππ, KK, ηη). For physical scattering between the ππ and KK thresholds, sheet (−, +, +), which we will also call sheet II, is the nearest sheet to the scattering axis. Between the KK and ηη thresholds, sheet (−, −, +), which we will also call sheet III, is closest, and above ηη threshold, it is sheet (−, −, −) which is closest. Some further discussion of sheet structure and pole positions in the context of two-channel scattering can be found in [42]. A relevant observation is that a single resonance typically appears as a pole 11 on several unphysical sheets ('mirror poles'), with shifts in position that are small if the resonance is dominantly coupled to only one scattering channel, but which can be large if coupled strongly to multiple channels 12 .
Because the amplitudes we considered in Section III are described by explicit forms, we can continue them to complex values of s, and search for pole singularities on all Riemann sheets. Uncertainties on the pole positions and couplings extracted from the residues can be estimated by propagating through the correlated uncertainties on the fitted amplitude parameters.

σ bound-state pole
Our discussion in Section III A concerning elastic ππ scattering leads us to expect that there is a bound-state pole singularity, lying below ππ threshold, that we may associate with a stable σ meson. Indeed in all successful descriptions of the lattice QCD spectrum, including all those presented in Section III B, we find the amplitude features such a pole singularity, and we observe very little movement in its position with variation of parameterization form. A best estimate for its position, including in the uncertainty a conservative measure of the rather small degree of variation with parameterization form is a t √ s 0 = 0.1316 (9). The main observed variation is a systematically lower mass (typically by roughly 0.0005) for parameterizations which use the ordinary phase-space for I rather than the Chew-Mandelstam form. The coupling to the ππ channel, a t c ππ = 0.092(4), shows very little variation under parameterization form, with the uncertainty dominated by the uncertainty on the energy levels.
In multichannel fits, couplings of this state to the KK and ηη channels can be obtained from t-matrix residues at the pole, but the very large extrapolation below the relevant thresholds renders these numbers largely meaningless. This pole position and coupling differs slightly from that presented in [49] -as discussed earlier, it follows from fitting a set of energy levels that are not identical to those presented in that reference. Owing to the more cautious estimation of systematic errors in this paper, the 11 actually a complex-conjugate pair of poles, but we typically speak of this as one pole. 12 an illustrative presentation of this effect in the context of the Flatté amplitude for two-channels will be given later in Section VI.
result above should be considered to supersede the one in [49].

f0 resonance pole
In Section III B we presented suspicions that the sharp dip in ππ → ππ and rapid turn on of KK amplitudes at the KK threshold might be due to a nearby resonancethis can be tested by analytically continuing the t-matrix into the complex plane, and searching for poles on all Riemann sheets. The consistent result of doing this for all successful parameterizations is the presence of a nearby pole on sheet II(−, +, +). This pole is partnered by a 'mirror' pole on sheet (−, +, −), which we expect is of the same origin and appears at roughly the same energy due to the small coupling to ηη. As an example, the amplitude described in Section III B, given by Eq. 3, is found to have poles at and while there are also poles found on sheets III(−, −, +) and (−, −, −), they are rather far into the complex plane such that they will not have a significant impact on physical scattering 13 .
In [42] we discussed the possible manifestation of a sheet II pole lying near the second threshold in the scattering amplitude for two strongly-coupled channels, which is almost exactly the scenario we are witnessing here. There we concluded that the scattering amplitude of the primary channel would have an asymmetric peak near the second threshold, which obviously differs from the 'dip' behavior we see in the ππ amplitude in, for example, Fig. 7. This dip is explained by the interference between the σ and the f 0 poles, something that was not considered in the simple discussion presented in [42]. Figure 15 shows the variation in pole location in the complex s-plane (upper panel) and complex k KK -plane (lower panel) for the 20 amplitudes presented in Section III B, where we show only the poles on sheets I(+, +, +), II(−, +, +), III(−, −, +) and IV(+, −, +). All these amplitudes feature a sheet II pole at roughly the same location, while a few also have sheet III or sheet IV poles located further from physical scattering. A handful have a second sheet II pole (visible in the lower panel as the dashed points) but again these are rather distant from physical scattering.
Focussing on the consistent sheet II pole, we may determine couplings to the ππ, KK and ηη channels by factorizing the residues of t(s) at the complex pole position. The result of doing so is plotted in Figure 16 observe ππ and KK couplings of comparable magnitudes, and a coupling to ηη that is somewhat smaller, albeit with some scatter under changes in parameterization.

The KK threshold region
The f 0 resonance pole described in the previous section dominates the amplitude in the energy region around the KK threshold, and in this region, the effect of the distant σ bound-state is limited to providing a smoothly varying 'background'. Given this, it is worthwhile to attempt to describe just that part of the spectrum lying above a t E cm = 0.17 using amplitudes that need not lead to an explicit σ bound-state pole. In Figure 17 we show 9 parameterizations which describe 41 levels in the energy region 0.17 < a t E cm < 0.24, all with χ 2 /N dof < 1.05. As in the previous section we note that the degree of coupling of the ππ, KK sector to the ηη sector is somewhat imprecisely determined, but that otherwise there is very little variation in amplitude with change in parameterization. There is again always a sheet II pole, but we note that it is systematically at a slightly lower mass than in the previous section. We observe there to be somewhat less scatter in the ππ and KK couplings, which show a small systematic shift in phase with respect to the previous section, but which have very similar magnitudes.
Note that some, but not all, of these amplitudes do feature a bound-state pole in roughly the position of the σ, but that this pole position is not precisely determined due to the energy levels below a t E cm = 0.17 not being included in the fit. In those amplitudes which do not feature a bound-state pole, the effect of the σ is being handled by smooth energy dependences that we might think of as 'background'.

Controlling resonant pole content with Jost functions
As presented in the previous two sections, it appears that the lattice QCD spectra are best described by amplitudes which feature a sheet II pole lying close to the KK threshold (as well as a bound-state σ pole at much lower energy), and if they feature poles on sheet III these are distant from physical scattering. This was all determined 'after-the-fact', as the K-matrix forms used do not provide explicit control over the distribution of pole singularities. The position of the poles follows from a potentially complicated interplay of parameter values that is only determined once the fit is complete. On the other hand, Jost functions [77][78][79][80] offer a parameterization of coupled-channel scattering in which the position of resonance poles in the multi-sheeted complex plane can be specified explicitly. We previously made use of such forms to describe two-channel πη, KK scattering in [42], and indeed their implementation is much simplified for twochannel scattering, compared to three-channel scattering. Since we have found in previous sections that ηη appears to be weakly coupled, we choose to eliminate it here by excluding those levels identified previously as having large overlap onto ηη operators (those colored blue in Figure 1). Furthermore, we avoid the complication of simultaneously describing the σ and the f 0 by restricting our attention to those levels in the energy region, 0.17 < a t E cm < 0.24, around the KK threshold. We choose to make use of a conservative set of 30 levels.
Before launching into a Jost function analysis, we first check that two-channel K-matrix analysis can successfully describe the spectrum with "ηη" levels excluded. As an example, we find that a parameterization with symmetric K −1 having independent linear behavior (a + bs) in each element is able to describe the spectrum with χ 2 /N dof = 26.8/(30 − 6) = 1.12. This amplitude is shown in Figure 18 where it is seen to closely resemble the ππ, KK part of previous successful three-channel fits. This suggests that it is indeed reasonable to consider ηη to be decoupled.
The two-channel S-matrix is expressed in terms of the Jost determinant function J by where k 1 and k 2 are the first and second channel cmframe momenta (k 1 = k ππ , k 2 = k KK in the current application). It is convenient to write this as a function of a single kinematic variable ω, defined by as then the S-matrix elements can be expressed as A convenient parameterization of D(ω), similar to that used in [81], features a product of zeroes, giving rise to Smatrix poles, multiplied by a smooth background function to describe the tail of the σ bound-state, (8) A constant term γ 0 does not appear since it cancels in the ratios in Eq. 7. The real-analytic nature of the S-matrix implies that D(ω) = D (−ω ), which fixes Re(γ b odd ) = 0 and Im(γ beven ) = 0. This form does not contain any additional singularities in the energy region around KK threshold beyond the poles whose positions are specified by the values of the complex parameters ω p .
The results of previous sections lead us to believe that a parameterization with just a single pole may be capable of describing the spectra in the limited energy region around the KK threshold. An implementation of Eq. 8 featuring one pole and with three terms in the background polynomial describes the spectrum with χ 2 /N dof = 28.8/(30 − 5) = 1.15. Figure 18 shows the resulting amplitude which is quite similar to those we have previously seen. The position of the pole is allowed to float in the fit, and ends up in the expected position on sheet II, as shown in Figure 19.
Attempting a description using Eq. 8 with two poles and two terms in the exponentiated background polynomial, we find the spectrum is described with χ 2 /N dof = 29.4/(30 − 6) = 1.23. Having allowed both pole positions to float in the fit, we find one ends up in the expected position on sheet II, while the other is poorly determined and appears to be distant on sheet III. The amplitude is plotted in Figure 18 and the poles are shown in Figure 19.
As in our previous three-channel analysis, we observe that restricting the energy region under consideration to be around the KK threshold, discarding energy levels at lower energy, tends to cause the sheet II pole to move to a slightly lower mass.
In Appendix A we consider what happens if we force the amplitude to feature a sheet III pole as well as the sheet II pole demanded by the spectrum. As we move the sheet III pole close to physical scattering, the amplitude comes to have a rapid energy dependence which is not supported by the lattice QCD spectra, leading to an unacceptably large χ 2 . For a more distant sheet III pole, the background part of the amplitude parameterization is able to largely compensate, and the χ 2 remains acceptable. Gray point shows the pole from a typical amplitude taken from Section IV A 3 where three-channel amplitudes described all levels in the region 0.17 < atEcm < 0.24. Remaining colored points correspond to the amplitudes described in the text and plotted in Figure 18 5. S-wave resonance pole summary All successful descriptions of the lattice QCD spectra either over a large energy region, or restricted to the region around the KK threshold, feature a pole on sheet II(−, +, +) at roughly the same position. This pole is found to have large couplings to both ππ and KK. A mirror pole on sheet (−, +, −) is also present. The lattice QCD spectra can tolerate in addition a fairly distant pole on sheet III(−, −, +), but if present it does not appear to be a dominant feature in the amplitude.
Our best estimate for the properties of the sheet II pole are a t √ s 0 = 0.2060(80) − i 2 0.032(12) a t |c ππ | = 0.125(25) a t |c KK | = 0.150 (20) a t |c ηη | = 0.090(35), where the quoted uncertainties take into account the variation over parameterization presented in this section.

B. D-wave resonances
The pole content of the D-wave amplitude can be guessed quite easily from the energy dependence displayed in Figure 14. In order to have such sharp peaks, there must be pole singularities on nearby sheets, and above all three thresholds (ππ, KK and ηη), the nearest unphysical sheet is (−, −, −). Hence we would expect there to be two poles on this sheet, with the one having lower mass being somewhat further from the real axis, corresponding to the larger width of the lower peak.
Indeed, when we examine the pole content of the amplitude described in Section III C, presented in Eq. 4, we find two poles on sheet (−, −, −). As shown in Figure 20, these poles (shown in orange) have 'mirrors' on other unphysical sheets. The lower mass pole, which dominantly couples to ππ, and which we will label " f a 2 ", has mirrors on sheet II and sheet III (not visible in the plot as it lies almost exactly underneath the (−, −, −) pole). The higher mass pole, which couples dominantly to KK, and which we will label " f b 2 ", has mirrors on sheet IV and sheet III (not visible in the plot as it lies almost exactly underneath the (−, −, −) pole). The mirror poles lie at almost the same position owing to the relatively small couplings to sub-dominant channels.
Focussing on the (−, −, −) poles, since these are the closest to physical scattering, we show the variation with parameterization presented in Section III C in Figure 21, which we observe to be extremely small. Figure 22 shows the couplings extracted from the residues of poles on the (−, −, −) sheet. Again, the variation with parameterization change is extremely small, and in all cases the ηη coupling is compatible with zero. As one would expect from Figure 14 where the uncertainties include the small variation with parameterization form described above.
In previous sections we have described a study of coupled ππ, KK, ηη scattering with isospin=0 (G-parity positive). With J P = 0 + we found two singularities are required: a bound-state σ pole appearing on the physical sheet at √ s 0 = 745 (5)  The σ dominates low-energy ππ elastic scattering, giving rise to an S-wave scattering length of m π a = −2.8 (3). As reported on in [49], as the light quark mass is reduced, this meson evolves from being a stable bound-state into a broad resonance.
The f 0 resonance, lying in the KK threshold region, has not previously been observed in a first-principles QCD calculation. It appears to share many of the properties of the experimental f 0 (980) resonance, and we suggest that what we are observing may well represent the evolution of this meson as the light quark mass increases. The absence of a nearby sheet III mirror to the f 0 sheet II pole appears to support a longstanding suggestion that the f 0 resonance may be dominated by KK-molecule configurations (see for example Refs. [11,81,82]). The logic is that a KK molecular state bound by long-range inter-meson forces would be a stable bound-state lying just below the KK threshold were it not for the kinematically open ππ channel into which it decays. The presence of such a decay moves the pole off the real energy axis into sheet II. This is to be compared to a compact state bound by confining interquark forces (whether qq or tetraquark or an even higher quark-gluon Fock state), which is expected to manifest itself as 'mirror' poles on both of sheets II and III.
Presently, the absence of a formalism describing scattering of more than two hadrons prevents us from considering higher-mass scalar meson resonances. Experimentally, these are seen to have dominant decays to four-meson final states. Fortunately, the formalism has been under rapid development [43][44][45][46][47][48], and it is hoped that a final result for at least three-body decays will be available shortly.
With J P = 2 + we isolated the presence of two narrow resonances lying well above the ππ, KK and ηη thresholds, but with negligible couplings to the ηη channel. Because they lie significantly above thresholds, it makes sense to speak of 'branching fractions' for their decay. We compute these using the approach outlined by the PDG [26], where the real and imaginary parts of the pole position, and the couplings extracted from the factorized residue at the pole are used in We note in passing that this approach does not guarantee that the sum of branching fractions be 100%. Our two determined resonances have the following properties: f a 2 : These resonances, computed with heavier than physical light quarks, may be compared to the experimental resonances, f 2 (1270) and f 2 (1525) [26]. The experimental f 2 (1270) has a total width ∼ 190 MeV and decays 84% of the time to ππ and only 5% to KK. The f 2 (1525) is narrower, Γ ∼ 80 MeV, and decays to KK with a 90% branch, and to ππ less than 1% of the time.
These states are often considered to be exemplars of the phenomenological 'OZI' rule of meson decays, which posits that decays proceeding through annihilation of existing quark-antiquark content are suppressed with respect to decays in which extra quarks are generated. In this case this would suggest that f 2 (1270) ∼ 1 √ 2 uū + dd with the decays to ππ and KK being 'OZI-allowed' through creation of extra light or strange quark-antiquark pairs respectively 14 , while f 2 (1525) ∼ ss will decay only to KK through creation of extra light quark-antiquark pairs, with ππ requiring the initial ss to annihilate -an 'OZIsuppressed' decay. This logic can be turned around and used to infer a resonance's quark content on the basis of its preferred hadronic decays.
It is interesting to observe that the results of our calculation (at m π ∼ 391 MeV) appear to support this picture with the f a while the f b 2 closely resembles the f 2 (1525). A somewhat non-rigorous suggestion for the quark-antiquark content of these two states can be obtained by examining the overlaps presented in Figure 4, in particular the 16 3 spectrum in the T + 2 irrep, where two finite-volume states are observed, far from any non-interacting energy levels, but close to the resonance masses for f a 2 and f b 2 . The lighter of the two is seen to have large overlap ontoqΓq operator constructions built from light-quarks, while the heavier state dominantly overlaps with those made from strange-quarks.
We note in passing that the presence of two narrow J P = 2 + resonances with this quark content was anticipated in the simpler calculation of the isoscalar meson spectrum presented in [61], where no meson-meson-like operators where included. Some discussion of the justification for expecting the results presented in [61] to be a reasonable guide to the narrow resonance content of QCD is presented in [71].
Experimentally, the f 2 (1270) is observed to have negligible coupling to ηη, and only a ∼ 10% branch to ππππ, which is not kinematically open at m π ∼ 391 MeV. The f 2 (1525) has a ∼ 10% branch to ηη -this small but significant coupling does not appear to be present at the heavier light quark mass considered in this paper, despite the phase-space for the decay being quite similar to the experimental case.
A plausible method to more rigorously study the internal quark-gluon structure of resonance states of the type we have considered here is to compute their formfactors, by coupling external currents into meson-meson scattering amplitudes. The possibility of studying formfactors of unstable states from lattice QCD was recently reviewed in [71], building on ideas first presented by Lel-louch and Lüscher [83] for the process K → ππ. The 'elastic' form-factors of resonant states could be extracted from amplitudes computed at real energies, by extrapolating to the complex pole positions [84,85]. Exploiting the flexibility of lattice QCD calculation, one could study flavor-and spin-dependence of these form factors, for example by introducingsΓs andūΓu +dΓd currents separately, for various choices of Γ. The distillation technology for such calculations has already been developed [62], and the first nontrivial test has been carried out in the resonant πγ → ππ amplitude [63,64] using the formalism first presented in [86] for transition processes. Fig. 23 presents a summary of the main results of the calculation reported on in this paper. The scalar and tensor ππ, KK amplitudes are plotted, along with the corresponding pole structure, and determinations of the relative strengths of coupling of the resonances to their decay channels.

VI. THE LIGHTEST SCALAR MESONS OF QCD AT mπ ∼ 400 MEV
With the results for isoscalar mesons presented in this paper, taken together with the results in Refs. [39,40,42] for isovector and strange mesons, we have what could constitute a complete nonet of scalar (σ, f 0 , a 0 , κ) and tensor (f 2 , f 2 , a 2 , K 2 ) mesons. It is appropriate at this stage to consider to what extent they have common properties that justifies associating them in this manner.
In the scalar sector, the two states which appear to be mostly closely connected are the f 0 and the a 0 . While their appearance in the relevant elastic amplitude (ππ → ππ or πη → πη) is superficially rather different, being a sharp cusp-like peak for the a 0 and a narrow dip for the f 0 , both effects appear close to the KK threshold, and the corresponding resonance poles are found at positions whose real parts are in close agreement, In addition, the couplings of the resonances to their decay channels agree within statistical uncertainties, Because isospin is an exact symmetry in our calculations, and the effects of electromagnetism are not included, there can be no mixing between the a 0 and the f 0 .
Where the pole singularities differ is in their sheet location and distance into the complex plane. The single relevant pole for the a 0 is located close to the real axis on sheet IV (where Im k πη > 0, Im k KK < 0) while the f 0 pole lies further into the complex plane on sheet II (Im k ππ < 0, Im k KK > 0). Considering the similarities in mass and couplings for the f 0 and a 0 this difference might be considered surprising, but it likely has a relatively simple explanation arising from the only major difference between these two cases: the amount of phase-space for the resonance to decay to the lowest threshold channel, where we have ρ ππ > ρ πη . We can illustrate the effect this has using the simple example of a two-channel Flatté amplitude, in which all elements of the t-matrix have a denominator, where g 1 , g 2 are real valued couplings to channels labelled 1, 2. In the case of a resonance having a relatively small width, the Flatté amplitude has pole singularities at where ρ i is the phase-space for channel i evaluated at s = m 2 0 . If we use the magnitudes of the couplings we extracted from the residues of poles, as in Eq. 10, as the values g 1 , g 2 , we have for both a 0 and f 0 , g 1 /g 2 ∼ 0.8. However, due to the difference in phase-space, we have g1 g2 2 ρ1 ρ2 being slightly less than unity for the a 0 , but somewhat larger than unity for the f 0 , explaining the difference in sheet location. In addition this analysis also offers an explanation for the smaller total width of the a 0 as being due to the near exact cancellation between ρ2 . Within this simple model, there is also an additional sheet III pole, but we observe that it will be much further into the complex plane due to the two terms summing in 1 + g1 g2 2 ρ1 ρ2 , which matches our results which suggest that any sheet III pole lies far into the complex plane.
As argued in the previous section, the dominance of a single pole close to the KK threshold may suggest an association with a KK molecular configuration. Given the apparent parallels between the a 0 and f 0 resonances above, it would appear that at this heavier-than-physical light-quark mass, these states may have a common source, being isovector and isoscalar manifestations of the same KK bound system.
From the significant volume dependence of the lowest lying energy level seen in Fig. 6, one might conclude that the σ could be a physically large object, of size comparable to the lattice volumes used [87], and therefore that it can be interpreted as a molecular ππ state, rather than a compact object bound by confining forces. Since the σ is bound for these values of the quark masses, one can directly apply Weinberg's compositeness criterion [88] to this state. This approach relates the S-wave effective range parameters (k cot δ 0 = a −1 + 1 2 rk 2 + . . .) to the probability, Z, of finding the σ in an elementary bare-particle state, where = 2m π − m σ is the binding energy of the σ, and where each equation potentially receives corrections of order of the range of the ππ interaction. If Z = 1 the σ is purely a compact state, while if Z = 0 it is purely a ππ molecule. Using the values obtained for a and r presented in Section III A, and a binding energy determined from the σ pole mass, Eqs. 11 and 12 give compatible estimates of Z ∼ 0.3 (1). This suggests that for quark masses where m π ∼ 391 MeV the σ can be understood as being predominantly a ππ molecule. While we have successfully determined a complete nonet of scalar mesons with m π ∼ 391 MeV, it remains to be seen how these states evolve with quark mass. To date, the only one of these states that have been studied using smaller values of the quark mass is the σ, where in Ref. [49] it was found that when the pion mass is approximately 236 MeV, the σ has become a broad resonance, resembling somewhat the experimental situation. This is consistent with the expectation that the σ should become increasingly unstable as its phase-space for decay to ππ opens up. With this same line of logic, one might also expect the κ pole, which at m π ∼ 391 MeV resides below threshold on the unphysical sheet, to become a complex-valued pole on the second sheet above threshold as the pion mass is decreased. Similarly, one can expect the phase space of the a 0 and f 0 resonances for decay to the lower channel to increase, and following the discussion above of these poles in the context of the Flatté parametrization, if the decay channel couplings have relatively mild quark-mass dependence, one would expect the poles for f 0 and a 0 to both come to reside on sheet II.
It is also interesting to consider what would happen to the scalar nonet if the light quark mass were increased. For instance, does the κ become a real bound state below πK threshold on the physical sheet for increasingly heavy quarks as is expected in unitarized chiral perturbation theory (UχPT) [89]? What happens if the light quark mass is increased until it is equal to the strange quark mass, meaning the theory has an exact SU(3) flavor symmetry? In this limit there is no splitting between the ππ, KK, ηη thresholds, and the scalar nonet is split into an SU(3) octet and a singlet. The evolution of the nonet in this limit has been previously studied using UχPT [90], where a particular trajectory was chosen which transformed from the physical point to the m π = m K = m η = 300 MeV point. A smooth evolution of the σ pole to a singlet pole lying well below threshold was observed, while all other (octet) states move to a common complex-value pole on the unphysical sheet.
The lightest set of tensor mesons determined at this value of the light-quark mass,  (14) MeV, suggests a rather simple interpretation in terms of qq-like states, lightly "dressed" by their meson-meson decays.
The masses and dominant decays 15 of these states suggest a picture where f a 2 ∼ 1 √ 2 uū + dd , a 2 ∼ 1 √ 2 uū − dd , K 2 ∼ us, f b 2 ∼ ss, with negligible flavor mixing in the isoscalar sector. The small mass difference between f a 2 and a 2 , despite them being constructed from the same quarks, can be ascribed to the small contribution of disconnected diagrams to the f a 2 where the quarks annihilate. The larger width of the f a 2 can be explained by it having allowed decays to ππ with a large phase-space, a channel which is not allowed for the G-parity negative isovector a 2 .

VII. SUMMARY
In this paper, we present the first study with firstprinciples QCD of low-energy isoscalar J = 0 and J = 2 coupled ππ, KK, ηη scattering amplitudes, and their resonance content. This, in conjunction with previous works [39,40,42,49], allows us to paint a complete picture of two different low-lying SU(3) nonets at m π ∼ 400 MeV.
The tensor mesons manifest as clear narrow bumps, and from their decay couplings and masses, one can conclude that they behave like a quark-model qq nonet.
The scalars prove to be far richer in structure. The f 0 and a 0 both appear right at KK threshold, but they manifest themselves differently -the a 0 appears as an asymmetric peak, while the f 0 appears as a dip in a broad enhancement -the source of this difference can be traced to the interference between the f 0 and the σ, and the absence of a σ-like state in the isovector channel. Nevertheless, the f 0 and a 0 are found to have very similar pole properties. The σ appears as a bound state, which is likely to be dominated by a ππ molecular configuration, while the κ emerges as a virtual bound-state.
The calculations presented in this paper and in [39,40,42,49] demonstrate that the finite-volume spectrum approach within lattice QCD can expose non-trivial resonance physics that can significantly inform our understanding of hadron spectroscopy within QCD. Using the Jost formalism of Section IV A 4, we can explore the role played by a sheet III pole supplementing the required sheet II pole describing coupled ππ, KK S-wave scattering. Fixing a sheet II pole at a t √ s 0 = 0.1941 − i 2 0.0389, we have scanned over possible positions for a sheet III pole, allowing two parameters in the exponentiated background polynomial to vary at each new position. Figure 24(a) illustrates the resulting variation in the χ 2 describing the lattice QCD spectrum, which indicates that a pole close to physical scattering is disfavored. Figure 24(b), which focusses on the case where Re(a t k KK ) = 0.04, makes clear why a nearby sheet III pole is problematic -it will tend to produce a rapid variation of the amplitudes on the real energy axis which the lattice QCD spectrum is not in agreement with. When a sheet III pole is more distant however, its effect can be largely compensated by the background function, and a reasonable χ 2 attained.
We note that attempts to supplement the sheet II pole with a second pole located on sheet IV lead to very poor descriptions of the spectra or amplitudes which violate the S ππ,ππ ≤ 1 unitarity bound. Our lattice QCD spectrum is clearly incompatible with such a pole distribution.