Pole position of the $a_1(1260)$ from $\tau$-decay

We perform an analysis of the three-pion system with quantum numbers $J^{PC}=1^{++}$ produced in the weak decay of $\tau$ leptons. The interaction is known to be dominated by the axial meson $a_1(1260)$. We build a model based on approximate three-body unitarity and fix the free parameters by fitting it to the ALEPH data on $\tau^-\to \pi^-\pi^+\pi^-\,\nu_\tau$ decay. We then perform the analytic continuation of the amplitude to the complex energy plane. The singularity structures related to the $\pi\pi$ subchannel resonances are carefully addressed. Finally, we extract the $a_1(1260)$ pole position $m_p^{(a_1(1260))}-i\Gamma_p^{(a_1(1260))}/2$ with $m_p^{(a_1(1260))} = (1209 \pm 4^{+12}_{-9})\text{MeV}$, $\Gamma_p^{(a_1(1260))} = (576 \pm 11 ^{+89}_{-20})\text{MeV}$}.

We perform an analysis of the three-pion system with quantum numbers J P C = 1 ++ produced in the weak decay of τ leptons. The interaction is known to be dominated by the axial meson a1(1260). We build a model based on approximate three-body unitary and fix the free parameters by fitting it to the ALEPH data on τ − → π − π + π − ντ decay. We then perform the analytic continuation of the amplitude to the complex energy plane. The singularity structures related to the ππ subchannel resonances are carefully addressed. Finally, we extract the a1(1260) pole position m

I. INTRODUCTION
The internal dynamics of the Quantum Chromodynamics (QCD) degrees of freedom manifests itself in the spectrum of hadron resonances. The mass of a resonance characterizes the energy of the excitation while its width reflect on the coupling to the decay channels. The meson spectrum has been qualitatively elucidated by the quark model [1] and recently, at least for some states, calculations based on first principles lattice QCD are becoming available [2,3]. For a majority of states, however, ab initio QCD calculations of their decay properties, e.g. decay widths, branching ratios, are not yet available. Pushing such calculations forward is important given the growing body of evidence for novel hadronic phenomena [4][5][6][6][7][8][9][10], e.g. the X, Y, Z states observed in heavy quarkonia [11][12][13]. Many of these new states are observed in decays to three-particle final sates.While hadron scattering involving two stable particles is rather well understood formally, the methodology for incorporating three and more particles is still being developed both in the infinite volume [14][15][16][17][18] and finite volume [19][20][21][22][23][24].
A large number of light meson resonances dominantly decay to three pions. This includes the enigmatic a 1 (1260) resonance, which is the lightest axial vector meson with J P C = 1 ++ . The properties of the a 1 resonance are difficult to assess, due to its large width that is affected by the three-pion dynamics. The ππ subchannel is dominated by the ρ resonance whose finite width is expected to be important for the extraction of the a 1 resonance properties. Indeed, a large part of the a 1 (1260) peak seen in the invariant mass distribution of three pions lays below the nominal ρπ threshold. However, the pole of the resonance was previously addressed in Lagrangian-based models [25,26], assuming a stable ρ-meson.
The J P C = 1 ++ three-pion state can be observed in the τ → 3π ν τ decay as well as in pion diffraction off a proton target π p → 3π p. There appears to be a discrepancy in the a 1 resonance parameters extracted from the two reactions [1,27]. The problem may be related to the presence of a large, coherent, non-resonant background, known as the Deck process in pion diffraction [28][29][30][31]. This process happens to dominate in the J P C = 1 ++ partial wave and directly influences the extraction of the a 1 (1260) resonance parameters in pion diffraction. Thus, an independent determination of the a 1 (1260) resonance properties is not only relevant for a better understanding of this state but also to constrain the Deck process, which contributes significantly to other partial waves including the ones with the exotic quantum numbers 1 −+ [31]. In this paper, we therefore focus on the τ − → π − π + π − ν τ decay with the aim of extracting the a 1 (1260) resonance parameters.
The paper is organized as follows. In Sec. II we present our model, we relate the differential width of τ − → π − π + π − ν τ to the three-pion scattering amplitude in the 1 ++ sector. In Sec. III we show how the model is constrained by the fit to ALEPH data. In Sec. IV we explore the analytic properties of our model for complex values of the 3π invariant mass squared, establishing the main singularities of the amplitude, and we determine the location of the a 1 (1260) pole. The studies of the systematics are described in Sec. V. Our conclusions are summarized in Sec. VI.

II. THE REACTION MODEL
We consider the reaction τ → 3π ν τ and derive an expression for the differential width which characterizes the 3π invariant mass spectrum [32][33][34][35]. The differential width is calculated by averaging (summing) over the τ (ν τ ) polarizations and integrating the matrix element squared over the final-state momenta, where m τ is the mass of the τ -lepton, m τ = 1776 MeV [1], the neutrino is considered massless, dΦ 4 is the four-body differential phase space, and λ x are the lepton helicities of the x = τ, ν. The process is dominated by the emission of a W boson by the leptonic current, 2 is the Cabibbo-favored weak coupling, p 3π , p τ , and p ν are the four-momenta of three-pion system and the leptons, u (ū) are the Dirac spinors of the τ (ν τ ), see Fig. 1. Because of G-parity conservation the π − π + π − system has positive C-parity. Hence, the vector currentūγ α u does not couple it, and can be removed. Since the W − is heavily off-shell, one should also consider the timelike polarization, which carries J P C = 0 −+ . However, the corresponding helicity amplitude is suppressed by the PCAC [33,36]. This enables us to treat the off-shell W − as purely axial. The polarization of the real W − provides a complete basis which we use to expand the hadronic current, where ε α * (Λ) 3π| is the helicity amplitude for the decay of the axial current to three pions. The squared matrix element summed and averaged over the ν τ and τ helicities, respectively, is The explicit evaluation of the expression is performed in the τ -rest frame where p τ · ε(0) = (m 2 τ − s)/(2 √ s), and p τ · ε(±) = 0.
Using the recursive relation for the phase space, we split it into the τ − → W − ν τ -phase space dΦ 2 , and the threepion phase space dΦ 3 : dΦ 4 = dΦ 2 dΦ 3 ds/(2π), where √ s is the invariant mass of the hadronic system. To obtain the differential width dΓ/ds, we integrate explicitly over the neutrino angles, Here, one power of the factor (m 2 τ − s) follows from the matrix element in Eq. (4), the other is given by the W − ν τ two-body phase space. The expression for the dΦ 3 is given in Appendix B. The integral is kept in the final expression to facilitate the further discussion on partial-wave expansion of the amplitude A Λ .
The helicity amplitude A Λ describes the coupling of the axial current to the three charged pions. The pions are labeled as follows, π − 1 π + 2 π − 3 (see Fig. 1). We use the isobar model to parametrize the dynamics and explicitly incorporate the π − 1 π − 3 Bose symmetry, 1: Diagram for the decay τ − → π − π + π − ντ . The momenta of the τ lepton and ντ are denoted by pτ and pν . The pions momenta are labeled by pi, i = 1, 2, 3. s is the invariant mass of the three pions.
where the isobar amplitude A (k) Λ includes only the subchannel interaction in a pion pair leaving the pion indexed k as a bachelor. In Eq. (6), we disregard the π − π − interaction since it is negligible compared to the dominant ρ-meson in the π + π − subchannel. The pion momenta are denoted by p i where i = 1, 2, 3 as shown in Fig. 1 and the subchannel invariant mass squared is denoted as σ k = (p i + p j ) 2 . Here and below we use the circular convention, i.e. the bachelor pion has index k such that the (ijk) are numbers (123), (231) or (312).
Each isobar amplitude receives different contributions, often referred to as decay channels [1]. The importance of different decay channels can be estimated by the relative branching fractions of the a 1 (1260) decay. The latest measurements were carried out by the CLEO experiment from τ decay [37,38] and by the COMPASS experiment in diffractive production [39]. The extraction of branching ratios is model dependent and is influenced by the production mechanism, however we get a rough estimate of their relative importance. The ρπ S-wave channel is dominant with a branching ratio of 60% − 80%. The second important channel, f 0 (500)π P -wave, was estimated to contribute less than 20%. The combined branching ratio to the remaining channels (ρπ D-wave, f 2 π P -wave, K * K S, D-waves) does not exceed 10%. We thus limit the analysis to the main ρπ S-wave channel. Including other decay channels would require the introduction of additional parameters for couplings and production strengths, which cannot be fixed by current publicly available data.
Therefore, we take the isobar amplitude to have the form, where C (k) = 1, µ i ; 1, µ j |1, 0 = ±1/ √ 2 is the Clebsch-Gordan coefficient relating the two pion with isospin projection µ i,j = ±1 to ρ 0 isospin states, thus, the sign depends on the index k. The a(s) denotes the dynamical part of the amplitude a 1 → ρπ S-wave in the canonical basis [40,41], f ρ (σ) is a parametrization for the ρ-meson decay amplitude, and N Λ (Ω k , Ω ij ) is the angular decay function for the decay chain a 1 → ρπ, ρ → ππ, The angles Ω k = (θ k , φ k ) are the polar and the azimuthal angles of the vector p i + p j in the 3π helicity frame, i.e. the center-of-mass (CM) frame with the axis orientation fixed by the production kinematics. The Ω ij = (θ ij , φ ij ) are the spherical angles of the pion i in the helicity frame of the isobar (ij). Detailed discussion on the decay function in Eq. (8) can be found in the Appendix B.
The line shape of the ρ-meson is given by the customary Breit-Wigner amplitude with dynamical width [39,42] where p(σ) = σ/4 − m 2 π is the pion break-up momentum, the function F 1 (pR) combines the threshold factor p(σ) and the customary Blatt-Weisskopf barrier factor with size parameter R = 5 GeV −1 . We use in the analysis m π = 139.57 MeV, m ρ = 775.26 MeV [1]. For convenience we fix N so that the phase-space integral ρ(s) defined below in Eq. (12a) approaches the two-body phase space asymptotic value, 1/8π, in the limit s → ∞, i.e.
The normalization for f ρ (σ) fixes the normalization of a(s).
Using Eqs. (6), (7) to substitute the amplitude A Λ in Eq. (5), we get the expression for the differential width in terms of the dynamic amplitude a(s). where ρ(s) is the effective ρπ phase space. We will consider two models for ρ(s)'s: The expression in Eq. (12a) strictly follows from Eqs. (6), (7), and (11). The label SYMM is introduced to emphasize the symmetrization between the decay channels, i.e. the ρπ channels k = 1 and 3. The relative minus sign comes from the symmetry of the isospin coefficient in Eq. (7). The integral in Eq. (12a) is the same for all helicities Λ due to the properties of the Wigner d-functions, therefore we set Λ = 0 for simplicity. The interference term is only significant at low energy, where the overlapping region of the two ρ-mesons contributes to a substantial fraction of the Dalitz plot. The ρ QTB (Quasi-Two-Body) in Eq. (12b) is a simplified phase space where the interference term is neglected. In this case, the integrals of the two decay chains squared are identical, which cancels the 1/2 factor in front. This model treats the ρ-meson as quasi-stable and the interaction between the ρπ as a two-body interaction. The simplification is suggested and discussed in Ref. [43] to treat the multiparticle final states. The same approximation is commonly used to account for 4π channel in the ππ/KK coupled-channels problem [44,45]). Finally, as shown in Fig. 2, the interference is rather small. Since this model is simpler, we would like to test it as an alternative. Our model for the decay amplitude is constrained by the approximate three-body unitary [29,43]. Turning general 3-body unitarity into some practically useful equations is cumbersome and not complete yet. A significant progress in this direction has been made in Refs. [14][15][16][17]. In particular, one can separate the genuine three-body unitary from the subchannel unitarity related to rescattering between different isobars. These processes modify the line shape of the subchannel amplitudes [46][47][48][49][50][51][52]. A good example is the ρπ-dynamics studied in the 1 −− sector in the decay of ω/φ [47,49], where the final-state interaction were found to shift and skew the ρ-meson peak. Conversely, in our models we focus on the 3-body resonance dynamics, and simplify the problem by neglecting the effects of the rescattering on the isobar line shapes. We introduce the ρπ elastic scattering isobar amplitude t(s), to impose the unitarity constraints for the amplitude a(s): where ρ(s) is the effective phase space given by Eq. (12a) or Eq. (12b). The factor of 2 in the left-hand-side of Eq. (13) is kept for convenience. The unitarity equations (13) can be satisfied by a certain choice of the parametrization. where C(s) is an analytic function constrained by condition Im iC(s) = ρ(s). To describe the amplitude dominated by a single resonance, we added a first order polynomial (m 2 − s)/g 2 to the denominator of t(s), which is equivalent to have the K-matrix with a single pole [41]. The numerator function α(s) is supposed to incorporate the singularities specific to the production process into the amplitude a(s). In the case at hand we use α = const. There are two common constructions for C(s) which both satisfy unitarity: 1. The models with C(s) = ρ(s) will be called non-dispersive. These models have left-hand singularities on the physical sheet inherited from the phase space, which are not motivated by physics.
where the subtraction constant l 0 is chosen such that the real part of iρ(s) is zero at the point (m ρ + m π ) 2 . The function iρ(s) has no singularities other than the unitarity cut as guaranteed by the Cauchy integral theorem. It is analogous to the Chew-Mandelstam function for the two-body scattering amplitude [29].
To summarize, the final expression for the differential cross section is.
Eq. (16) follows from Eq. (11). The constant c absorbs all energy-independent numerical factors; m, g, and c are real parameters which are fitted to data. The four models we are going to test are summarized in Table I. Our primary model is SYMM-DISP, which is the one that incorporates the most of physical arguments. The SYMM model contains additional left-hand singularities with respect to SYMM-DISP. The QTB and QTB-DISP models do not include the interference between the two decay chains, but are much simpler to calculate on the real axis and continue to the complex plane. The C(s) is calculated using the same ρ(s) as in the numerator of Eq. (16), which is either ρ QTB or ρ SYMM as given in Table I.

III. FIT RESULTS AND RESONANCE PARAMETERS
The largest public dataset for τ → 3π ν τ was collected by the ALEPH experiment in 2005 [53]. 1 The distribution dΓ/ds is binned in 0.025 GeV 2 bins and normalized by the measured branching ratio. We fit 103 data points in the range 0.38 GeV 2 ≤ s ≤ 2.94 GeV 2 . We minimize the χ 2 -function taking into account the covariance matrix provided in Ref. [53], where D is a vector of ALEPH data points, M (c, m, g) is a vector of the model predictions calculated for the centers of the bins. The matrix C stat is the covariance matrix of the statistical errors. The systematic uncertainties are smaller than the statistical ones by a factor 5, and can be neglected. Nonzero correlations among different bins are introduced by the unfolding procedure. It is worth noticing the 3π spectrum does not show the expected random noise. As discussed in the follow up analysis of the ALEPH [54], the problem appears because the errors of the unfolding procedure were not correctly propagated. Hence, the absolute value of χ 2 we obtained does not have a strict statistical meaning. However, we assume that for the model characterization based on relative χ 2 values, the problem should not be critical. The gradient minimization is performed using the NLopt optimizer and the ND MMA algorithm [55] with the automatic differentiation provided by the ForwardDiff.jl-package [56]. The minimum we find is always stable and isolated, as checked by repeating the minimization from different starting values. Fits to the ALEPH dataset are shown in Figs. 3, and the fit parameters and χ 2 values are shown in Table I. The non-dispersive models are not consistent with the data, with χ 2 at least three times worse than we have obtained for the dispersive models. In particular, they fail to reproduce the line shape around the peak and in the threshold region, and we do not consider them any further. On the other hand, the dispersive models show a good agreement with data, obtaining χ 2 /n.d.f. = 94/100 and χ 2 /n.d.f. = 61/100 for the SYMM-DISP and QTB-DISP, respectively.
In the next section we will perform the analytic continuation of the amplitude to the second sheet and search for the a 1 (1260) resonance pole. For comparison with the PDG [1], we first provide the customary Breit-Wigner parameters, that can be extracted on the real axis. We remind the reader that these are expected to be reaction-dependent, and do not provide an unambiguous characterization of the resonance. We define the Breit-Wigner mass squared m 2 BW as the value of s when the denominator of the amplitude t(s) in Eq. (14)

IV. ANALYTIC CONTINUATION THE POLE POSITION
Once the amplitude is fixed on the real axis, its analytic structure is unambiguously defined and can be explored. Unitarity introduces a branch cut along the real axis from the 3π threshold to infinity, which opens a non-trivial Riemann topology or sheet structure. The first Riemann sheet is the one containing the physical values of the amplitude slightly above the real axis. By construction, the amplitudes in the dispersive models contain no other singularity on the first sheet than the unitarity cut. Resonance poles are expected to lie on the second sheet, which is connected to the physical axis from below. The unitarity condition Eq. (13b) gives us a relation on the real axis that can be used to continue the amplitude in the complex s-plane. The real-axis relation followed from Eq. (13) reads is the discontinuity across the cut, s is real, is an infinitesimal positive number, and the Roman subscript indicates the Riemann sheet. Thus, t −1 II (s) = t −1 I (s) + iρ(s) and the pole positions are determined by t −1 II (s) = 0. The first sheet amplitude, t −1 I (s), is straightforward to calculate in the complex plane using the dispersive integral in Eq. (15). Continuation of the discontinuity, however, is more challenging since it is not explicitly analytical expression, as Eq. (12a) contains a modulus operator. Therefore, we need to find an analytic function which coincides with the discontinuity on the real axis. All singularities of the discontinuity −iρ(s) will be present in the second sheet amplitude according to Eq. (18). Among those, we expect the reflection of the ρπ unitarity cut, which is pushed into the second sheet due to the unstable nature of the ρ-meson.
For the continuation to the complex s-plane, we need to evaluate f ρ (σ) and f * ρ (σ) in Eq. (12a) and Eq. (12b) for complex argument σ. Along the physical axis f ρ (σ) = f A. Analytic continuation of the QTB-DISP model We start with the QTB-DISP model, whose analytic continuation is simpler than the one of the SYMM-DISP model. The discontinuity across the unitarity cut is given by −iρ QTB in Eq. (12b). The angular integrals in the phase space can be solved analytically due to the properties of the Wigner D-functions. We obtain where we used the definition λ i = λ(σ i , m 2 π , m 2 π ), λ si = λ(s, σ i , m 2 π ), with λ being the Källén triangle function λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + yz + zx). Using Eq. (19), we replaced |f ρ (σ 1 )| 2 by the analytic expression The corresponding ρπ cut locations are shown in the right panel of Fig. 4. The path C (hook) σ rotates the ρπ cut such that it opens up a larger area of the closest unphysical sheet and is used in the following for finding poles and illustration purposes.
The amplitude t(s) in the complex s-plane for the QTB-DISP model is shown in the left panel of Fig. 5. Naively, one would expect a single pole in the complex plane, originating from the single K-matrix pole, g 2 /(m 2 − s), present in Eq. (14). In contrast to this expectation, two poles are observed. Furthermore, both are rather close to the physical region. The correspondence between the K-matrix poles and the complex poles can be established by varying the coupling g. In the limit g → 0 the complex poles should approach the real axis at the position of the corresponding K-matrix poles. We find that the deep pole approaches the real axis at s = m 2 = (1223 MeV) 2 (see Table I with the fit results), while the left pole goes to s = 0. Due to these observations, we identify the deep pole with a 1 (1260)-pole label, i.e. corresponding to a resonance, and the left pole with a "spurious"-pole, i.e., an artifact from our parametrization in Eq. (14). This exercise also helps us to understand the origin of the spurious pole: it is the 1/s singularity in 2 The branch points are connected by cuts. Since the integral is calculated numerically it is important to make sure that the integration path does not cross any cut between the integration end points. To illustrate the cut choice shown in Fig. 4, we write For real values of s, this expression has two short branch cuts on the real axis: one between 0 and σ th , and the other between the points ( √ s ± mπ) 2 . When s is complex the first s-independent cut remains, while the second one splits into two straight cuts to the right with the branching points ( √ s ± mπ) 2 as shown in Fig. 4.
For the error estimation we used the bootstrap technique [57,58]: 1000 sets of pseudo data were generated using the original data and the covariance matrices, with the correlations taken into account in the Gaussian approximation. By refitting the pseudo datasets, we collect samples of the parameters, which we use to estimate their uncertainties. The distributions of the mass and width of the pole obtained from the bootstrap are Gaussian to a good approximation. The fit results and the calculated error ellipses are shown in Fig. 9. The mean values of the bootstrap sample for the pole positions differ from the real data fit results by < 0.2σ which indicate a good consistency and negligible bias of the bootstrap method [58].

B. Analytic continuation of the SYMM-DISP model
The evaluation of the discontinuity given by Eq. (12a) for complex s is more complicated since the angular integrals cannot be solved completely, see Appendix B. We start by casting ρ SYMM (s) in the form: where the first term in the sum is the phase-space factor in the QTB-DISP model, the second term is the interference contribution given by Eq. (B16). Substituting f ρ → f .
The function W (a, b, c) is a multivariable polynomial defined in Eq. (B14). Omitting constant factors, the function f ρ (σ) is given by A right-hand cut is introduced by i 4m 2 π − σ. In addition, there are two branch points: one at σ = 0 from the phase space in the width Γ(σ), and another one at σ = 4m 2 π − 4/R 2 due to the Blatt-Weisskopf factor in the numerator. The break-up momentum singularity σ − 4m 2 π in the numerator of f (σ) is canceled by the same factor which arises from the angular function (see Eq. (24)). The parametrization of f ρ (s) in Eq. (25) contains 5 poles, as one can count by the order of the polynomial which would give zeros of the denominator. They correspond to the ρ-meson poles at (m ρ ± iΓ ρ /2) 2 , and three spurious poles lying far away from the physical region as shown in Fig. 4. The integration endpoints of the σ 3 variable, σ ± 3 (σ 1 , s), describe the border of the Dalitz plot for fixed value of s (Fig. 6, left panel), As soon as s becomes complex the endpoints depart from the real axis and move into the complex plane. The trajectories of the σ ± 3 as functions of σ 1 moving from 4m 2 π to ( √ s − m π ) 2 are non-trivial. As shown in Fig. 6, while σ 1 moves along the C (hook) path (see Eq. (21c)), the σ − 3 circles around the branch point 4m 2 π . When σ 3 crosses the unitarity cut, the sheet, on which it is evaluated, must be changed. However, if the σ 1 path goes exactly through the point (s − m 2 π )/2, σ − 3 just touches the branch point 4m 2 π , (indeed, σ − 3 ((s − m 2 π )/2, s) = 4m 2 π ). In that case we are allowed to stay on the same sheet. Therefore, there are two ways to calculate ρ INT for a complex argument (see Appendix C for more details): 1. ρ (1) INT : We choose a special path in σ 1 , the σ ± 3 always stay on the same sheet and can be connected with a straight (undistorted) path.

ρ
(2) INT : We let σ − 3 circle around the branch point, changing sheets of f (σ 3 ) appropriately. When σ 1 = σ th , the integration limits σ ± 3 coincide. For certain values of σ 1 , σ − 3 changes the sheet and, therefore, when σ 1 is in its upper limit σ lim , the positions of σ ± 3 coincide, but they are on the different sheets. The integration path must be detoured around the branch point as shown in Fig. 7.
The first option provides a unique continuation of Eq. (B16), however, the integration contour is bound to pass through (s − m 2 π )/2 which is non-analytic point of the integrand (see Appendix C). The integrand in the second option stays analytic on the integration contour, but in the limit of real s, the function ρ  of Fig. 6. .
The difference is practically negligible as shown in Fig. 2. The impact on the fit parameters and the values of the amplitude in the complex plane is a few orders of magnitude smaller than the statistical uncertainties. For the following discussion we use ρ INT (s) for the reason that the ρπ-cut can be rotated in the same way as before by using C (hook) path in σ 1 . Interestingly, an analogous problem appears in relation to the Khuri-Treiman equations (see Appendix in Ref. [14], Sec. IV in Ref. [59]). Ref. [60] gives arguments in favor of the first option.
As soon as the discontinuity is known for the whole complex plane, the amplitude on the unphysical sheet can be computed according to Eq. (18). The contour plot on the right panel of Fig. 5 presents the closest unphysical sheet of the amplitude, which is smoothly connected to the real axis. We find two poles and identify them as the a 1 (1260) resonance pole and the left "spurious" pole as shown in The statistical errors are obtained from a bootstrap analysis as described above in Sec. IV A. The combined results are presented in Fig. 9.

V. SYSTEMATIC UNCERTAINTIES
The description of three-particle resonances is a difficult problem because of the complicated structure of finalstate interactions, which induces an interplay between different decay channels. The latter manifests itself in the modification of the isobar line shape and the presence of interference terms. The importance of three-body effects is readily seen in the difference of SYMM-DISP and QTB-DISP pole positions, cf. Eqs. (29) and (22). Knowing that the interference between two ρπ decay channels must be present, we now focus on systematic studies of SYMM-DISP, keeping QTB-DISP for a mere comparison. The largest systematic uncertainty is the dependence of the a 1 (1260) pole position on the line shape of the subchannel resonance ρ. In principle, we know that final-state interactions shift and skew the ρ peak. The scale of the ρ-meson mass shift can be estimated from the studies of ω/φ decays using Khuri-Treiman equations [47,49]. Fig. 3 of Ref. [47] suggests a shift of the real and imaginary parts of the isobar amplitude of the order of 10 MeV before and after final-state interactions are taken into account. To estimate the effect on the a 1 (1260) pole position, we vary the parameters of f ρ (σ) in Eq. (9), i.e. the mass m ρ , the width Γ ρ and the Blatt-Weisskopf radius R, performing a χ 2 scan over each parameter, while keeping the others at their nominal values (Fig. 8). The new pole position obtained for the parameter value which minimizes the χ 2 for each scan is then used to estimate the systematic error for the pole position of the main fit. The results of these studies are summarized in Table II (see fit studies #2-7, were #4 was introduced as an additional intermediate point outside of the minimum). The a 1 (1260) pole position is extracted, the results for the pole mass and width are represented in Fig. 9 by open ellipses. The values m, g and χ 2 for fits described in Sec. V. For scans over parameters mρ, R and Γp we present the values of m, g and χ 2 obtained in the the minimum in the profile χ 2 plots shown in Fig. 8. We perform an additional test of the influence of heavier resonances, as the a 1 (1640), by excluding the region s > 2 GeV 2 from the fit. The fit quality does not change substantially, but get slightly worse due to the reduction of the degrees of freedom (see #1 in Table II). The values for the pole position are shown in Fig. 9 and included to the systematic error of our final result.

VI. CONCLUSIONS
In this paper we have presented a new analysis of the lightest iso-vector axial-vector resonance a 1 (1260) decaying to three charged pions. Despite the fact that the corresponding J P C = 1 ++ partial wave dominates the hadronic weak decay of τ leptons as well as diffractive reactions of high-energy pions, the parameters of the a 1 (1260) are still poorly known. While the latter reactions suffer from an irreducible background due to non-resonant processes, the system of three pions produced in τ decay provides a very clean access to axial-vector resonances. Compared to a twoparticle system, however, the system of three interacting particles exhibits additional phenomena, such as 3-particle rescattering or interference between different decay chains. These 3-body effects are taken into account using reaction models constraining the dynamics in the total invariant mass, however, without imposing subchannel unitarity. We have considered four analytic models of an isolated resonance decaying to three pions via the ρπ channel. All these models satisfied approximate three-body unitary, but differ by the left-hand singularities and the treatment of the interference between the two ρπ decay channels. Using the τ − → π − π + π − ν τ data from ALEPH [53], we found that the dispersive models, having no left-hand singularities on the physical sheet, fit the data clearly better.
In order to find the pole position corresponding to the a 1 (1260) resonance, we have explored the analytic structure of the amplitude and performed its analytic continuation into the complex plane of the three-pion invariant mass squared, a challenging, and technically demanding task, requiring us to use a prescription for the integration paths in the two-pion invariant mass squared. We have searched for the singularities in the closest unphysical sheet, and have identified a pole as the a 1 (1260) resonance. The mass and width of the a 1 (1260) are given in terms of its pole position in the main SYMM-DISP model: The dominant source of systematic errors is the sensitivity to the details of the subchannel interactions. The simplified QTB-DISP model, which neglects the interference between the two ρπ-channels, results in a significantly different pole position and a larger systematic uncertainty. This analysis can be extended by further advancing the theoretical framework and constraining the model by fitting the Dalitz decay variables. This will be possible when the data from BelleII or BESIII become available. In addition, the results from this analysis will help to better constrain the non-resonant background in diffractive reactions, as measured by the COMPASS experiment. Performing the analytical continuation in Sec. IV we have shown that, in addition to the expected a 1 (1260) pole, there is a spurious pole rather close to the physical region. At first, the spurious pole looks surprising, however, it is clearly present in every Breit-Wigner-like model of a resonance decaying to particles of different masses. Indeed, the denominator of the Breit-Wigner amplitude with energy-dependent width decaying to two scalar particles in an S-wave reads: When m 1 = m 2 , the equation D BW (s) = 0 has 4 complex roots, which we can identify by the order of the polynomial which gives those roots: Since all coefficients of the polynomial are real, the poles appear in conjugated pairs above and below the real axis. The two Breit-Wigner poles below the real axis are analogous to the a 1 (1260) and the spurious pole. To demonstrate this further, we draw the complex plane of the 1/D BW (s) function with m = 1.2 GeV, g = 7.8 GeV, m 1 = m ρ , m 1 = m π in Fig. 10. We find that the spurious pole has no influence on the physical region as long as the resonance is far from threshold and rather narrow. Both poles become important for the real axis physics when the studied resonance is close to threshold or/and wider. The spurious pole is a feature of Breit-Wigner-like models. It is generated by the 1/s singularity of the phase space in Eq. (A1), and Eq. (20). In order to remove it, we try to exclude the 1/s factor from the dispersive term. Following the studies of QTB-DISP, we consider a new model for scattering and production amplitudest(s) = t(s)/s andâ(s) = a(s)/s, and modify the unitarity equations accordingly.
2 Imt(s) =t * (s) (sρ QTB (s))t(s), 2 Imâ(s) =t * (s) (sρ QTB (s))â(s), where sρ(s) is free of the 1/s singularity. The parametrization which satisfies the unitarity constraints iŝ where the index k gives the number of parameters in the function K −1 k (s), the models are labeled sQTB-DISP (k) . The function sρ(s) has a ∼ s 1 asymptotic behavior, therefore the dispersive integral must be subtracted twice. The integrand is thus the same as in Eq. (15), but the integral is multiplied by an extra factor of s as in Eq. (A4). To make the dispersive integral independent of the subtraction points we must consider a polynomial of order k ≥ 2. We  consider three forms of functions K k (s), The K 2 (s) and K 4 (s) are inspired by the K-matrix approach with one and two poles, respectively, while K 3 (s) is a special two-pole model which exactly coincides with QTB-DISP when h = 0. In Fig. 11 we show the continuation of the sQTB-DISP (2) model, fitted to data. The spurious pole is no longer present. However, the quality of the fit is not acceptable: the best χ 2 /n.d.f. is equal to 979/100. When we increase the freedom by taking the model sQTB-DISP (3) the fit quality significantly improves to yield a χ 2 /n.d.f. = 67/100. Quite spectacularly, the picture of the complex plane is changed back: the place of the spurious pole is taken by the explicit pole introduced in the K-function (see the right plot of Fig. 11). The next relaxation of the setup in sQTB-DISP (4) overfits the data and gives χ 2 /n.d.f. = 42/100. However, the positions of the resonance and spurious poles do not change much.
The position of the spurious pole was investigated for all systematic studies we performed in Sec. V as shown in Fig. 12. Appendix B: Analytical simplification of the phase-space integral In this Appendix we demonstrate how the integrals in the phase-space factor ρ(s) Eq. (12a) can be simplified using the properties of the Wigner D-functions.  Table II. We start by explicitly defining the angles in the decay functions N 0 (Ω k , Ω ij ) given by Eq. (8). The three-pion center-of-mass (CM) frame is oriented by the direction of W in τ decay (W helicity frame). The momentum vector of the τ defines the xz plane, a.k.a. the production plane. Ω k = (θ k , φ k ) denotes the polar and azimuthal angles of the vector p i + p j in the CM-frame. The Ω ij = (θ ij , φ ij ) are the spherical angles of the pion i in the helicity frame of the isobar (ij). This helicity frame is obtained from the CM frame by active rotation R −1 (Ω k ) and boost along the z-axis. Equivalently, we can notice that the boost does not change azimuthal orientation, therefore, the y-axis direction e y in the helicity frame can be found by e z × e z , where e z is the original orientation of the CM z-axis.
We write the phase-space differential through the two pairs of spherical angles.
(B17)  1. X (1) : we draw the σ 1 path directly through the branch point (s − m 2 π )/2 (the point P in Fig. 13 can be aligned with the branch point (s − m 2 π )/2). The point σ 1 = (s − m 2 π )/2 is special because when the path goes through it, the σ − 3 does not circle the branch point but just touches it.
2. X (2) : we go analytically under the cut taking any arbitrary path. X (2) corresponds to the function which we would obtain connecting the points σ ± 3 properly, i.e. avoiding 1/ σ 3 − 4m 2 π cut. The two options give two different analytical functions. Additional discussions on the subject can be found in Ref. [14,59].