Using a new analysis method to extract excited states in the scalar meson sector
Jacob Finkenrath Constantia Alexandrou Joshua Berlin Mattia Dalla Brida Theodoros Leontiou Marc Wagner






§ INTRODUCTION



In this work we explore AMIAS, a rather new analysis method to extract energy differences and
amplitudes from lattice QCD correlators and correlation matrices. The basic idea is to consider
a very large number of fits to these correlators using Monte Carlo methods, which will lead to 
probability distributions for the fit parameters, that are the energy differences and amplitudes.
This avoids the necessity of identifying plateau regions for effective energies. In particular, for unstable 
systems and resonances identifying plateaus can be very challenging, since the signal-to-noise ratio
is typically rather poor for large temporal separations. Moreover it overcomes the limitation of extracting at most
N energy eigenstates, when N interpolating fields are used (as it is e.g. the case, when doing a standard GEVP analysis)

We apply AMIAS to analyze correlators computed for studying the a_0(980) meson, which has quantum
numbers I(J^P) = 1(0^+) and mass m_a_0(980)≈ 980  MeV <cit.>. We
use interpolating fields with both two-quark and a four-quark structures formed by an up-quark u and
an anti-down quark d̅ and except for one case a strange and an anti-strange quark s  and s̅. The four-quarks
are arranged as a meson-meson interpolating field or in a diquark-anti-diquark combination,
which can probe different possible tetraquark structures of a_0(980). We consider six interpolating fields,


    𝒪^1 𝒪^q q̅ = 1/√(V_s)∑_x(d̅( x) u( x)) 
    
          𝒪^2 𝒪^K K̅,  point = 1/√(V_s)∑_x(s̅( x) γ_5 u( x)) (d̅( x) γ_5 s( x)) 
    
           𝒪^3 𝒪^η_sπ,  point = 1/√(V_s)∑_x(s̅( x) γ_5 s( x)) (d̅( x) γ_5 u( x)) 
    
           𝒪^4 𝒪^Q Q̅ = 1/√(V_s)∑_xϵ_a b c(s̅_b( x) (C γ_5) d̅_c^T( x)) ϵ_a d e(u_d^T( x) (C γ_5) s_e( x)) 
    
           𝒪^5 𝒪^K K̅,  2part = 1/V_s∑_ x, y(s̅( x) γ_5 u( x)) (d̅( y) γ_5 s( y)) 
    
          𝒪^6 𝒪^η_sπ,  2part = 1/V_s∑_ x, y(s̅( x) γ_5 s( x)) (d̅( y) γ_5 u( y)) ,


where V_s denotes the spatial lattice volume and C the charge conjugation matrix. All of them couple to a_0(980)
and other states with the same quantum numbers. For example, the interpolating fields 𝒪_5 and 𝒪_6 
mainly generate the two-meson states K + K and  π + η, respectively, which 
are expected to have masses close that of the a_0(980). Note that the interpolating fields 𝒪_2
and 𝒪_3 represent two mesons located at the same point in space with only total momentum zero, but
arbitrary relative momenta involved (a structure resembling a 4-quark bound state). In contrast to that
𝒪_5 and 𝒪_6 correspond to two mesons with both total and relative momentum zero.

We compute the full six by six correlation matrix including both connected and disuconnected contributions,
which we neglect in our previous studies <cit.>.
Moreover here we have increased statistical accuracy of the correlators
C_j k(t) = ⟨𝒪^j(t) 𝒪^k †(0) ⟩,
included propagation of strange quarks within a timeslice and analyzed the correlators with AMIAS. 
We use an ensemble of around 500 gauge link configurations generated with 2+1 dynamical Wilson clover 
quarks and Iwasaki gauge action generated by the PACS-CS Collaboration <cit.>. The lattice size 
is 64 × 32^3 with lattice spacing a ≈ 0.09  fm and pion mass m_π≈ 300  MeV.





§ CORRELATORS ON A PERIODIC LATTICE



A correlator on a periodic lattice with time extent T can be expanded according to


    C_j k(t) = ⟨𝒪^j(t) 𝒪^k †(0) ⟩ = ∑_m,nexp{-ℰ_m (T-t)}⟨ m | 𝒪^j | n ⟩exp{-ℰ_n t}⟨ n | 𝒪^k † | m ⟩/∑_m exp{-ℰ_m T}


with energy eigenstates | m ⟩ and corresponding energy eigenvalues 
ℰ_0 ≤ℰ_1 ≤ℰ_2 … (| 0 ⟩ = | Ω⟩ denotes the vaccum). 
This correlator can also be expressed in a more convenient form,


    C_j k(t) = ⟨𝒪^j(t) 𝒪^k †(0) ⟩ = ∑_m,n c^j_m,n (c^k_m,n)^∗exp{-(ℰ_m + ℰ_n) (T/2)}cosh{Δℰ_n,m (t - T/2)}/∑_m exp{-ℰ_m T}


with c_m,n^j = ⟨ m | 𝒪^j | n ⟩ and Δℰ_n,m = ℰ_n - ℰ_m. For 
example, if j = k, if 𝒪^j probes the sector, which contains energy eigenstate | 1 ⟩, if the quantum 
numbers of | 1 ⟩ are different from those of the vacuum | Ω⟩, and if the correlator is not contaminated
by multi-particle states (see discussion below), (<ref>) reduces to


    C_j j(t) = ⟨𝒪^j(t) 𝒪^k †(0) ⟩≈ 2 |c^j_0,1|^2 exp{-Δℰ_1,Ω (T/2)}cosh{Δℰ_1,Ω (t - T/2)}


for sufficiently large t and T. A standard technique to determine the energy difference Δℰ_1,Ω = E_1 - E_Ω, 
the mass of state | 1 ⟩, from the asymptotic t behavior of C_j j(t) (Eq. (<ref>)) is to fit


    C_j j(t) = A cosh{-Δℰ_1,Ω (t - T/2)}


to the lattice QCD results for C_j j(t) with fitting parameters Δℰ_1,Ω and A. Alternatively, one can also solve the equation


    C_j j(t)/C_j j(t-a) = cosh{E_eff(t) (t - T/2)}/cosh{E_eff(t) (t-a - T/2)}


with respect to E_eff(t), where E_eff(t) ≈Δℰ_1,Ω. In other words, a plateau-like 
behavior of E_eff(t) indicates the mass ℰ_1,Ω. In practice, however, the temporal
extent T of the lattice is limited and the effective energy E_eff(t) is often very noisy for large t, 
rendering a reliable extraction of Δℰ_1,Ω difficult.

When multi-particle states are present in the investigated sector, the determination of low-lying masses is even more difficult.
We sketch this for a simple non-interacting two-meson system, which is generated by an interpolating field 𝒪 = 𝒪^(1)⊗𝒪^(2). Energy eigenstates of such a system can be written as | n ⟩ = | n_1 ⟩^(1)⊗  | n_2 ⟩^(2), where | n_1 ⟩^(1) and | n_2 ⟩^(2) (n_1 , n_2 = 0,1,2,…) are eigenstates of the two Hamiltonians describing the individual mesons. Energy eigenvalues of this two-meson system are ℰ_n_1,n_2^(1+2) = ℰ_n_1^(1) + ℰ_n_2^(2), in particular the lowest energy eigenvalue is ℰ_0,0^(1+2) = ℰ_0^(1) + ℰ_0^(2). Note that in the expansion (<ref>) there is a term, which is proportional to cosh{(ℰ_0^(1) - ℰ_0^(2)) (t - T/2)}, i.e. only decaying with the mass difference of the two mesons. In the region of t ≈ T/2 this unwanted term can be more dominant than the signal term, which is proportional to cosh{ℰ_0,0^(1+2) - ℰ_Ω) (t - T/2)}. Thus, one has to be rather careful, when analyzing correlators of sectors containing multi-particle states. This is illustrated by Fig. <ref>,
where we show the correlator using the interpolating field 𝒪^q q̅ (Eq. (<ref>)). Clearly, 
there are two effective mass plateaus, a E_eff(t) ≈ 0.60 for t/a  7, corresponding to 
the mass of a π + η or a K + K two-meson state, and a E_eff(t) ≈ 0.25 for t/a  8,
which does not correspond to the mass of any state, but to the mass difference m_η - m_π of two single-meson states.










§ AMIAS



Lattice QCD results for the correlators C_j k(t) = ⟨𝒪^j(t) 𝒪^k †(0) ⟩ 
with 𝒪^j, j = 1,…,6 defined in Eqs. (<ref>) to (<ref>) can be parameterized and fitted
the expression given in Eq. (<ref>). Since statistical accuracy is limited, it is sufficient to consider a rather small number 
of energy eigenstates, i.e. ∑_m,n , ∑_m →∑_m,n = 0^N-1 , ∑_m = 0^N-1 with N  10.

With AMIAS <cit.> one can
determine probability distribution functions (PDFs) for the fit parameters Δℰ_n,m (energy differences)
and c^j_m,n (amplitudes). AMIAS is able to deal with a rather large number of parameters by using Monte Carlo
techniques. In contrast to e.g. the GEVP, it is not necessary to identify plateau regions or to specify temporal 
fitting ranges.

The PDF for the complete set of fit parameters is given by


    P(Δℰ_n,m , c^j_m,n) = 1/Z e^-χ^2/2


with appropriate normalization Z and


    χ^2 = ∑_j,k∑_t/a=1^(T-1)/a(C_j k^lattice(t)- C_j k(t))^2/(w_j k(t))^2 ,


which is the well-known χ^2 used in χ^2 minimizing fits (C_j k^lattice(t) denote lattice
QCD results for correlators with corresponding statistical errors w_j k(t)). To obtain the PDF for a 
specific fit parameter, one has to integrate in (<ref>) over all other parameters. In particular, the
probabilty for parameter 𝒜_j to be inside [a,b] (𝒜_j represents either one of the energy 
eigenvalue differences Δℰ_n,m or amplitudes c^j_m,n) is


    Π(𝒜_j ∈ [a,b]) = ∫_a^b d𝒜_j  ∫_-∞^+∞∏_k  j d𝒜_k  
               e^-χ^2/2/∫_-∞^+∞∏_k d𝒜_k  
              e^-χ^2/2 .


This multi-dimensional integral can be computed with standard Monte Carlo methods. We implemented 
a parallel-tempering scheme combined with the Metropolis algorithm as described in detail in 
Ref. <cit.>. The parallel-tempering scheme prevents that the algorithm is stuck in a 
region around a local minimum of χ^2 and, thus, guarantees ergodicity of the algorithm.

We have again analyzed the correlator of 𝒪^q q̅ shown in Fig. <ref> (left), this time using AMIAS with fit function


    C(t) = ∑_n=1^2 A_n cosh{Δℰ_n (t - T/2)} .


The resulting PDFs for the four parameters Δℰ_1, Δℰ_2, A_1 and A_2 are shown in Fig. <ref>.
As in Fig. <ref> (right) two energy differences can clearly be identified, a Δℰ_1 ≈ 0.25 
(the mass difference m_η - m_π) and a Δℰ_2 ≈ 0.60 (the mass of a π + η or K + K two-meson state).










§ ANALYSIS OF THE 6 × 6 CORRELATION MATRIX



We now use AMIAS to analyze the 6 × 6 correlation matrix with interpolating fields (<ref>) to (<ref>).
The fit function is obtained by restricting (<ref>) to a finite number of terms,


    C_j k(t) = ∑_n=1^N c^j_n (c^k_n)^∗cosh{Δℰ_n (t - T/2)} .


In a first analysis we use lattice QCD data, where the propagation of strange quarks within the same 
timeslice has been neglected. The interpolating field
𝒪^q q̅ then decouples and, hence, cannot be considered in a tetraquark study 
of a_0(980) (cf. also Refs. <cit.>, where the same data has been used). 
Moreover, we do not consider the diquark-antidiquark interpolating field 𝒪^Q Q̅. In 
Fig. <ref> (top) we show the four lowest masses obtained with N = 8 terms in 
(<ref>)[AMIAS is able to determine the optimal N automatically. For details we refer to <cit.>.], 
corresponding to the expected two-particle states π + η and K + K with both mesons at rest
(Δℰ_2 and Δℰ_3) as well as with one quantum of relative momentum (Δℰ_4 and Δℰ_5).
Note that AMIAS also finds the energy eigenvalue difference m_η - m_π (Δℰ_1, not shown in Fig. <ref>),
which has already been discussed in previous sections.



Since (c^j_n)^∗ = ⟨ n | 𝒪^j † | Ω⟩, the amplitudes extracted with AMIAS are the coefficients of 
the expansions of the trial states 𝒪^j † | Ω⟩ in terms of the extracted energy eigenstates | n ⟩, i.e.

    𝒪^j † | Ω⟩≈∑_n=1^N (c^j_n)^∗ | n ⟩ .


More interesting, however, are the opposite expansions, i.e. the extracted energy eigenstates in terms of the trial states,


    | n ⟩≈∑_j=2,3,5,6 v^j_n 𝒪^j † | Ω⟩ .


It is easy to show that the matrix formed by the coefficients v^i_n is the inverse of the matrix formed by the coefficients (c^j_n)^∗ (j = 2,3,5,6, n = 2,…,5) up to 
exponentially small corrections, i.e. ∑_j v^j_m (c^j_n)^∗ = δ_m,n. Alternatively, 
one can also use the AMIAS samples for Δℰ_n and c^j_n and Eq. (<ref>) to reconstruct a correlation matrix for each sample. Solving standard GEVPs for each reconstructed 
correlation matrix yields eigenvectors with components identical to v^j_n. The coefficients |v^j_n| are shown in Fig. <ref> (bottom) (not the PDFs, but the most likely values).
Clearly, Δℰ_2 and Δℰ_3 correspond to two-particle states π + η and K + K with 
both mesons at rest, since the coefficients v^j_n show almost exclusively overlap to the trial states
𝒪^η_sπ,  2part† | Ω⟩ and 𝒪^K K̅,  2part† | Ω⟩. Δℰ_4 
and Δℰ_5 are close to the expected two-particle states with one quantum of relative momentum and the overlaps
shown by v^j_n are consistent with this interpretation. There is no indication for any additional state in the region of 1000  MeV, 
which could correspond to a_0(980). These findings are consistent with our previous study using ETMC gauge link configurations <cit.>.

When including the propagation of strange quarks within the same timeslice, the interpolating field 𝒪^qq̅ couples 
to the four-quark interpolating fields given in Eqs. (<ref>) to (<ref>). This introduces, however, a lot of statistical noise and 
it is rather difficult to resolve energy differences precisely. In Fig. <ref> we show preliminary AMIAS results for
the full 6 × 6 correlation matrix, where N = 12 terms in Eq. (<ref>) have been used. The two lowest energy eigenstates exhibit
rather clear signals and correspond to the π + η and K + K two-meson states. Even though higher excitations are less prominent,
we find strong indication for an additional state with mass Δℰ_3 ≈ 0.6 / a ≈ 1300  MeV, which
is below the expectation for the two-particle states with one quantum of relative momentum and, hence, might be a candidate for the a_0(980).








§ CONCLUSIONS


With AMIAS one can extract energy differences and amplitudes without the necessity to identify pleateau regions or to specify temporal 
fitting ranges. We have shown that AMIAS can successfully be used to analyze correlators with strong contributions from several energy 
eigenstates as well as rather noisy correlation matrices. An example for the latter is e.g. the full 6 × 6 correlation matrix 
for the a_0(980) channel. Another major advantage of AMAIAS might be, that it does not require all elements of a correlation matrix, 
i.e. particularly noisy elements or elements, which are very time-consuming to compute, can be omitted. This we plan to address in a future publication.







Acknowledgements:
M.W. acknowledges support by the Emmy Noether Programme of the DFG (German Research Foundation), 
grant WA 3000/1-1. This work was supported in part by the Helmholtz International Center for FAIR within the framework 
of the LOEWE program launched by the State of Hesse.
This work was cofunded by the European Regional Development Fund and the Republic of Cyprus through the Research
Promotion Foundation (Project Cy-Tera NEA YΠOΔOMH/ΣTPATH/0308/31) by the grand cypro914.




99

Agashe:2014kda 
  K. A. Olive et al. [PDG Collaboration],
  
  Chin. Phys. C 38, 090001 (2014).
  
  

Abdel-Rehim:2014zwa 
  J. Berlin et al.,
  
  PoS LATTICE 2014, 104 (2014)
  [arXiv:1410.8757 [hep-lat]].
  

Berlin:2015faa 
  J. Berlin et al.,
  
  PoS LATTICE 2015, 096 (2015)
  [arXiv:1508.04685 [hep-lat]].
  

Berlin:2016zci 
  J. Berlin et al.,
  
  PoS LATTICE 2016, 128 (2016)
  [arXiv:1611.07762 [hep-lat]].
  

Aoki:2008sm 
  S. Aoki et al. [PACS-CS Collaboration],
  
  Phys. Rev. D 79, 034503 (2009)
  
  [arXiv:0807.1661 [hep-lat]].
  

Alexandrou:2008bp
  C. Alexandrou, C. N. Papanicolas and E. Stiliaris,
  
  PoS LATTICE 2008, 099 (2008)
  [arXiv:0810.3982 [hep-lat]].
  

Papanicolas:2012sb 
  C. N. Papanicolas and E. Stiliaris,
  
  arXiv:1205.6505 [hep-ph].
  

Alexandrou:2014mka 
  C. Alexandrou, T. Leontiou, C. N. Papanicolas and E. Stiliaris,
  
  Phys. Rev. D 91, 014506 (2015)
  
  [arXiv:1411.6765 [hep-lat]].
  

Alexandrou:2012rm 
  C. Alexandrou et al.,
  
  JHEP 1304, 137 (2013)
  
  [arXiv:1212.1418 [hep-lat]].
  





