Axial resonances a1(1260), b1(1235) and their decays from the lattice

The light axial-vector resonances $a_1(1260)$ and $b_1(1235)$ are explored in Nf=2 lattice QCD by simulating the corresponding scattering channels $\rho\pi$ and $\omega\pi$. Interpolating fields $\bar{q} q$ and $\rho\pi$ or $\omega\pi$ are used to extract the s-wave phase shifts for the first time. The $\rho$ and $\omega$ are treated as stable and we argue that this is justified in the considered energy range and for our parameters $m_\pi\simeq 266~$MeV and $L\simeq 2~$fm. We neglect other channels that would be open when using physical masses in continuum. Assuming a resonance interpretation a Breit-Wigner fit to the phase shift gives the $a_1(1260)$ resonance mass $m_{a1}^{res}=1.435(53)(^{+0}_{-109})$ GeV compared to $m_{a1}^{exp}=1.230(40)$ GeV. The $a_1$ width $\Gamma_{a1}(s)=g^2 p/s$ is parametrized in terms of the coupling and we obtain $g_{a_1\rho\pi}=1.71(39)$ GeV compared to $g_{a_1\rho\pi}^{exp}=1.35(30)$ GeV derived from $\Gamma_{a1}^{exp}=425(175)$ MeV. In the $b_1$ channel, we find energy levels related to $\pi(0)\omega(0)$ and $b_1(1235)$, and the lowest level is found at $E_1 \gtrsim m_\omega+m_\pi$ but is within uncertainty also compatible with an attractive interaction. Assuming the coupling $g_{b_1\omega\pi}$ extracted from the experimental width we estimate $m_{b_1}^{res}=1.414(36)(^{+0}_{-83})$.


Introduction
Ab-initio lattice simulations of light mesons are complicated by the fact that the majority of these states are hadronic resonances and therefore decay under the strong interaction. This is the case also for all vector (J P = 1 − ) and axial-vector (J P = 1 + ) mesons. The vector resonance ρ is the only hadronic resonance that has been treated on the lattice by several collaborations [1][2][3][4][5][6], taking into account its unstable nature and extracting its resonance mass as well as the width. First steps in this direction were also made recently for the vector resonance K * (892) [7], where mass and width were extracted.
In the present paper we aim at simulating the ground state axial-vector resonances a 1 (1260) and b 1 (1235) with I G (J P C ) = 1 − (1 ++ ) and 1 + (1 +− ), respectively, taking into account their strong decays a 1 → ρπ and b 1 → ωπ. The experimental mass and width Γ b 1 = 142 ± 9 MeV of the b 1 (1235) are known rather precisely. The determination of the width for the broader a 1 from the lattice is important in view of its large uncertainty Γ a 1 = 250 − 600 MeV assigned by the Particle Data Group [8]. While most experiments and phenomenological extractions agree on the mass of the a 1 leading to a PDG value of m PDG a 1 = 1.230 (40) GeV [8], determinations of Γ a 1 from diffractive processes where the extraction of the resonance parameters has considerable model dependence [9], deviate substantially from an analysis of data obtained from τ → a 1 ν [10]. As an example for data obtained from a diffractive process, a recent COMPASS measurement published in 2010 [11] provides a much smaller error-bar Γ a 1 = 367 ± 9 +28 −25 MeV.
Our simulation is the first attempt at extracting the s-wave scattering phase shifts for ρπ or ωπ channels in lattice QCD. We are interested in the energy region covering the lowest resonance states, i.e., up to O(1500) MeV. In addition to ρπ and ωπ there are further two-meson decay channels to consider. The a 1 can also couple to f 0 (500)π and f 0 (980)π in p-wave, which goes to three pions. The lowest possible state with three pions combining to total spin 1 needs pions with non-vanishing momentum units like π(0)π(1)π(−1) and thus -for our parameters -has an energy clearly above the range considered here. Another possible coupled channel for a 1 and b 1 is KK * in s-wave; its threshold would be close to 1.4 GeV, but it cannot be produced fromqq interpolator in our N f = 2 simulation and it is Zweig-suppressed, therefore we neglect the corresponding interpolator. The b 1 channel also couples to η 2 ρ, its the contribution will be relevant only for energies above O(1500) MeV and we omit this interpolator. 1 To summarize, we assume that the elastic scattering ρπ and ωπ dominates both channels in the energy region of interest and rely on the same assumption when extracting Γ a 1 →ρπ and Γ b 1 →ωπ from the experimental data. We note that the inelastic scattering with coupled channels has not been treated in the lattice simulation yet, while the analytic frameworks for this challenging problem were developed, for example, in [12][13][14][15].
We also assume that ρ and ω are stable in our case as detailed below. To determine the phase shift in the a 1 channel, for example, we first extract the discrete energy levels of the system with I G (J P C ) = 1 − (1 ++ ) and total momentum zero. These levels are related to the resonant state a 1 as well as the two-particle states ρ(p)π(−p) with p ≃ 2π L n. Both, qq as well as meson-meson interpolators are implemented. The energy levels render the scattering phase shifts via Lüscher's relation [16][17][18][19] and the energy-dependence of the phase shift allows the determination of the resonance mass as well as the width.
The decay b 1 → ωπ was simulated on the lattice earlier only by McNeile and Michael [20], where the Wick contractions with back-tracking loops were omitted. The width Γ[b 1 → ωπ] was extracted based on the amplitude method [21] and reasonable agreement with experiment was found (see figure 3 of [20]). The a 1 has never been simulated taking into account its strong decay, with exception of our preliminary results [22].
While in quantum field theory with dynamical quarks one expects to find contributions of the intermediate meson-meson states even in correlators of qq, in practice such signals were not clearly observed. The masses of a 1 and b 1 were previously determined by a number of lattice collaborations within the so-called single-hadron approach. In this approach onlȳ qq interpolators are used. One assumes that the coupling to two-particle interpolators of the type V π (V = ρ, ω) is negligible and that the mass equals the observed ground state energy level m = E (for P = 0). It turned out to be a better strategy to include mesonmeson interpolators explicitly in the set of source and sink operators and to study the meson-meson system this way. Indeed in this more complete, coupled system level shifts as compared to the non-interaction case can be observed. The assumption that two-particle interpolators may be neglected is particularly questionable in our situation, since V (0)π(0) with E ≃ m V + m π < m a 1 ,b 1 is the ground state of the system. An additional complication for the channels of interest is that the scattering particles ρ and ω decay under the strong interaction to 2π and 3π, respectively. We assume stable vector mesons and restrict the simulation of the axial channels to the total momentum P ≡ |P| = 0, where the possible effect of vector meson decay is least significant. Analytic frameworks for the scattering of unstable particles (with some emphasis on s-wave scattering) in the finite volume were recently formulated in [23], while the related challenging problem of three particles in a finite box was analytically considered in [14,[24][25][26]. In Ref. [23] the approach was applied to study the effect of the ρ → ππ decay to the ρπ scattering. Its effect on the discrete spectrum at zero momentum is found to be negligible on our lattice with L ≃ 2 fm and m π ≃ 266 MeV in the region of interest, where a 1 resides [23] 2 . The ρ(0) behaves nearly like a stable particle since its p-wave decay to the lowest two-particle state π(1)π(−1) has sizable energy E π(1)π(−1) ≃ 1.3 GeV [3]. The effect of ω → 3π decay to ωπ scattering has not been analytically studied in detail, but the effect is suppressed again since at least two pions have momentum when ω decays at rest, and it is suppressed due to the small width of ω.
The a 1 and b 1 mesons have been extensively studied also by non-lattice methods. Weinberg derived m a 1 /m ρ ≃ √ 2 based on spectral functions of currents [27]. More recently, axial mesons emerged as dynamically generated resonances in the study of scattering between pseudoscalar and vector mesons based on chiral Lagrangians [28][29][30]. Unitarized Chiral Perturbation Theory models render poles related to axial mesons dynamically although these are not present in the original Lagrangian [31]. This approach was used also to determine the dependence of their masses and widths as a function of number of colors (N c ). It was found that this dependence significantly deviates from the expectation M ∝ N 0 c and Γ ∝ 1/N c for pure qq structure [32]. The chiral properties and issues related to chiral restoration were considered for example in [33,34].
After introducing our lattice setup, we present the masses of the scattering particles π, ρ and ω in section 3. The interpolators and the energy spectra in the a 1 and b 1 channels are discussed in section 4, while the Wick contractions are relegated to the appendix. The analysis related to the scattering phase shifts and the resonance parameters is given in two subsequent sections followed by our conclusions.

Lattice setup
The simulation is based on one ensemble of N f = 2 gauge configurations with clover Wilson u/d quarks, generated by the authors of [35,36]. The dynamical and valence u/d quarks have the same mass corresponding to m π = 266(4) MeV. The parameters of the ensemble are shown in table 1, while more details are given in [3]. Due to the limited data for just a single ensemble, our determination of the lattice spacing a reported in [3] results from taking a typical value of the Sommer parameter r 0 . The uncertainty associated with this choice might lead to a small shift of all dimensionful quantities. m π a m ρ a m ω a 0.1673(16) 0.5107(40) 0.514 (15)  The sea and valence quarks obey periodic boundary conditions in space. Valence quark propagators periodic and anti-periodic in time are combined into so-called "P + A" propagators, which effectively extends the time direction to 2N T = 64 [3]. The errors of our results correspond to single-elimination jack-knife errors. We use translation invariance in t to sum over correlators from sources in every time slice.
A rather small volume V = 16 3 × 32 (L ≃ 2 fm) has a practical advantage since it enforces the scattering particles ρ(p) and ω(p) with p ≡ |p| < 2π/L to be nearly stable, which is due to the sizable difference between energies of decay products π(0) and π(2π/L) at our L. As a consequence, the finite volume energies of ρ(p) and ω(p) with p < 2π/L are relatively close to the resonance mass. The given volume also enables us to use the powerful full distillation method [37], which allows for the computation of all contractions for the correlation matrix with qq and V π interpolators.
3 Masses of π, ρ, ω The masses of the scattering particles π, ρ, ω are needed for the position of the thresholds and we collect them in table 2. We use m π as determined in [3].
The ρ mass in table 2 was extracted as m ρ = E ρ (p = 0) in [3] and is indeed found close to m res ρ in all simulations [1][2][3][4][5][6]. The ω energy is calculated using quark-antiquark interpolators O s=n 1−5 given by Eq. (21) of [3] with five different Dirac/space structures, while flavor is replaced withūu +dd. This approach is consistent with our approximation of treating ρ(p) and ω(p) with p ≪ 2π/L as stable; this holds well for the narrow ω and is commonly applied also for the broader ρ. The distillation method enables straightforward calculation of the disconnected contributions to ω and the final correlators are averaged over all initial time-slices and three polarizations. The disconnected contribution is small and we find m ω ≃ m ρ as expected. The value of m ω follows from a 2-exponent fit in range t = 3 − 12 using interpolators O n 1,2,3,5 and t 0 = 2. This is a conservative choice with a comparatively large uncertainty and is fully consistent with the result of other possible choices.
where x i = (x i , t) and ∇ denotes the covariant derivative. The ρ and π mesons are separately projected to zero momentum in O ρπ . We do not implement the interpolator ρ(1)π(−1) (the argument ±1 indicates momenta ±2π/L) since effects of ρ(1) → π(1)π(0) n t 0 interp. fit range p cot(δ) √ s  Table 3. Energies and phases in the a 1 channel with I G (J P C ) = 1 − (1 ++ ) and P = 0, where a −1 ≃ 1.59 GeV. Both levels were obtained using a 1 exponential fit. The p give the eigen-momenta of the interacting system determined from the energy levels according to Eq. (5.1). The ground state is below ρπ threshold, so p and δ are imaginary, while p cot δ is real.  Table 4. Similar as table 3 but for b 1 channel with I G (J P C ) = 1 + (1 +− ). Both levels were obtained using a 2 exponential fit. The second level is consistent with threshold energy m π + m ω due to relatively large uncertainties of E 2 and m ω ; the corresponding phase with huge errors is consistent with δ ≃ 0 and therefore not provided. The uncertainty in m ω has negligible effect on δ for the second level.
decay are non-negligible [3]. Our analysis will therefore concentrate on the lower energy region E < E ρ(1)π(−1) ≃ 1.7 GeV. Similarly, for the b 1 channel with J P C = 1 +− and |I, where, again, ω and π are separately projected to zero momentum. All quark fields q in (4.1), (4.2) are smeared according to the distillation method [37], thus effectively replaced by Nv where the v (k) denote the Laplacian eigenvectors of the time slice.
The energy spectrum E n is extracted from the correlation matrix averaged over all initial times t i . The Wick contractions for both channels are presented in figures 3 and 4 of the appendix. In particular in the b 1 channel a rather large number of diagrams appears. Expressions for various elements of the correlation matrix in terms of those Wick contractions are also provided in the appendix. We evaluate all the Wick contractions using the distillation method, which handles efficiently also those with backtracking quark lines. The variational method with the generalized eigenvalue equation is applied to extract the discrete spectrum E n [16,[38][39][40]. The resulting eigenvalues λ n (t) ∼ e −En(t−t 0 ) give the effective energies E eff n (t) ≡ log[λ n (t)/λ n (t + 1)] → E n . The spectrum E n is extracted using correlated fits to λ n (t).
The resulting spectrum E n is shown in figure 1, where effective energies are plotted for the cases when O V π is included or excluded from the correlation matrix. The horizontal lines indicate the position of the threshold m V + m π (which has sizable uncertainty in the b 1 channel), and the energy of the non-interacting V (1)π(−1) system.
We concentrate on the spectrum obtained including O V π , which is shown in the first and third pane of figure 1 and listed in tables 3 and 4. The lowest levels (circles) are near the m V + m π threshold, as expected, and its dominant component is the V (0)π(0) two-particle interpolator. The second level (squares) arises due to the presence of a 1 (1260) or b 1 (1235) resonances in these channels. The next level would be expected close to V (1)π(−1), but we do not expect to see it since the corresponding interpolator is not implemented in (4.1) and (4.2).
The third levels in figure 1 in both channels are noisy and unreliable, so we refrain from presenting quantitative results. Both levels correspond to masses close to 2 GeV or above, so we see no indication for the possible existence of a 1 (1420), which was introduced to explain recent preliminary data by COMPASS in a 1 → f 0 π channel [41]. The third level in b 1 channel might be related to the observed b 1 (1960).

Resonance parameters for a 1 (1260)
The position of the 1 ++ ground state below m ρ + m π threshold indicates that the energy of ρ and π is smaller if they are in the box together then if they are in the box alone. The negative energy shift is consistent with an attractive interaction in the resonant a 1 channel.
We proceed to extract the ρπ phase shifts and the resonance parameters of a 1 (1260). Outside the interaction region the mesons are considered as free particles and the energy levels E are related to the momenta p of the two-particle state ρ(p)π(−p) through where we employ the continuum dispersion relation which applies well for the small momenta p < 2π/L of interest [3]. The s-wave phase shifts δ at these values of p are given by the well known Lüscher relation [19] tan δ(p) = √ π p L 2 Z 00 1;   Table 5. The resulting Breit-Wigner masses m res together with the couplings g for a 1 → ρπ and b 1 → ωπ, which are related to the Breit-Wigner width Γ ≡ g 2 p/s. For the resonance masses the second uncertainty given stems from the systematic uncertainty in extracting the first excited state reliably. This uncertainty has negligible effect on the extracted coupling. The experimental values for the couplings g are derived from the measured total widths [8] since the branching ratios to V π have not been measured, but are expected to be largely dominant. The lattice value for the resonance mass of b 1 (1235) is obtained assuming experimental g exp b1ωπ . All results are for our value of m π ≃ 266 GeV.
which applies above and below threshold for elastic scattering. The ground state below the m ρ + m π threshold renders imaginary p and real p cot δ in table 3. The first excited level gives δ ≈ 90 • , thus it is located close to the a 1 (1260) resonance mass and m a 1 ≈ E 2 holds.
The Breit-Wigner parametrization gives Γ a 1 (s) in terms of the phase space and the coupling g a 1 ρπ . We obtain which applies in the vicinity of the resonance. Assuming that (like it is the case for the ρ) linearity (5.4) is a good approximation down to the threshold and slightly below it, we interpolate linearly in s between the two values of p cot δ/ √ s of table 3, as shown in figure   2. From the zero and the slope we obtain m res a 1 and g a 1 ρπ .
The resulting parameters of the a 1 (1260) resonance are compared to experiment in table 5. The value of m res a1 at m π = 266 MeV is slightly higher than that of the experimental resonance a 1 (1260). This first lattice result for g a 1 ρπ is valuable since there is still considerable uncertainty on the total width Γ a 1 and on the a 1 → ρπ branching ratio. We provide the upper limit for g exp a 1 ρπ resulting from the total width Γ exp a1 = 250 − 600 MeV [8] 3 , which agrees with our g a 1 ρπ within the large experimental and theoretical uncertainties. Our lattice result for g a 1 ρπ is also close to the value g phen a 1 ρπ ≈ 0.9 GeV obtained using an Unitarized Effective Field Theory approach [31] (converted to our convention).
The scattering length a ρπ l=0 in table 5 is obtained using the effective range fit p cot(δ) = 1 a 0 + 1 2 r 0 p 2 through two energy levels. Our result at m π = 266 MeV is close to a ρπ l=0 (m phys π ) ≈ 0.37 fm obtained from Unitarized Effective Field Theory [42], while the corresponding experimental value is not known. However the exact position of the central value for the ground state E 1 in figure 1 with respect to the threshold m ω + m π shows a slight disagreement with the expectation. It is expected to be slightly below threshold due to an attractive interaction in the resonant channel. Since Γ b 1 < Γ a 1 one expects a smaller size of the energy shift ∆E 1 = E 1 −m V −m π in the b 1 channel than in the a 1 channel, rendering a lattice extraction challenging. We estimate the expected energy shift a∆E 1 ≃ −0.01 based on g exp b 1 ωπ , Breit-Wigner dependence (5.4), the Lüscher relation (5.2) and the value of m res b 1 (6.1) determined below; note that this shift is smaller than the uncertainty of the ground state energy level in the b 1 channel and comparable to the uncertainty in am ω . A correlated analysis reveals that our ground state E 1 is compatible with m ω + m π within the large uncertainties, although the central value leads to E 1 m ω + m π . We stress that this discrepancy with expectations is not statistically significant.
If upon improving the statistics the ground state level still stays above threshold, this would be hard to understand. The three-pion state, for example, could hardly explain such behavior since the lowest state π(0)π(1)π(−1) with J = 1 has energy ≃ 1.6 GeV, which is significantly above the threshold on our lattice. Due to the discussed uncertainty of the ground state energy level we determine the resonance b 1 (1235) mass using the Breit-Wigner relation (5.4), Br[b 1 → ωπ] ≃ 1, which has not been measured but is expected to be valid to a good approximation. The resulting resonance mass at our m π is somewhat higher than the experimental one.

Conclusions and outlook
We presented the first lattice simulation of light axial resonances a 1 and b 1 taking into account their dominant strong decays modes ρπ and ωπ. The interpolating fields qq as well as V π (V = ρ, ω) were used for this purpose. We find an energy level near the V (0)π(0) threshold and excited levels close to the resonance positions of a 1 (1260) and b 1 (1235). The s-wave phase shifts for V π scattering were extracted using the Lüscher relation and they are fitted using the Breit-Wigner formula. The resonance parameters are compared to experiment in table 5.
The resonance masses m res a 1 and m res b 1 at our m π ≃ 266 MeV are somewhat higher than the experimental values. The a 1 → ρπ coupling g a 1 ρπ , which parametrizes the corresponding decay width, agrees with experiment within sizable error bars. The analogous b 1 → ωπ coupling was not extracted since the central value for the ground state in the b 1 channel is found slightly above m ω + m π , while one expects it to be slightly below threshold in an attractive channel. The ground state is still consistent with m ω + m π and an attractive interaction within the sizable uncertainty. Future lattice studies with improved statistical accuracy could clarify this issue by resolving smaller energy shifts which should lead to a determination of g b 1 ωπ from Lattice QCD.
Our analysis assumes stable scattering particles ρ and ω, which is a reasonable approximation for our simulation with zero momentum, L ≃ 2 fm and m π ≃ 266 MeV. However, the effects of vector meson decays might have to be taken into account in future simulations with physical pion masses, larger box sizes and possibly non-zero total momenta. Proper treatment of this problem appears considerably more challenging. An analytic framework for the scattering of unstable particles in finite volume has been recently proposed in [23] while the related problem of three particles where two of them can resonate was considered in [14,[24][25][26]. These approaches can serve as a guideline for future lattice simulations.