Methodology to determine the spin-parity of muon-philic $X$ boson in $J/\psi \rightarrow \mu^- \mu^+ X$ decay

The present anomaly in muon anaomalous magnetic moment can be explained by the presence of a muon-philic $X$ boson which could be a scalar particle or a vector particle with mass less than twice the mass of muon. The muon-philic $X$ boson could interestingly not be a parity eigenstate as well. If there exists such a boson, irrespective of its parity, it can be directly observed in the decay $J/\psi \to \mu^- \mu^+ X$ where $X$ remains invisible. We show that by using the angular distribution or the distribution of events in the square Dalitz plot, along with two well defined dimensionless ratios, one can clearly distinguish among the various spin-parity possibilities. This would constitute an important probe of both the existence and the nature of this new physics possibility.


I. INTRODUCTION
It is well known that the anomalous magnetic moment of muon, a µ = g µ − 2 /2 where g µ is the g-factor of muon, is sensitive to contributions from various new physics possibilities via quantum loops (see Ref. [1] for a recent review on contributions from various specific new physics possibilities on a µ ). Therefore, a difference between the experimental measurement of a µ and its standard model (SM) prediction is considered to be a good probe of new physics. The precise experimental measurement of muon anomalous magnetic moment as reported by the E821 experiment a exp µ at Brookhaven National Laboratory (BNL) [2] is found to be larger than the existing most precise theoretical prediction for the same from the SM a SM µ [3]: This 3.7σ discrepancy can be considered as a tantalizing hint of the presence of some new physics. It is expected that this discrepancy might increase in the near future once results from the E989 experiment at Fermilab [4] are available which would improve the existing experimental measurement by a factor of four. The J-PARC New g-2/EDM experiment at KEK, Japan [5] is also going to probe this discrepancy with more precision. If the discrepancy goes beyond the accepted discovery threshold of 5σ it would only imply that there is definitively some new physics contributing via the quantum loops. However, to pin point the nature of the new physics it is pertinent that we directly probe the various new physics possibilities in other experiments. In one of the simplest new physics possibilities one introduces a new, electrically neutral, spin-0 or spin-1 boson, say X, with exclusively muon-philic interaction so as to avoid possible constraints from various other experimental studies. The X boson can have positive or negative parity, i.e. it could be a scalar, pseudo-scalar, vector or axial-vector particle, or it could even not be a parity eigenstate. To denote these possibilities in a unified manner let us write X ≡ X s ± * Email at: manimala@iopb.res.in † Email at: dibyakrupa.s@iopb.res.in where s = 0, 1 denotes the spin and ± signify the parity. It is well known that for a scalar (X 0 + ), pseudo-scalar (X 0 − ), vector (X 1 + ) or axial-vector (X 1 − ) boson coupling exclusively to muon via the following interaction Lagrangians, with mass m X 2m µ where m µ is the mass of muon (so that X is invisible), we get the following contributions to ∆a µ at leading-order [1] from Fig. 1, where g s± are the various coupling constants introduced in Eq. (2). It is clear that when X is not an eigenstate of parity, it simply implies that we could have both scalar and pseudoscalar couplings for spin-0 case and both vector and axialvector couplings for spin-1 case. In such cases, one would have to add the corresponding contributions to ∆a µ to explain the observed anomaly in Eq. (1). There exist possible UV complete models that can provide interactions as required in our Eq. (2), see Ref. [1]. Inevitably, these UV complete models introduce additional interactions which put constraint on the relevant parameter space from other processes. For the scalar case, one example is the type-X two Higgs doublet model [6][7][8]. As was recently shown in [9], a light scalar can accommodate both muon g − 2 anomaly as well as the KOTO anomaly in K L → π 0 νν decay [10,11]. Similarly, for the vector case, an example of UV complete model could be the Contribution from the muon-philic boson X to ∆a µ . model which conserves the difference between muon and tau lepton number, in addition to being anomaly free [12][13][14].
In a recent study [15] it was shown that such a model could accommodate both muon g − 2 anomaly and the anomaly reported by the Atomki experiment in 8 Be * [16,17]. The mass of the new scalar or vector particles is in the MeV region which is also the region of interest to us. However, in this paper we are not concerned with the UV completion of the simplistic muon-philic model that we have considered. Instead we will focus on how the existence of the X boson can be studied and how its spin-parity can be determined. This would offer a direct probe of the nature of the new physics under our consideration. It is clear from Eq. (3) that the pseudo-scalar and axialvector contributions have opposite sign compared with contributions from scalar and vector cases. This would imply that pseudo-scalar and axial-vector contributions would push the theoretical value farther from the experimental value, and can not explain the observed anomaly of Eq. (1) when considered separately. Therefore, purely from the point of view of muon anomalous magnetic moment one usually considers only the scalar and vector new physics scenarios. In a recent work of one of the authors [18] it was explored how the decay J/ψ → µ − µ + X can indeed be used to probe this sort of new physics in a conclusive and independent manner at the BES III experiment. This decay mode can be compared to the analogous process e + e − → µ − µ + X proposed in context of Belle II [19] which has a large number of competing background processes that must be carefully considered. The X boson in the final state is electrically neutral and does not decay (in the simplest muon-philic scenario as exemplified by Eq. (2) above) leaving no track in the detector. As explained in detail in [18] the most dominant SM background for the decay J/ψ → µ − µ + X arises from final state radiation which can be readily handled experimentally by employing a cut on the missing energy distribution. The feasibility of observing the decay process is exemplified by the fact that around 30 to 300 events for the scalar case and about 300 to 2000 events for the vector case are expected at BES III after the aforementioned missing energy cut is applied [18].
It is important to note that our paper differs significantly from Ref. [18] by the facts that (1) we consider all the four spin-parity possibilities of X as well as the possibility of X not being a parity eigenstate, instead of only scalar and vector cases as in [18], (2) we show that for the case of X being not a parity eigenstate the parameter space allowed by muon g−2 is widened enhancing the branching ratios for the decay J/ψ → µ − µ + X making their experimental observation more feasible, and (3) we provide a completely different methodology to ascertain the various spin-parity possibilities of X using the difference in distribution patterns in conventional Dalitz plot, angular distribution, square Dalitz plot as well as by using two experimentally accessible dimensionless ratios.
Our paper is organized as follows. In Sec. II we find out the region of parameter space allowed by the existing anomaly in muon anomalous magnetic moment, considering the simple muon-philic interactions of Eq. (2). This is followed by an analytical analysis of the decay J/ψ → µ − µ + X for all the various spin-parity possibilities of X in Sec. III. Along the way, in Sec. IV, we numerically show the differences in patterns of distribution of events in different kinds of Dalitz plots and angular distributions for the various spinparity possibilities, and discuss how they can be used in conjunction with two dimensionless experimentally measurable ratios to distinguish the various spin-parity possibilities. Finally we conclude in Sec. V highlighting the main ideas and the future prospect of our study. g 0+ from ∆a µ (2σ allowed) g 1+ from ∆a µ (2σ allowed) FIG. 2. The coupling constants g 0+ and g 1+ as allowed by ∆a µ at 2σ level for m X 2m µ . Here X is considered to be a parity eigenstate, either a scalar or vector.
By ascribing the anomaly in muon anomalous magnetic moment to muon-philic interactions alone we can readily find out the allowed range of values for the coupling constants g s± where s = 0, 1. Allowing the current discrepancy in Eq. (1) at 2σ level and considering Eq. (3) and taking the X boson as a particle with definite parity, the scalar coupling constant g 0+ and the vector coupling constant g 1+ are allowed to vary in the green colored bands of Fig. 2.
If, instead, we relax the condition that the X boson must A part of the allowed parameter space for g 1+ and g 1− .
FIG. 3. The region of parameter space allowed by ∆a µ at 2σ level when X is not a parity eigenstate, shown for a few chosen values of m X 2m µ . The allowed region for the specific m X value is the region bounded by the corresponding colored and dashed lines. be a parity eigenstate, the coupling constants g 0+ and g 0− are together allowed for spin-0 case, and for spin-1 case both g 1+ and g 1− are allowed. Because of the fact that we have two parameters to fit one data, the allowed region naturally widens, as can be seen from Fig. 3. Interestingly, it is the negative sign in front of Eqs. (3b) and (3d) which enable this widening in parameter space. In Fig. 3 the region allowed by ∆a µ at 2σ level is shown for the specific m X values of 10 MeV, 110 MeV and 210 MeV as the region bounded by the corresponding colored and dashed lines. It is interesting that one can consider much larger values than what is allowed if X were to be a parity eigenstate. It is also clear from Fig. 3 that the special cases of g 0+ = g 0− are also allowed for most of the m X values.
Since the branching ratio of any process that would probe X would be, in general and at least, proportional to square of the coupling constant, a larger value of the coupling constant would imply bigger branching ratio. Therefore, the scenario where X is not a parity eigenstate could be probed better in experiment. Nevertheless, as we have seen from Figs. 2 and 3 the study of muon anomalous magnetic moment alone can not decipher the nature of the new physics. In order to find out the nature of new physics we need to supplement exploration of ∆a µ with some other study, such as the decay J/ψ → µ − µ + X. By analyzing the Dalitz plot distributions of J/ψ → µ − µ + X or its angular distribution in the center-ofmomentum frame of µ − µ + one should, in principle, distinguish between the various new physics possibilities. In the following we develop this methodology in detail.

III. STUDY OF THE DECAY
The decay J/ψ → µ − µ + X s ± takes place via the Feynman diagrams as shown in Fig. 4. The electrically neutral X boson can arise from either of the muon legs, and considering only the muon-philic interactions of Eq. (3) it is stable and remains invisible in the mass regime m X 2m µ which is the region of our interest. Considering the X boson to be a parity eigenstate, it was shown in Ref. [18] that this decay can be studied at BES III. As noted in the previous section, from Fig. 3 we find that when X is not a parity eigenstate the coupling constants can have much larger values and, therefore, it is evident that this possibility is experimentally interesting.
Since we are considering a three-body decay, the allowed phase-space is fully described by only two variables, such as two energies, or two invariant mass-squares, or one invariant mass square and an angle. Below we consider a few of these pairs of variables and find out the expressions for corresponding distributions.

A. The normal Dalitz plot distribution in terms of two invariant mass-squares
Denoting the 4-momenta of J/ψ, µ − , µ + and X by p J , p − , p + and p X respectively, let us define the following three invariant mass-squares s, t and u, The three invariant mass-squares are not independent since s+t+u = m 2 J +m 2 X +2m 2 µ , where m J is the mass of J/ψ. Taking the invariant mass-squares t and u as the two parameters, the Scalar (X 0 + ) and Pseudo-scalar (X 0 − ) cases Vector (X 1 + ) and axial-vector (X 1 − ) cases 4. Feynman diagrams for the decay J/ψ → µ − µ + X. The notation for 4-momenta of all the particles is also shown here. differential decay rate for the decay J/ψ → µ − µ + X s ± can be expressed by where α is the fine-structure constant, f J = 0.407 GeV [20] is the decay constant of the J/ψ, and |A s ± | 2 are some function of t, u and the masses involved, as given in Appendix A. All differences among the scalar, pseudo-scalar, vector and axial-vector possibilities are present in the expressions for |A s ± | 2 as shown in Eq. (A1). The distribution expressed in Eq. (5) (which is the normal Dalitz plot distribution) is unique in the sense that the invariant mass-squares t and u do not depend on the frame of reference.
This must be noted that when X is not a parity eigenstate, the differential decay rate is easily obtained, analytically, by simply adding the two contributions one would have if X had definite parity, i.e.
From Eq. (6) it is clear that, in general, The decay rate Γ (J/ψ → µ − µ + X s ) also gets further enhanced due to the fact that when X is not a parity eigenstate the allowed parameter space has a much wider spread as seen from Fig. 3 vis-à-vis Fig. 2. For example, the value of g 0+ could be about 5 times larger when X is not a parity eigenstate than when X is purely scalar, and thus one can expect that the number of J/ψ → µ − µ + X 0 events for the generic spin-0 case would be larger than 25 times the number of J/ψ → µ − µ + X 0 + events for the scalar case. Similar arguments also hold for the spin-1 case. Instead of the t vs. u Dalitz plot considered above, one could also think of another normal Dalitz plot using E − and E + , where E ± denotes the energy of µ ± in the rest frame of J/ψ. These two Dalitz plot distributions are related to one another by arising from the simple observation that in the rest frame of J/ψ we have For our analysis either of these normal Dalitz plots would serve the same purpose. However, for our numerical study in Sec. IV we have chosen the t vs. u Dalitz plot due to its Lorentz invariant nature. To arrive at another useful distribution we need to consider the Gottfried-Jackson frame of reference, which is essentially the center-of-momentum frame of the muon pair. Let the angle subtended by the direction of flight of µ − in this frame with respect to the direction of flight of J/ψ (equivalently that of X s ± as well) be θ, see Fig. 5. In this frame of reference we find that, with the Källén function being given by λ (x, y, z) = x 2 + y 2 + z 2 − 2 (xy + yz + zx) .
Using Eqs. (9) and (10) we can make a change of variables and write down the following expression for the differential decay rate in the Gottfried-Jackson frame in terms of s and cos θ, Note that while going from Eq. (5) to (11) the effect of time dilation has been taken into consideration, and the frame dependency of Eq. (11) essentially stems from terms proportional to various powers of cos θ after substituting t and u using Eq. (9).

C. The square Dalitz plot distribution
Another alternative distribution, favored by experimentalists, is the one called square Dalitz plot. The idea of square Dalitz plot first appeared in a BaBar paper [21], and recently it has been used in a LHCb study [22] as well. For our case, the expression for distribution of events in the square Dalitz plot can be obtained from Eq. (11) by making the following change of variables, such that both m and θ vary between 0 and 1. It is easy to show that where This square Dalitz plot provides a clearer view of the peaks in the differential decay rate as compared to the normal Dalitz plot and as we will see in Sec. IV they play a vital role in our study.

D. Information in the angular distribution
Before we study the patterns in the various distributions mentioned above, it is helpful to do a detailed analysis of the angular distribution, Eq. (11), in particular. We find that the angular distribution for all the spin-parity possibilities can be described in a unified manner by using the orthogonal Legendre polynomials in the following manner, P 2 (x) = 1 2 3x 2 − 1 , P 4 (x) = 1 8 35x 4 − 30x 2 + 3 denote the Legendre polynomials of order 2 and 4 respectively, and the angular coefficients T s± , U s± and V s± are as given in Appendix B.
Due to the orthogonal nature of the Legendre polynomials in Eq. (15), all the individual terms T s± , U s± and V s± explicitly shown in Eq. (B1) of Appendix B can be obtained via the following integrations following the method of moments [23,24], Extraction of these angular coefficients by Eq. (17) help us to define the following two ratios, which are independent of the coupling constants g s± . From Eqs. (B1i) and (B1j) in Appendix B it is clear that for spin-0 case (both scalar and pseudo-scalar) R 0± 2 = 0 while from Eqs. (B1k) and (B1l) we get R 1± 2 0 for spin-1 case. Therefore, R s± 2 can easily distinguish between spin-0 and spin-1 cases. Nevertheless, the ratio R s± 1 also has distinct features for the various spin-parity possibilities. We shall numerically show in the next section that these two ratios can indeed be used to distinguish the various spin-parity possibilities.
It should be noted that when X is not a parity eigenstate, the angular distribution would have the following form, where We can also define two ratios similar to the ones in Eq. (18), which, however, do depend on the coupling constants g s± . Interestingly when X is not a parity eigenstate, we only need to know its spin, which can be easily inferred from R s 2 alone, since for spin-0 case R 0 2 = 0 and for spin-1 case R 1 2 0. It should be noted that although the angular distribution of Eq. (11) is defined in the Gottfried-Jackson frame, the ratios as defined in Eqs. (18) and (21) are valid in all frames of reference, since they are function of the invariant masssquare s alone.
It is possible to do the integration over s and define the following ratios as well, Once again, for the spin-0 case we expect R 0± 2 = 0 and for the spin-1 case R 1± 2 0. There is one easy prescription to numerically measure the ratios R s± 1,2 . For this we can consider the t vs. u Dalitz plot, or the s vs. cos θ distribution of events. For every event we calculate the values of Y = a − m 2 µ 2 − b 2 cos 2 θ 2 , P 2 (cos θ) and P 4 (cos θ). Then we sum over all the values of Y, Y P 2 (cos θ) and Y P 4 (cos θ) to obtain Y , Y P 2 (cos θ) and Y P 4 (cos θ) respectively. Finally we use the following expressions to measure the values of R s± 1,2 , Within experimental accuracy, if R s± 2 = 0, then we have a spin-0 X boson, else the X boson has spin-1.

A. Canonical branching ratios and expected number of events at BESIII
Before we analyze the distribution patterns and ratios to distinguish the various spin-parity possibilities of X, it would be helpful to know whether such decays can indeed be probed experimentally. For this study we consider the BESIII experiment as a natural choice because at BESIII a large number (∼ 10 11 ) of on-shell J/ψ mesons will be produced at rest [25][26][27] which would help in inferring the 4momentum of the invisible X in a much easy and clean manner. Let us denote the canonical branching ratios for the decays J/ψ → µ − µ + X s ± by Br (J/ψ → µ − µ + X s ± ) and they are defined as the usual branching ratios sans the coupling constants i.e.
where Γ J/ψ is the experimentally measured total decay rate of J/ψ. In Fig. 6 we have shown the variation of the canonical branching ratios with the mass of X for its various spinparity possibilities. The distinct trend of rise in the canonical branching ratio for the axial-vector case in the limit m X → 0 is easily discernible from the m −2 X factor in Eq. (A1d). It must be noted that in evaluation of the canonical branching ratios shown in Fig. 6 we have included a cut on the missing energy at E missing = 140 MeV as suggested in Ref. [18] to completely eliminate the SM background events arising from final state radiation in the form of J/ψ → µ − µ + γ, especially the soft photon ones since at BESIII photons with energy less than 20 MeV can not be detected at the detector [28]. Details of the background study can be found in Ref. [18]. The canonical branching ratios are helpful in predicting the number of signal events, i.e. number of J/ψ → µ − µ + X s decays, one would expect at BESIII if the muon-philic X boson were to describe the muon g − 2 anomaly. In Fig. 7 we show the expected number of such signal events across the parameter space allowed by anomalous magnetic moment of muon (region shown before in Fig. 3), for three representative masses m X = 10 MeV, 110 MeV and 210 MeV. The high number of events shown in Fig. 7 are possible only when one considers the X boson to not be a parity eigenstate. This is due to the widening of the parameter space allowed by the presence of additional coupling constants g 0− and g 1− for the spin-0 and spin-1 cases respectively.
It is nevertheless interesting and important to numerically study and find out how to distinguish the individual spinparity possibilities irrespective of the muon anomalous magnetic moment consideration. For this we focus on a numerical study of the various distributions and the ratios that were presented, analytically, in the previous section.

B. Patterns in distributions, ratios and discussion
To study how the patterns of distribution differ for the different spin-parity possibilities of X as well as how these patterns evolve when we go from lower values of m X to higher values of m X , we once again consider the three representative masses, m X = 10 MeV, 110 MeV and 210 MeV. We will analyze the t vs. u Dalitz plot distribution of Eq. (5), the s vs. cos θ angular distribution of Eq. (11) and the m vs. θ distribution also known as the square Dalitz plot distribution as described by Eq. (13). To make a meaningful comparison across the various spin-parity possibilities, we have normalized the distributions by their maximum values which helps the patterns become more visually acute and easy to grasp. Due to this normalization by the maximum value of the distribution, all the distributions shown do not depend on the size of the coupling constants g s± .
To put the distribution patterns into proper perspective we have also included distribution of 1000 simulated Monte Carlo events generated after applying a cut on the missing energy at 140 MeV as mentioned before.
a. Normal Dalitz plot distribution: Let us first consider the t vs. u Dalitz plots, see Fig. 8. It is clear that for m X = 10 MeV, all the four distributions, corresponding to the four distinct spin-parity possibilities, are extremely flat in most of the allowed region except at the kinematic boundary. If we go to m X = 110 MeV, the distribution peaks are found to be indeed close to the boundary, but are more prominent than before. Here we can see that for spin-0 case there is no peak close to origin, while for spin-1 case the peak sits nearer to the origin. Considering m X = 210 MeV we find that the distribution still retains the characteristic locations of the peaks as mentioned before, but the scalar and pseudo-scalar distributions can now be apparently dis-tinguished, whereas its still difficult to distinguish between vector and axial-vector cases.
b. The s vs. cos θ distribution: Considering the s vs. cos θ distributions as shown in Fig. 9, we find that for m X = 10 MeV, there are not enough useful distinguishing features. However, if we consider m X = 110 MeV and 210 MeV, it is possible to distinguish between all the four spin-parity cases, albeit the fact that the distribution peaks are still close to the boundary of the plots.
c. Square Dalitz plot distribution: Considering the square Dalitz plot distributions as shown in Fig. 10 for m X = 10 MeV, 110 MeV and 210 MeV, we find that the differences between all the four spin-parity possibilities which were kind of hidden in other distributions, do show up very clearly. In fact for the low mass of m X = 10 MeV the differences are much clearly visible, when we compare, for example, the patterns of vector and axial-vector cases of m X = 210 MeV. The square Dalitz plot is indeed a very powerful discerning distribution due to the fact that the distribution peaks which are located at the boundaries in other distributions correspond to internally observed peaks here. d. Ratios: Finally looking at the ratios R s± 1,2 , we find from Figs. 11 that using R s± 2 = 0 as a criteria one can easily distinguish the spin-0 and spin-1 cases. For m X = 10 MeV case we find that for axial-vector scenario the R 1− 2 is extremely tiny but non-zero. In this case one can clearly use the curves for R 1± 1 to distinguish the various spin-parity possibilities without any ambiguity. Therefore, both R s± 1 and R s± 2 should be used to distinguish among the four spin-parity possibilities.
It would indeed be more practical to employ both the ratios and the distributions, especially the square Dalitz plot distribution to decipher the spin-parity of the X boson in the decay J/ψ → µ − µ + X.

V. CONCLUSION
The present anomaly in muon anomalous magnetic moment can be considered as a tantalizing hint of the presence of some new physics. The simplest new physics that can contribute involves purely muon-philic interactions and one needs only to hypothesize that there exists a spin-0 or spin-1 boson (say, X) taking part in these interactions. If the X boson is a parity eigenstate, then only scalar and vector possibilities are allowed from the muon anomalous magnetic moment. Interestingly the allowed paramater space for the coupling constants widens if one considers X to be not a parity eigenstate.
Irrespective of the parity of the X boson, and considering its mass m X 2m µ as suggested by the muon anomalous magnetic moment and other experimental studies, the electrically neutral X boson can be definitively probed by the decay J/ψ → µ − µ + X. This decay not only provides an excellent avenue to test this specific muon-philic new physics possibility, but it also can be utilized to ascertain whether X is a parity eigenstate or not, and if it is a parity eigenstate we can probe both its spin and parity. Due to the widening of parameter space for parity non-eigenstate scenario, decay with such a X would have larger branching fraction than the one which is a parity eigenstate.
Since the decay J/ψ → µ − µ + X is a three-body decay, there are a multitude of distributions that can be studied. We find that the square Dalitz plot distribution as well as the angular distribution of events in the center-of-momentum frame of the muon pair, do indeed exhibit the distinguishing features among all the four spin-parity possibilities, viz. scalar, pseudo-scalar, vector and axial-vector cases. We also provide two dimensionless ratios, obtained from the angular distribution, which are also useful to ascertain the spinparity of X. The search for decay J/ψ → µ − µ + X or similar decay such as φ → µ − µ + X would constitute an important strategy in our search for definitive signature of new physics.