A mathematical solution for the parameters of three interfering resonances *

: The multiple-solution problem in determining the parameters of three interfering resonances from a ﬂt to an experimentally measured distribution is considered from a mathematical viewpoint. It is shown that there are four numerical solutions for a ﬂt with three coherent Breit-Wigner functions. Although explicit analytical formulae cannot be derived in this case, we provide some constraint equations between the four solutions. For the cases of nonrelativistic and relativistic Breit-Wigner forms of amplitude functions, a numerical method is provided to derive the other solutions from that already obtained, based on the obtained constraint equations. In real experimental measurements with more complicated amplitude forms similar to Breit-Wigner functions, the same method can be deduced and performed to get numerical solutions. The good agreement between the solutions found using this mathematical method and those directly from the ﬂt veriﬂes the correctness of the constraint equations and mathematical methodology used.


Introduction
One of the main aims during the physics analysis of experimental data is determination of the parameters of several resonances by fitting the cross sections or measured mass spectrum with possible interference between the resonances considered. In some cases, although the fitted results with interference are not taken as nominal results, the interference still needs to be considered as an estimate of the systematic uncertainty.
In particle physics, we usually take Breit-Wigner (BW) function to represent resonance amplitude. A typical task is determination of the BW parameters from the fit to the measured distributions in experiment, such as cross sections. The measured physical quantities are usually in proportion to the modulus of the total amplitude squared, for example, |BW 1 +BW 2 e iφ | 2 for two interfering resonances and |BW 1 +BW 2 e iφ 1 +BW 3 e iφ 2 | 2 for three interfering resonances, where φ, φ 1 , and φ 2 are the relative phases between resonances. Due to this squaring operation in the amplitudes, to connect with the measured physical quantities, we can find multiple solutions in extracting amplitudes from the fit to the experimental measurements. Often it occurs that these multi-solutions have the same goodness-of-the-fit, and resonance mass and width, but the relative phases are different. This indicates that different solutions have different coupling strengths to decay channels, which would result in different interpretations in physics. Therefore, for a fit with interfering resonances, we need to make sure that all the solutions have been found. If there are multiple solutions, but only one is reported, the experimental results may be incomplete or even biased.
Recently, more and more experimental analyses, especially the studies of the vector charmonium-like Y states, have indicated this. For example, in Ref. [1] two or three coherent resonances plus an incoherent background shape are used to fit the π + π − ψ(2S) invariant mass distribution. Correspondingly, two or four solutions are found with identical resonance mass and width but different couplings to electron-positron pairs. Another example is presented in Ref. [2], where two solutions are found in the branching fraction measurement for the φ→ωπ 0 process and the study of ρ−ω mixing.
In real physics analyses, all the multiple solutions are found via a fitting process. Due to background statistical fluctuations or limited statistics, not all the solutions can be found easily in some cases. Therefore, from the mathematical point of view, a natural question is raised: if a particular solution has been found, whether other solutions can be derived from it. For the above question, the authors in Refs. [3,4] proved that if we use two coherent BW functions to fit the measured distribution, there should be only two different solutions, and they can be derived from each other using analytical formulae and a numerical method. As pointed out in Ref. [4], in the case of three resonances with constant widths there are four solutions with the same likelihood function minimum, but an analytical solution of this problem appeared too hard due to technical difficulties.
In this paper, we discuss the multiple-solution problem in determining the resonant parameters of three interfering resonances from a mathematical viewpoint. Although explicit analytical formulae cannot be derived, we provide some constraint equations between the four solutions. We also provide a mathematical method to get additional solutions from the obtained one.
This work is organized as follows. After the Introduction, we present a general mathematical model for the amplitudes of three coherent resonance states in Section 2. If the three resonances are described by the normal BW functions, the analytical expressions for the relationship between the four solutions are deduced and obtained. An effective approach is developed to obtain the algebraic equations of the relationships between the four solutions. In Section 3, the relations between the four solutions are also deduced for relativistic BW forms. In Section 4, two numerical examples produced by toy Monte Carlo (MC) are utilized to cross check and confirm our results. When the form of the resonance amplitude is extremely complex, we can take a similar numerical procedure to obtain other unknown solutions from the known one. Finally, in Section 5, a short discussion is given.

Mathematical methodology for three simple BW amplitudes
In the light of two distinct features, (1) all solutions have the same goodness-of-fit and (2) different solutions have identical resonance masses and widths but different couplings to electron-positron pairs, we construct a general mathematical model for multiple solutions based on three interfering amplitude functions.
A sum of three quantum amplitudes can be described by a complex function e(x,z 1 ,z 2 ,z 3 ) with the form where x is a measured variable, g(x), h(x), and f (x) are complex functions of x, and z 1 , z 2 , and z 3 are complex numbers. Our purpose is to find different parameters z 1 , z 2 , and z 3 satisfying Since the global phase does not work on amplitude squared operations, we can reduce the dimension of {z 1 ,z 2 ,z 3 } parameter space to a {d,z α ,z β } parameter space, where d is a real number. The module of the amplitude squared of e(x,z 1 ,z 2 ,z 3 ), |e(x,z 1 ,z 2 ,z 3 )| 2 , can be rewritten in a more convenient form by defining Here . Considering |g(x)| 2 is only a product factor and is independent of z α , z β , and d, we remove it in the following derivation. What we need to do now is to find different z α , z β , and d values , and (R z β , I z β ) as real and imaginary parts of F (x), H(x), z α , and z β , respectively, and using them to represent For the sake of brevity, the specific form of dependence of R F (x), I F (x), R H (x), and I H (x) on x is removed here. Without loss of generality, we take d = 1 as an initial solution for convenience. The next task is to find all the possible z α , z β , and d values to make To be more specific about our work, we consider that g(x), h(x), and f (x) are widely accepted nonrelativistic BW functions, as an example: where M is the mass and Γ is the width of a resonance, respectively. Using the above forms of g(x), h(x), and f (x), the real and imaginary parts of F (x) and H(x) become , 043001-2 , respectively. After some algebra, we obtain the interesting relations below: A similar expression can be obtained for E(x,z α ,z β ).
] is a constant for x space. We notice that the term (R F R H +I F I H ) and the linear combination of R F , I F , R H , and I H have the same number of x terms with the same power. It is the same for the term ( The factors {c 1 ,c 2 ,c 3 ,c 4 ,c 5 } and {c 6 ,c 7 ,c 8 ,c 9 ,c 10 } follow Eq. (11): .
Then we can get with A=2(R zα R z β +I zα I z β ) and B =−2(R zα I z β −I zα R z β ). We know that R zα , I zα , R z β , and I z β are functions in parameter space {d,z α ,z β }. If we want to make E(x,z α ,z β )/d = E(x,z α ,z β ) hold for any x, then the corresponding coefficients of the functions in parameter space should be equal, which immediately leads to the following equations: All what we need is to solve Eq. (13) to obtain the values of R zα , I zα , R z β , I z β , and d . Unfortunately, there are no explicit analytical expressions for them. So we cannot prove that there must be four solutions. Such a conclusion agrees with that in Ref. [4]. However, by using Mathematica [5] to input Eq. (13) and an initial solution, we get exactly four numerical solutions quickly. The numerical solutions can be taken as cross checks and references compared with those from the fits. This definitely saves a lot of time and energy. We need to point out that Eqs.

Mathematical Methodology for three relativistic BW amplitudes
Here we take another form for f (x) , g(x), and h(x), i.e., relativistic BW amplitudes that are usually used in e + e − reactions to extract the parameters of the Y resonance: where s is the e + e − center-of-mass square; M R is the mass of the resonance R; Γ R and Γ R e + e − are the total width and partial width to e + e − , respectively; B R is the branching fraction of the resonance R decays into a final state; and P S is the n−body decay phase space factor, which increases smoothly from the mass threshold with √ s [6]. Notice that Eq. (13) is independent of the forms of amplitudes, while its coefficients will change. With some algebra, we can obtain the coefficients for other forms of amplitudes.
With Eq. (14), the F (x) and H(x) are changed to In this situation, R F , I F , R H , and I H are changed. So we need resolve parameters a f , b f , c f , a h , b h , c h , {c 1 ,c 2 ,c 3 ,c 4 ,c 5 }, and {c 6 ,c 7 ,c 8 ,c 9 ,c 10 } using Eqs. (6) and (10), respectively. We obtain and , 043001-4 c 3 =

Simple BW amplitudes
In order to verify our deduction on constraint equations and mathematical program in obtaining numerical solutions, let us take a random example for the case of three simple BW amplitudes with interference. The parameter values of the three BW functions as one solution are set as The modulus of the amplitude squared of three interfering resonances is |BW g (m)+BW f (m)e iφ f +BW h (m)e iφ h | 2 and the BW amplitudes use the formats shown in Eq. (5). That is to say, z α = e iφ f = 1/2 + √ 3/2i and z β = e iφ h = −1/ √ 2+1/ √ 2i for the above solution. Using the above probability density function and the first set of input solutions, toy MC is used to generate a data sample of 100,000 events. The generated distributions with dots with error bars are shown in Fig. 1. A binned extended maximum likelihood fit is applied to such a distribution with three interfering resonances to extract the parameters of the resonances. Four sets of solutions are found. The fitted results are summarized in Table 1 and the corresponding fitted plots are shown in Fig. 1 in solid lines. Using the aforementioned method, we can also obtain another three sets of solutions numerically. We found the numerical solutions are exactly repeated by fitting. For those with little difference, they are consistent within 0.5σ, where σ is the error from the fit. A comparison of the results is shown in Table 1.
It is obvious that, for the case of three nonrelativistic BW amplitudes with interference, if one solution is known from the fit, the other three can be derived readily and numerically by solving Eq. (13).

Relativistic BW amplitudes
For the case of relativistic BW amplitudes with interference, the values of the parameters as one solution are set as The modulus of the amplitude squared of three interfering resonances is |BW g (m)+BW f (m)e iφ f +BW h (m)e iφ h | 2 and the BW amplitudes use the formats shown in Eq. (14), where for the phase space factor we assume the reaction process is e + e − → π + π − J/ψ. That is to say z α = B f Γ f e + e − /B g Γ g e + e − e iφ f = i and z β = B h Γ h e + e − /B g Γ g e + e − e iφ h = −1/ √ 2+1/ √ 2i for the above solution, where the values of B R Γ R e + e − are set as 1 for R=g, f, and h.
According to the above probability density function and the first set of input solutions, a data sample of 100000 events is generated by using toy MC. Similarly,     using the method mentioned earlier, another three sets of solutions can be found numerically, which are exactly repeated by fitting with the maximum likelihood method. The comparison of the results is shown in Table 2 and the fitted plots are shown in Fig. 2.

Discussion
As we found, when we need to describe a measured distribution using three interfering resonances |g(x) + z α f (x) + z β h(x)| 2 /d , F (x) = f (x)/g(x) and H(x) = h(x)/g(x) satisfy the relation of Eq. (6). If f (x), h(x), and g(x) are widely used BW functions, it has also been proved that such relation is exactly satisfied. In the case of three interfering resonances there are already four equivalent solutions with the same likelihood function minimum. Although the explicit analytical formulae cannot be derived between different solutions, Eq. (13) can be utilized to derive the other three solutions numerically from the solution obtained by fitting. If three resonant amplitudes take simple or relativistic BW functions, two data samples generated by toy MC are used to cross check and verify our results. For other complicated BW functions, the relations Eqs. (6), (10), (12), and (13) still hold for F (x) and H(x). For other forms of BW functions, with the coefficients obtained by Eqs. (6) and (10), the other solutions can be derived numerically by using the method mentioned earlier. The obtained numerical solutions agree well with those from the fit, which justifies our method. We believe that with the help of finding other solutions numerically, it is easy to find all the solutions in real fits to the experimental distribution as long as the initial values of resonant parameters are set correctly.