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.


Introduction
Signals with missing energy at hadron colliders commonly arise in many new physics models solving the dark matter problem of the Universe. Even in the Standard Model (SM), the neutrinos are undetectable, so recorded as missing energy. Though being ubiquitous, the missing energy has been a challenging object that prevents the full reconstruction of particles involved in the decay process. As will be discussed in Sec. 2, it is because of the fundamental absence of inverse projection from the phase space of visible particle momenta to the full phase space. Still, there are a plethora of useful methods and algorithms proposed for determining the mass spectrum of the underlying dynamics from missing energy events [1,2].
A mathematical insight on the missing energy kinematics led to the invention of the algebraic kinematic method [3]. It was realized that the phase space of visible particles might possess identifiable singularities, from which one could extract the mass spectrum information for given event topology. Then, an optimized one-dimension variable called the singularity coordinate has been proposed to capture the singular behavior in an effective way. The distribution of the singularity coordinate becomes singular, i.e., having a sharp peak or distinct edge, when the input masses equal to the true values. Though powerful and insightful, it has not been widely considered as applicable for practical new physics searches, partly due to the lack of concrete prescriptions with more examples. Another obstructing factor is that the singularity coordinate is an implicit function of the mass spectrum involved in the decay process. In the absence of good ansatz, it is necessary to perform multi-dimensional fitting, which is impractical at the stage of discovery. Furthermore, as the value of the singularity coordinate does not relate directly to physical parameters, it is not trivial to interpret the singularity coordinate in terms of the physical quantities. The situation has recently improved due to the studies in Ref. [4], where the singularity method is reexamined and expanded for various event topologies.
Motivated by the results in Ref. [4], we examine the singularity method for the antler decay topology [5,6,7], and propose a derived singularity variable. Many new physics resonances decaying into the final state with missing energy can be represented by the antler decay topology. 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, it can be useful for measuring the resonance mass accurately. A brief overview of the singularity method and the derivation of the new singularity variable are presented in the next section.
In Sec. 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 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 ) ∈ R n | g i (p 1 , . . . , p n ) = 0 for all 1 ≤ i ≤ m} . (1) However, since not all final-state momenta are measurable at collider experiments, we cannot always 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 humanmade 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. 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 frame hard, or even impossible, on an event-byevent 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 (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, n j=1 ∂g i ∂p j (p j − q j ) (i = 1, . . . , m).
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 → ν [8] and h → W W → 2 + / E T processes [9] 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 X(Q) and a pair of invisible particles via unstable intermediate states [5,6,7], as shown in Fig. 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 two-Higgs doublet model (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 directly be inferred, unlike simple invariant or transverse masses. 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 Fig. 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. For simplicity, 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, Using the condition of energy-momentum conservation Q = p 1 +p 2 +k 1 +k 2 , the kinematic constraints for the antler decay topology (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 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 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 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. 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 Fig. 2, we show the ∆ AT distribution for pure phase space, without assuming the underlying model. 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 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 , 3 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 2 Here we have performed some row and column operations to simplify the expression. 3 In Ref. [4], MX and MB are unknown, while MA is to be determined. Then, M antler as the solutions to the singularity variable has been introduced. We think that our consideration is more common in new physics searches at hadron colliders. Furthermore, as can be seen in Sec. 3, we perform a more realistic study beyond the phase space as well. variables, we provide the coded implementations in Haskell [10] and C ++ [11]. 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 Fig. 2. Note that the two distributions overlap only at the true M X value. It means that if a unique solution for Eq. (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 distributions when M AT = M true X leads us to conclude that 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 Sec. 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,9], 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 The corresponding ∆ AT is still essentially a quadratic polynomial. In Fig. 3, we have shown the ∆ 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 [12]. The formula was used for illustration purposes only. We will properly use the parton distribution function in the next section. added to each panel. By ignoring the longitudinal momentum, the relation in (11) is no longer valid. We, however, find that the distributions still have the endpoints near the resonance mass, although they are not as sharp as those of the true distributions and slightly smeared. Thus, we expect that the M (0) 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 M AT variable will further enable us to measure the resonance mass accurately. We will employ one example of such an estimator 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 [13,14] favor the alignment limit, where one of the neutral Higgs bosons has the SM-like couplings to the SM vector bosons. In the alignment limit, the non-SM-like Higgs bosons H, A, and H + interact among themselves more strongly, while the branching ratios of H → W W , ZZ, and hh are all suppressed. Moreover, if the heavy Higgs bosons are mass-degenerate as in the decoupling limit of the supersymmetric model, the dominant decay mode is H → tt unless 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. [15] 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 fully-hadronic [16] and semi-leptonic decay modes [17], but the CMS collaboration has recently performed the search using the di-leptonic final state [18]. In the CMS analysis, the invisible neutrino momenta have been obtained by directly solving the kinematic constraints, following the method in Ref. [19]. We note that though the singularity variable is also based on the kinematic constraints, it is computationally cheaper as well as powerful.
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. 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 [20], interfacing with LHAPDF 6 [21] for parton distribution functions. We use the NNPDF parton distributions [22], 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.
The parton-level distributions for the M there is a combinatorial ambiguity of pairing the visible particles into two sets. In Fig. 4, we have used the correct pairing. We will address this issue shortly.
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 [23,24,25,26]. 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. 5 The M T 2 variable can be used in the presence of two invisible particles in the final state [28,29]. 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 the singularity condition becomes the quartic polynomial of Q 0 = (M 2 AT + Q T 2 +(Q maos Therefore, the total number of M 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. We show the M AT distributions using the MAOS method in Fig. 5. One can see that the endpoints of the M AT distributions are more pronounced than M (0) 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 [23], 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. [30] and the improved version [31] 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 [32]. In our simulation, the anti-k T clustering algorithm [33] 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 [34]. We set the cone sizes for isolating the electrons and muons to be 0.4. The jet energy scale correction is fixed 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 [18]: • 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.
• 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 Fig. 6, we show the M AT distributions for the detector-level event data. Here we have used the algorithm in [30] 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 . By checking the distributions at jet level without the detector effects, we have found that the peak shift is mainly due to the jet reconstruction. The momentum smearing effect and basic selection cuts do not distort the overall shape of the distributions. For M max AT , the peak is located at slightly above the m H value by the same reason. It is worth scrutinizing the variables in more detail at jet level, but it is beyond the scope of our study.
The simulation results show that one can deduce the heavy Higgs mass by investigating the edge of the M min AT and the peak of M max AT distributions. The signal distributions remain to be well separated from the background, even when taking the detector effects into account. 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. We have identified the singularity condition using the prescription in Ref. [4] and derived a onedimensional 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.