A singular way to search for heavy resonances in missing energy events

The phase space of visible particles in missing energy events may have singularity structures. The singularity variables are devised to capture the singularities effectively for given event topology. They can greatly improve the discovery potential of new physics signals as well as to extract the mass spectrum information at hadron colliders. Focusing on the antler decay topology of resonance, we derive a novel singularity variable whose distribution has endpoints directly correlated with the resonance mass. As a practical application, we examine the applicability of the singularity variable to the searches for heavy neutral Higgs bosons in the two-Higgs doublet model.

The novel singularity variable has a direct correlation with the mass scale of the resonance, so it can greatly help distinguish the new physics signal from backgrounds. And, along with the inclusion of supplementary approximation for the unknown longitudinal momentum of the resonance, the endpoint shape of the singularity variable distributions becomes more pronounced, so we expect that it can be useful for the mass measurement as well. A brief overview of the singularity method and the derivation of the new singularity variable are presented in the next section.
In section 3, we evaluate the performance of the singularity variable in the case of a heavy Higgs boson decaying to a top pair. We compare the signal distributions to the dominant SM backgrounds and take the detector effects into account, to check the viability of the singularity variable in more realistic collider searches. Since the unknown longitudinal momentum of the Higgs boson produced at hadron colliders is not negligible, we improve the singularity variable by employing an approximation scheme.

Kinematic singularity of the antler decay topology
The phase space is the hypersurface of the final-state particle momenta, subject to the kinematic constraints like energy-momentum conservation and on-shell mass relations. In other words, the particle momenta reside in the solution space of coupled polynomial equations of total degree two. The set of all solutions of a system of polynomial equations is called an affine variety in mathematics: Π(g 1 , . . . , g m ) = {(P 1 , . . . , P n ) ∈ E n | g i (P 1 , . . . , P n ) = 0 for all 1 ≤ i ≤ m} , (2.1) where P i are the four-momenta of final state particles, and g i are polynomials for kinematic constraints such as on-shell mass relations. E is the four-dimensional pseudo-Euclidean space, R 1, 3 . We often confront the situation where we cannot fully reconstruct the phase space. It is because the final states of collider events may contain invisible particles such as neutrinos or dark matter candidates, which escape human-made detectors, as well as visible particles. The missing energy events are one of the typical types of signals predicted by many new physics models beyond the SM that we eager to discover at the LHC and future colliders. We decompose the phase space into {(p i , k j )}, where p i are the four-momenta of visible particles, and k j are those of invisible particles. In terms of the phase space, what detectors at collider experiments are doing is the projection of the full phase space {(p i , k j )} onto the space of visible momenta {(p i )}, up to finite detector resolution and acceptance, for the missing energy events. The projection makes the reconstruction of the center-of-mass JHEP07(2020)089 frame hard, or even impossible, on an event-by-event basis because it is not invertible. The only invertible projection is the identity mapping, which is unachievable by construction.
One important property of the affine variety is the presence of singular points or singularities. In ref. [3], it was noted that the projected visible phase space could possess singularities, where the tangent plane fails to exist, even though the full phase space including the invisible momenta was regular at the points. To formulate it, we introduce the Jacobian matrix J, the m × n matrix of partial derivatives: where g i are the defining polynomial of the phase space in (2.1). If the point Q is regular, the Jacobian matrix provides the linear approximation of the phase space near the point. It is the tangent space at Q, 3) The rank of J equals the dimension of the phase space. At the singularity, the Jacobian matrix has a reduced rank, lower than the rank at the regular point. In the case where m = n, the determinant of the Jacobian matrix vanishes at the singularity. In ref. [3], it is considered the restricted matrix composed of the derivatives for the invisible momenta, rather than the full Jacobian matrix. Then, a singularity coordinate has been proposed to exploit the singularity structure of the visible phase space. 1 It is an implicit function of the mass spectrum involved in the hypothesized decay process for given events. Once the hypothesis was correct, the singularity coordinates maximize the singular features at the true mass values, thus enabling us to determine the mass spectrum. Though being a simple and elegant idea based on mathematical constructions, one has failed to find practical examples except for the applications to the W → ν [12] and h → W W → 2 + / E T processes [13] at hadron colliders. It is due to the lack of more concrete examples and programmed implementations for practitioners. Recently, a set of worked-out examples for various event topologies have been thoroughly investigated and visualized in ref. [4], which deepens the understanding of the singularity variables and makes them more approachable.
Inspired by the studies in [4], we here concentrate our attention on the antler decay topology and propose a novel kinematic variable derived from the singularity variable. In the antler decay diagram, a singly produced heavy resonance decays into a pair of visible particles and a pair of invisible particles via unstable intermediate states [5][6][7], as shown in figure 1. It represents a typical decay pattern of heavy resonances in many new physics models such as the decays of the heavy Higgs bosons in the 2HDM, as well as the h → W W → 2 + / E T process in the SM. A technical subtlety of the singularity coordinate is that it is an implicit variable whose relation with the mass spectrum cannot JHEP07(2020)089 directly be inferred. It differs from invariant and transverse mass variables, whose peak or edge directly provides the mass information. On the other hand, the peak position of the singularity coordinate is not related to the mass value. What the singularity coordinate promises is that the distribution exhibits singular behavior when the hypothesis on the mass spectrum is correct. One can attempt the maximum likelihood estimation for the mass spectrum using template distributions, but it is not suitable at the stage of discovery. It is clear that kinematic variables would be more useful if they enable us to directly deduce the mass scale of heavy resonances by identifying the positions of the peak or edge of their distributions. To derive such a kinematic variable, we review the singularity variable for the antler decay topology at first.
As shown in figure 1, we consider a heavy resonance X, decaying into two visible v 1, 2 and two invisible particles B 1, 2 through unstable intermediate states A 1, 2 , Notice that the diagram also covers the case of a pair of cascade decay chains. One can define v i to be the collection of visible particles at the chain i. In this case, the invariant mass of the visible particle system m i is not constant, but an event variable. We take into account the case where all the particles are on mass-shell and the decay chains to be symmetric: B i corresponds to the neutrino or a dark matter candidate, yielding the missing transverse energy, In the cases where the intermediate states are off mass-shell, or their decay widths are large, the on-shell mass relations are not valid. It results in a decreased number of kinematic constraints. While the singularity method is still valid, we may not use the condition of vanishing Jacobian matrix, but have to find the reduced rank condition of the Jacobian.

JHEP07(2020)089
Using the condition of energy-momentum conservation Q = p 1 + p 2 + k 1 + k 2 , the kinematic constraints for the antler decay topology (2.4) are given as: And, the elements of the Jacobian matrix are obtained by taking partial derivatives on the constraint equations, for k 1j = (k 10 , k 1T , k 1L ). Since it is a square matrix, the singularity condition implies that its determinant is zero. To eliminate the invisible momentum components, we adopt a simple trick used in ref. [4], where det JηJ T with η = diag(1, −1, −1, −1) is taken instead of det J. Then, we find that the singularity condition is equivalent to the vanishing of 2 Here the subscript "AT" stands for the antler decay topology. Note that the transverse momentum of the heavy resonance can be determined by measured quantities because Thus, ∆ AT is a function of Q 0 = (M 2 X + Q T 2 +Q 2 L ) 1/2 and Q L if M A and M B were known a priori or determined by ansatz. Given Q L , ∆ AT is a quartic polynomial equation of Q 0 with nonzero coefficients, in general. If the resonance was produced at rest, i.e., Q 0 = M X and Q L = 0, ∆ AT is a quartic equation of M X in the form like M 2 X (aM 2 X + bM X + c). The trivial solution M X = 0 is the side effect of the determinant trick, which is unphysical. Therefore, ∆ AT reduces to a quadratic equation of M X in essence.
In order for a numerical study, we set M X = 800 GeV and M A = 173 GeV, while B i is massless. The invariant masses of visible particles m i are not vanishing, but are set to be varying event-by-event between 0 and 153 GeV. 3 These numbers have been chosen to match those used in the study of the next section. We also assume that the decay widths of the resonance and intermediate particles are negligible. In the left panel of figure 2, we show the ∆ AT distribution for pure phase space, ignoring the spins and parities of all the particles involved in the decay process. In order to see the impact of the singular feature on the masses, we have used the true spatial momentum of the resonance, while the mass M X is taken to be an input parameter. In other words, we consider the case where the JHEP07(2020)089 resonance was produced at rest. For M X = M true X , we can see a remarkably sharp peak at ∆ AT = 0, which confirms that the singular feature is maximized at the true mass values. Note that ∆ AT is negative since det JJ T > 0, while det η < 0 for the correct choice of the mass parameters. Meanwhile, the distributions for M X = M true X do not exhibit such a sharp peak: the event number densities around ∆ AT = 0 are not particularly noteworthy.
The observation suggests that the solutions of the singularity variable ∆ AT have a direct and strong correlation with the resonance mass. We call the solutions M AT , 4 (2.10) Although there exist general analytic solutions to the polynomial equations up to fourth order, it does not give us further insight, nor particularly useful. For practical use of the ∆ AT and M AT variables, we provide the coded implementations in Haskell [14] and C ++ [15]. As mentioned earlier, since ∆ AT is essentially a quadratic polynomial equation if X is at rest, there are up to two degenerate solutions to the equation. We sort the solutions into the smaller and the larger, and label them as M min AT and M max AT , respectively. The M AT distributions are shown in the right panel of figure 2. Note that the two distributions overlap only at the true M X value. It means that if a unique solution for eq. (2.10) exists for given events, it corresponds precisely to the resonance mass. One can see that the singularity of ∆ AT has transformed into the peak of the M AT distribution at M true X , and both M min AT and M max AT contribute to the peak. Moreover, the non-overlapping of the two JHEP07(2020)089 This relation implies that we can directly deduce the mass scale of M X from the edge of M min AT and the threshold of M max AT distributions. In practical applications, we expect that both variables can serve as important cuts to enhance the signal-to-background ratio, in particular at the stage of discovery. We will see this aspect more clearly in section 3 with a more concrete example in the presence of background.
Up to now, we have assumed that the longitudinal momentum of the resonance is known a priori. However, it is practically unattainable for hadron collider events. In practice, the resonance is boosted in the longitudinal direction by the net momentum sum of colliding partons, which follows parton distribution functions. The amount of the longitudinal boost is unknown on an event-by-event basis. It is an intrinsic nature of hadron colliders, one of the biggest obstacles in particle object reconstructions. In refs. [4,13], the longitudinal momentum of resonance has been neglected by assuming that the colliding gluons have nearly the same amount of energy. In our numerical simulation with parton distribution functions, we found that the longitudinal momentum could be so large that it could not be neglected, even in the case of gluon-gluon collision. However, in the absence of a good estimator or approximator for the unknown longitudinal momentum, we have no other choice but to take an ad hoc solution. We ignore it by setting it to be zero,à la transverse mass variables. We denote the ad hoc solution as M The corresponding ∆ AT is still essentially a quadratic polynomial. In figure 3, we have shown the ∆ AT variable can still give us a hint of the mass scale of the resonance at the stage of bump hunting. If a good estimator for the longitudinal momentum is supplemented, the relation (2.11) can be approximately restored so that the endpoint and peak shapes of the M AT distributions become more pronounced. We will employ one example of such an estimator in the next section. 5 To simulate the distribution of the longitudinal momentum, we used a simple inverse power formula, f (QL) ∝ QL (1 + QL/Q * ) −n with Q * = 100 GeV and n = 5. The formula is inspired by the one used in the measurements of transverse momentum distributions of strange mesons [16]. The formula was used for illustration purposes only. We will properly use the parton distribution function in the next section.

Searching for heavy Higgs bosons decaying into a top pair
As an application of the M AT variable, we consider heavy neutral Higgs bosons in the 2HDM. The current results on the SM Higgs measurements [17,18]  the ratio of Higgs vacuum expectation values tan β is very large in the type II model. For a review on the heavy Higgs boson decays in the alignment limit of the 2HDM, see ref. [19] and references therein.
Though the top-pair process is an important channel to search for the heavy Higgs boson, it must overcome the SM tt background, which has a huge cross-section and possesses exactly the same final state as the Higgs signal. Furthermore, the di-leptonic final state 2b+ 2 + / E T contains two neutrinos that prevent from reconstructing the center-of-mass frame of the Higgs boson. The ATLAS collaboration at the LHC is still focusing on the fullyhadronic [20] and semi-leptonic decay modes [21], but the CMS collaboration has recently performed the search using the di-leptonic final state [22]. In the CMS analysis, the invisible neutrino momenta have been obtained by following the method in ref. [23]. They directly solve the series of equations for the kinematic constraints and assign a weight to each solution to characterize how likely it is to occur in tt production based on the expected true b-lepton invariant mass spectrum. Then, kinematic quantities, such as the invariant mass of the tt system, are calculated as a weighted average. We note that though the singularity variable is based on the same kinematic constraints, we do not attempt to obtain the invisible neutrino momenta. Instead of solving all the kinematic constraints, we introduce the singularity condition as given in (2.10) and solve only the one polynomial equation. Therefore, it is computationally cheaper than the kinematic reconstruction method used in the CMS analysis.
For the numerical study, we have chosen a benchmark point: the CP-conserving type II model with M H = M A = M H + = 800 GeV, tan β = 3, and the Higgs mixing angle cos(β − α) = 0.01. The benchmark point is safe from the existing LHC bounds set by the SM Higgs decay width and coupling measurements and direct searches for the heavy Higgs bosons [24]. In our estimation, the branching ratio of the H → tt process is about 93%, while the subleading process is H → bb with the branching ratio of 3%. The total decay width of the heavy Higgs boson is 2.7 GeV. We have generated parton-level event samples for both signal and background using Pythia 8 [25], interfacing with LHAPDF 6 [26] for parton distribution functions. We use the NNPDF parton distributions [27], and set the proton-proton collision energy to be 13 TeV. In our estimation, we find that the gluon-fusion process dominates the Higgs production, but the bottom-fusion process is non-negligible as well. The latter contributes to about 15% of the total cross section at leading order. It is because of the enhanced Higgs coupling to the bottom quarks by tan β. In the case of the type I model, the bottom-fusion process are suppressed, so the gluon-fusion contribution will be predominant. As mentioned in the previous section, a good estimator for the longitudinal momentum of the heavy Higgs boson can improve the singularity variables. It will be particularly in need when attempting the mass measurement after discovery. One possible option worth to examine is the M T 2 -assisted on-shell (MAOS) method [28][29][30][31]. The MAOS method provides an approximation for the longitudinal momenta of invisible particles by solving the on-shell relations, given the solution of the M T 2 variable for the transverse momenta. 6 The M T 2 variable can be used in the presence of two invisible particles in the final state [33,34]. The endpoint of the M T 2 distribution is m t for both signal and background because the top quarks were produced on mass-shell. From the two quadratic on-shell relations, the MAOS method yields up to four possible solutions for the unknown longitudinal momenta. By setting the longitudinal momentum of the resonance or the tt system to be From the solution of the polynomial equation in the above, we have Note that the total number of M maos AT for given event is now up to 16. As we find no plausible criterion to choose a particular solution, we shall use all the real solutions for M AT , discarding the complex ones. As the number of solutions has been increased, we define AT , and lie at the heavy Higgs boson mass. It is known that the accuracy of the longitudinal momentum approximation can be improved by imposing an M T 2 cut [28], so it can make the endpoint shape more distinct. Therefore, we expect that the M AT with the MAOS method can serve as a useful variable to measure the resonance mass accurately with the M T 2 cut.
We here briefly leave a comment on the combinatorial ambiguity on the pairing of visible particles. There are two possible ways to pair one b quark and one charged lepton in each event. One can resolve the ambiguity by using various kinematic variables. In particular, the algorithm proposed in ref. [35] and the improved version [36] provide a good efficiency 80% for tt events. In practice, the combinatorial ambiguity does not interfere with the endpoints of M min AT and M max AT distributions. One computes M AT for all possible pairing and then take the minimum or maximum. The positions of the endpoints are intact by taking the extrema. Nonetheless, we have checked that the algorithms work well for both signal and background events, and the combinatorial ambiguity does not affect the overall shape of the M AT distributions when using the algorithms.
Finally, we examine the M AT variable including jet reconstruction, object isolation, event selection cuts, and various detector effects. The parton-level event samples have been processed by Pythia for parton shower and hadronization. Then, we employ FastJet 3 for reconstructing jets [37]. In our simulation, the anti-k T clustering algorithm [38] with a distance parameter R = 0.4 is chosen for the jet reconstruction. The object reconstructions and detector effects, such as flavor tagging, fake rates, and momentum smearing, have been performed by the fast detector simulation program DELPHES 3 [39]. We set the cone sizes for isolating the electrons and muons to be 0. to unity. Except for the settings described in the above, we use the default CMS card provided by DELPHES.
We have applied the basic selection cuts similar to those used in the CMS analysis to the detector-level events [22]: • Isolated electrons and muons satisfy p T > 20 GeV and |η| < 2.4. Only events with two oppositely charged leptons are selected. The leading lepton must have p T > 25 GeV.
• The invariant mass of isolated leptons is required to be larger than 20 GeV to suppress low-mass resonance events. And, in order to veto events with the Z boson, we reject events with 76 < m < 106 GeV.

JHEP07(2020)089
• We require that at least two jets with p T > 30 GeV and |η| < 2.4, and at least one of the jets are b tagged.
• The missing energy is required to be larger than 40 GeV.
In figure 6, we show the M AT distributions for the detector-level event data. Here we have used the algorithm in [35] to resolve the combinatorial ambiguity. In the M min AT distributions, the edge structure is clearly visible around the m H value, and the background distributions are populated near the threshold. Meanwhile, the peak has been shifted towards below m H . For M max AT , the peak is located at slightly above the m H value. One of the possible reasons for the peak shift is that some of the decay products of the b quark could be missed since b jets tend to be wider and have higher multiplicities than light jets. Consequently, the b jet may not catch the momentum of the initiating b quark well enough. It can be improved by investigating flavor-dependent jet energy scale corrections and comparing various jet reconstruction algorithms and the inherent parameters associated with the algorithms. On the other hand, we find that the momentum smearing effect and basic selection cuts do not distort much the overall shape of the distributions. It is worth scrutinizing the variables in more detail at jet level, but it is beyond the scope of our study. We note that the M AT variables are not affected by initial state radiations since the variables are derived from the Lorentz-invariant kinematic constraints, given in (2.6), in which the radiations do not participate.
The simulation results show that one can deduce the heavy Higgs mass by investigating the M min AT and M max AT distributions. The M min AT distribution has the edge structure at the heavy Higgs mass, even when taking into account the jet reconstruction and detector effects. The MAOS approximation method for the longitudinal momentum of the heavy Higgs boson makes the edge structure more clearer. One can see that the MAOS method has improved the shapes of the M AT distributions by comparing the upper and the lower panels in figure 6. The threshold of the signal M min AT distribution can be buried by the background distribution, but the peak of the distribution may give us the information of the upper bound on the heavy Higgs mass. Furthermore, the signal distributions remain to be well separated from the background. This observation suggests that the M AT variable can also be useful for setting the signal or control regions when searching for heavy resonances.

Conclusions
The algebraic singularity method was proposed from the observation that the projected visible phase space can have singularities in the presence of missing energy. The singularity variables have been devised to capture such singular features, and they implicitly provide the mass spectrum information of intermediate resonances and invisible particles in the final state. Recently, it has been outlined the prescriptions for deriving the singularity variables for various event topology with missing energy [4]. It makes the singularity method more accessible for practical applications.
In this article, we focused on the antler decay topology, where a heavy resonance decays into the final state of two visible and two invisible particles through intermediate states.

JHEP07(2020)089
We have identified the singularity condition using the prescription in ref. [4] and derived a one-dimensional variable that has a strong correlation with the resonance mass. The M AT variable is the solution to the polynomial equation of the singularity condition. We have confirmed that the phase-space distributions of the minimum and maximum of the M AT have endpoints at the correct resonance mass value, thus enabling us to measure the mass as well as to discover the resonant signal.
As a practical application, we have studied the signal of the heavy Higgs bosons decaying into a top pair, one of the typical signals in the 2HDM. Although the ignorance of the longitudinal momentum of the heavy Higgs boson smears the endpoint structure of the M AT distributions, we find that the signal distributions are well separated from the background, and they can give us a hint for the mass scale of the heavy Higgs boson. Moreover, by employing an approximation scheme for the longitudinal momentum, the endpoint structure can become sharper. We have used the MAOS method to exemplify our proposal and found that it successfully restores the shape of the M AT distributions. The feature remains intact even in the presence of the detector effects and selection cuts, on the whole. From these observations, we expect that the M AT can serve as the main variable of cut-based analyses, or an important input feature of multivariate studies, for heavy resonance searches at hadron colliders.