Using the (Modified) Matrix Element Method to constrain $L_\mu - L_\tau$ Interactions

TL;DR

This paper demonstrates that the matrix element method significantly enhances the sensitivity to a $Z'$ gauge boson in $L_-$ interactions at the LHC, especially in $Z o 4$ and $2+ ot{E}_T$ channels, compared to traditional methods.

Contribution

The study applies the matrix element method to constrain $L_-L_$ interactions, showing improved sensitivity over cut-and-count techniques for low-mass $Z'$ bosons.

Findings

Sensitivity to $g_{Z'} \u2265 0.002 g_1$ for $M_{Z'} < 20$ GeV

Sensitivity to $g_{Z'} \u2265 0.005 g_1$ for $20 < M_{Z'} < 40$ GeV

Order of magnitude improvement over traditional methods

Abstract

In this paper, we explore the discriminatory power of the matrix element method (MEM) in constraining the $L_\mu-L_\tau$ model at the LHC. The $Z'$ gauge boson associated with the spontaneously broken $U(1)_{L_\mu-L_\tau}$ symmetry only interacts with the second and third generation of leptons at tree level, and is thus difficult to produce at the LHC. We argue that the best channels for discovering this $Z'$ are in $Z \to 4\mu$ and $2\mu+\displaystyle{\not}E_T$. Both these channels have a large number of kinematic observables, which strongly motivates the usage of a multivariate technique. The MEM is a multivariate analysis that uses the squared matrix element to quantify the likelihood of the testing hypotheses. We find that with $300 \, \text{fb}^{-1}$ of integrated luminosity, we are sensitive to the couplings of $ g_{Z'} \gtrsim 0.002 ~ g_1$ and $M_{Z'} < 20 \text{ GeV}$, and $g_{Z'} \gtrsim 0.005 g_1$ and $20\, \text{GeV} <M_{Z'} < 40 ~ \text{GeV}$, which is about an order of magnitude improvement over the cut-and-count method for the same amount of data.

Tables

Table 1. Table 1 : Determining α a = ( M Z ′ a , ϵ a ) superscript 𝛼 𝑎 superscript subscript 𝑀 superscript 𝑍 ′ 𝑎 superscript italic-ϵ 𝑎 \alpha^{a}=(M_{Z^{\prime}}^{a},\epsilon^{a}) using a combination of likelihood maximization and the ℒ − 1 superscript ℒ 1 \mathcal{L}^{-1} distribution.

$M_{Z^{\prime }}^a=M_{Z^{\prime }}^{\mathrm{gen}}$	${ϵ}^a=0.05$
using likelihood maximization
$M_{Z^{\prime }}^a$ arbitrary	${ϵ}^a>0.05$
likelihood maximization fails, because of small ${ϵ}^{\mathrm{gen}}$

Equations

$$P(p^{\mathrm{par}}|\alpha )=\frac{1}{\sigma }\int d{x}_{{}_1}d{x}_{{}_2}\frac{f(x_{{}_1})f(x_{{}_2})}{2s{x}_{{}_1}{x}_{{}_2}}|\mathcal{M}(p^{\text{par}};\alpha ){|}^2{\delta }^4(p_{\text{initial}}^{\text{par}}-p_{\text{final}}^{\text{par}}),$$

$${\delta }^4(p_{\text{initial}}^{\text{par}}-p_{\text{final}}^{\text{par}})=$$

$${\delta }^2(p_{\text{final}}^{\text{par, transverse}})\times$$

$$\delta (\frac{\sqrt{s}}{2}(x_{{}_1}-x_{{}_2})-p_{\text{final}}^{\text{par, longitudinal}}).$$

$$P(p^{\mathrm{rec}}|\alpha )=\frac{1}{\sigma }\int d\Phi (p_{\text{final}}^{\text{par}})d{x}_{{}_1}d{x}_{{}_2}\frac{f(x_{{}_1})f(x_{{}_2})}{2s{x}_{{}_1}{x}_{{}_2}}|\mathcal{M}(p^{\text{par}};\alpha ){|}^2{\delta }^4(p_{\text{initial}}^{\text{par}}-p_{\text{final}}^{\text{par}})W(p^{\mathrm{rec}},p^{\mathrm{par}}),$$

$$\mathcal{L}(p_i^{\mathrm{rec}};\alpha )=\frac{P(p_i^{\mathrm{rec}}|\text{new phyiscs }(\alpha ))}{P(p_i^{\mathrm{rec}}|\text{null hypothesis})},$$

$$\sum _i\log [\mathcal{L}(p_i^{\mathrm{rec}};\alpha )]=\sum _i\log \left[\frac{P(p_i^{\mathrm{rec}}|\text{new phyiscs }(\alpha ))}{P(p_i^{\mathrm{rec}}|\text{null hypothesis})}\right].$$

$$\mathcal{L}\ni -\frac{1}{4}(Z^{\prime }{)}_{\alpha \beta }(Z^{\prime }{)}^{\alpha \beta }+\frac{1}{2}{M}_{Z^{\prime }}^2{Z}^{{}^{\prime }\alpha }{Z}_{\alpha }^{\prime }-ϵg_1{Z}_{\alpha }^{\prime }\left({\bar{\ell }}_2{\gamma }^{\alpha }{\ell }_2+\bar{\mu }{\gamma }^{\alpha }\mu -{\bar{\ell }}_3{\gamma }^{\alpha }{\ell }_3-\bar{\tau }{\gamma }^{\alpha }\tau \right),$$

$$\mathcal{L}(p^{\text{rec}};M_{Z^{\prime }},ϵ)=\frac{\frac{1}{\sigma }\frac{f(x_{{}_1})f(x_{{}_2})}{2s{x}_{{}_1}{x}_{{}_2}}|\mathcal{M}(p^{\mathrm{rec}};M_{Z^{\prime }},ϵ{)}_{q(x_{{}_1})q(x_{{}_2})\to 4\mu )}{|}^2}{\frac{1}{\sigma }\frac{f(x_{{}_1})f(x_{{}_2})}{2s{x}_{{}_1}{x}_{{}_2}}|\mathcal{M}(p^{\mathrm{rec}};SM{)}_{q(x_{{}_1})q(x_{{}_2})\to 4\mu )}{|}^2}.$$

$$\mathcal{L}(p^{\mathrm{rec}};M_{Z^{\prime }},ϵ)=\frac{|\mathcal{M}(p^{\mathrm{rec}};M_{Z^{\prime }},ϵ){|}^2}{|\mathcal{M}(p^{\mathrm{rec}};\mathrm{SM}){|}^2}$$

$$\mathcal{M}(p^{\text{rec}};M_{Z^{\prime }}^a,{ϵ}^a)={\mathcal{M}}_{\mathrm{SM}}+\frac{({ϵ}^a{)}^2}{M_{\mu \mu }^2-(M_{Z^{\prime }}^a{)}^2-i{\Gamma }_{Z^{\prime }}^a{M}_{Z^{\prime }}^a}(\cdots ),$$

$$\frac{\text{number of signal events}}{\sqrt{\text{number of background events}}}.$$

$$S$$

$$B$$

$$pp\to$$

$$(p_T^{\tau }{)}^{\mathrm{approx}}=\sqrt{\frac{M_Z^2}{2(1+\cosh (\Delta \eta ({\mu }^{+},{\mu }^{-})))}-m_{\tau }^2}.$$

$$|{\mathcal{M}}_{\mathrm{signal}}{|}^2=|{\mathcal{M}}_{SM}+{\mathcal{M}}_{Z^{\prime }}{|}^2$$

$$\supset |{\mathcal{M}}_{Z^{\prime },\text{ on-shell}}{|}^2+2|{\mathcal{M}}_{WW}^{\ast }{\mathcal{M}}_{Z^{\prime }}|,$$

$$S$$

$$B$$

$$p^{\nu s_i}=\left(\sqrt{{\left(p_T^{\nu s_i}\right)}^2+m_{\nu s_i}^2}\cosh {\eta }_{\nu s_i},p_T^{\nu s_i}\cos {\varphi }_{\nu s_i},p_T^{\nu s_i}\sin {\varphi }_{\nu s_i},\sqrt{{\left(p_T^{\nu s_i}\right)}^2+m_{\nu s_i}^2}\sinh {\eta }_{\nu s_i}\right),$$

$$p_T^{\nu s_i}$$

$$m_{\nu s_i}$$

$${\eta }_{\nu s_i}$$

$${\varphi }_{\nu s_i}$$

$$|{\mathcal{M}}_{\tau \tau }^{mod}{|}^2$$

$$-\frac{8m_{\nu s_1}[(p^Z\cdot p^{{\mu }_1})(p^Z\cdot p^{{\mu }_2})+2(p^Z\cdot p^{\nu s_1})(p^Z\cdot p^{{\mu }_2})]}{m_{\tau }{M}_Z^2}-\frac{2m_{\nu s_1}[(p^{{\mu }_1}\cdot p^{{\mu }_2})+2(p^{\nu s_1}\cdot p^{{\mu }_2})]}{m_{\tau }}$$

$$-\frac{4m_{\nu s_2}[(p^Z\cdot p^{{\mu }_1})(p^Z\cdot p^{{\mu }_2})+2(p^Z\cdot p^{{\mu }_1})(p^Z\cdot p^{\nu s_2})]}{m_{\tau }{M}_Z^2}-\frac{2m_{\nu s_2}[(p^{{\mu }_1}\cdot p^{{\mu }_2})+2(p^{{\mu }_1}\cdot p^{\nu s_2})]}{m_{\tau }}$$

$$+\frac{1}{M_Z^2}[4(p^Z\cdot p^{{\mu }_1})(p^Z\cdot p^{\nu s_2})+2(p^Z\cdot p^{{\mu }_1})(p^Z\cdot p^{{\mu }_2})+8(p^Z\cdot p^{\nu s_1})(p^Z\cdot p^{\nu s_2})]$$

$$\frac{4}{M_Z^2}(p^Z\cdot p^{\nu s_1})(p^Z\cdot p^{{\mu }_2})+2(p^{{\mu }_1}\cdot p^{\nu s_2})+(p^{{\mu }_1}\cdot p^{{\mu }_2})+4(p^{\nu s_1}\cdot p^{\nu s_2})+2(p^{\nu s_1}\cdot p^{{\mu }_2}))$$

Full text

I Introduction

The highest priority of LHC– run II is to find physics beyond the standard model (BSM). In the pursuit of optimal sensitivity to BSM physics, it is important to develop a diverse set of techniques that improve our sensitivity to BSM signal against the SM background. The most straightforward technique is the cut-and-count method, where one determines the suitable region of phase space by studying the kinematic distributions of the simulated events (e.g. the invariant mass of two detected particles, or the angle difference between them) and look for an excess over the SM expectation. While simple, this method can miss correlations among kinematic variables. To improve our discriminatory power, one can select events based on more complicated variables that take into account correlations. These complicated variables can be a linear or a non-linear combination of the kinematic variables used in the cut-count method. The techniques that consider combinations of kinematic variables are generally called “Multivariate Analysis” (MVA) Hocker:2007ht . Some of the main methods developed in this category are “Neural Network” Denby:1987rk ; Bellantoni:1991ei ; Bortoletto:1991ty ; Babbage:1993ry ; Becks:1993jf and “Boosted Decision Tree” Roe:2004na ; Yang:2005nz . However, these MVA methods usually require a phase of computer training and techniques that are not transparent to theorists. Another example of MVA that is calculated based on the theoretical assumptions for a given process is the Matrix Element Method (MEM) Fiedler:2010sg ; Artoisenet:2008zz ; Mattelaer:2011ywa ; Mertens:2014iya ; Bianchini:2016yrt ; Volobouev:2011vb ; Biedermann:2016lvg .

The MEM was originally developed in the Tevatron experiment and was successfully applied to measurement of the top quark mass and single top electroweak production Abazov:2004cs ; Abazov:2005zz ; Kondo:1991dw ; Dalitz:1991wa ; Abulencia:2007br ; Aaltonen:2008bd ; Abazov:2008ds ; Abazov:2009ii ; Aaltonen:2010jr ; Fiedler:2010sg ; D0:2016ull ; Bianchini:2016yrt ; Cortiana:2015rca ; Schouten:2014yza ; Gainer:2013iya ; Abazov:2011fc . In the MEM, given a reconstructed event and a theoretical hypothesis, the probability that the event matches the hypothesis is quantified based on the value of the partonic matrix element for the hypothesis processes evaluated using the final state four momenta. Repeating this procedure using multiple hypothesis processes, e.g. a background process and a signal process, and comparing the results, one can quantify how ‘signal-like’ or ‘background-like’ a given ensemble of events is. By definition, the MEM contains all of the kinematic information of the hypothesis process so it captures all correlations. Moreover, the MEM has a clear physical meaning and there is a transparent link between the theoretical assumptions and event reconstruction.

Despite the successes of the MEM in the SM measurements, this method has not been applied extensively in BSM searches. Some of the main difficulties of the MEM are the following Fiedler:2010sg ; Volobouev:2011vb :

• •

The squared matrix element depends on all of the momenta – both initial and final – in the event. In circumstances where one or more momenta is not determined uniquely, we must integrate over all possible values. The more integrations required for each event, the more time and computational resources required. Events with missing energy, a hallmark of many new physics scenarios, automatically fall into this category of events with unconstrained momenta.

• •

BSM hypotheses necessarily introduce new parameters ($\alpha=$ masses, couplings of new particles), which need to be determined before we can calculate the matrix element. Without a separate experiment telling us what $\alpha$ to use, our best option is to choose the $\alpha$ such that they maximize our reach for detecting the signal given the background. However, finding the best-suited values for these parameters can be difficult, as shown in Section IV.1.

Ideally, we would like to develop techniques to overcome the aforementioned challenges for any BSM hypothesis. However, as a first step, we focus on dealing with some of these challenging for a specific model. The BSM hypothesis we will consider in this paper is one of the simplest extensions of the SM, the $L_{\mu}-L_{\tau}$ model He:1991qd . In the $L_{\mu}-L_{\tau}$ model, the difference between muon number and tau number – an accidental symmetry of the SM – is gauged and spontaneously broken, generating a massive $Z^{\prime}$ that interacts only with the second and third generation of leptons at tree level He:1991qd . As the $L_{\mu}-L_{\tau}\,Z^{\prime}$ is hadrophobic and does not interact with electrons at tree level, it is difficult to produce at conventional experiments and therefore is challenging to constrain. In fact, tree-level $Z^{\prime}$ production at the LHC always involves four second/third generation lepton (either charged or neutral) final states: the initial partons create a pair of leptons via $W^{\pm}/Z/\gamma$ exchange, one of which radiates a $Z^{\prime}$ that subsequently decays to a second lepton pair. At the LHC, identifying muons is much easier and more accurate than identifying taus, thus in this work, we will only consider combinations of muons and $E_{T}{/}\hskip 9.95863pt$(due to $\nu_{\mu}$ or $\nu_{\tau}$) as potential final states. In particular, we will study the two processes of i) four muons ($4\mu$), and ii) two muons and missing energy ($2\mu+\not{E}_{T}$). Of course, as all of the final state particles are SM particles, there will be interference between the four muon events produced via $Z^{\prime}$ and SM four muon events. We will put this interference to use in sections IV.1 and IV.2.

The four muon final state is fully reconstructable, and it has 12 independent observables. The existence of this many kinematic variables begs for using an MVA. The MEM is a suitable choice because it optimally Fiedler:2010sg uses all of the available kinematic variables to distinguish signal from the background. Applying MEM to the process $2\mu+\not{E}_{T}$ with 6 observables is also more lucrative, however, the existence of $E_{T}{/}\hskip 9.95863pt$makes its computation more challenging. Fortunately, the main SM background to this process ($pp\to\tau\tau\to 2\mu\,\not{E}_{T}$) has a very specific topology, and thus the missing momenta for this process can be estimated from the momenta of muons. In section IV.2, we show that by using the squared matrix element of the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ background as a discriminatory variable (as opposed to the conventional likelihood function) with the guessed missing momenta derived from its topology (i.e. without having to integrate over missing momenta), we are able to sufficiently differentiate our signal from the remaining SM backgrounds as well. Due to the deviation from the canonical MEM procedure, we will refer to this approach as the “modified MEM”.

The organization of the rest of the paper is as follows. In the next section, we explain the Matrix Element method and how it is used to discriminate the signal from the background. In section III, we introduce the model and discuss some of the constraints on its parameters from LHC and other experiments. In section IV, we explore how the bounds can be improved at the LHC. Finally, a discussion about the results and concluding remarks are made in section V.

II Matrix Element Method

At the LHC, we are searching for BSM physics against the SM background. Although we could look for generic departure from the SM, our sensitivity is greater if we look for a particular BSM physics hypothesis. Therefore, we are usually dealing with two hypotheses: a specific new physics model (signal) and the null hypothesis (the SM background). To statistically analyze these two test hypotheses, the MEM uses their likelihood ratio, where the likelihoods are calculated based on the squared matrix element ($|\mathcal{M}|^{2}$) of a process at tree level111 Some papers have studied the MEM at NLO as well Baumeister:2016maz ; Martini:2015zkx ; Martini:2015fsa ., and evaluated using the final state momenta of individual events. Since all the kinematic information of a process is contained in the matrix element, $|\mathcal{M}|^{2}$ is a powerful variable to discriminate between the two test hypotheses. If the empirical LHC events are inconsistent with either of the conjectures, $|\mathcal{M}|^{2}$ should also favor neither of them.

The function $|\mathcal{M}|^{2}$, depends on the momentum of initial and final state partons ($p^{\rm par}$). The squared matrix element of the signal $|\mathcal{M}|^{2}_{\text{signal}}$ will additionally depend on the new model parameters ($\alpha$): $|\mathcal{M}|^{2}_{\text{signal}}=|\mathcal{M}(p^{\rm par};\alpha)|^{2}$, and $|\mathcal{M}|^{2}_{SM}=|\mathcal{M}(p^{\rm par})|^{2}$. In the special case where the final states of the signal exactly match the SM background (e.g, no BSM particles in the final state), $|\mathcal{M}|^{2}_{\text{signal}}$ must include the interference term with the SM as well.

If the LHC could detect the $p^{\rm par}$ of all of the final state particles, the likelihood $(P)$ that an observed event is due to a particular hypothesis would be defined as

[TABLE]

where $x_{i}$ and $p_{\text{initial}}^{\text{par}}$ are intimately related: $p_{\text{initial,i}}^{\text{par}}\equiv\frac{\sqrt{s}}{2}(x_{i},0,0,\pm x_{i})$. The factors $f(x_{i})$ are the parton distribution functions (PDF) of the initial states, the collider center of mass energy of collision is represented by $\sqrt{s}$, and $\sigma$ is the total cross section with which $P(p^{\rm par}|\alpha)$ is normalized to ensure $\int P(p^{\text{par}}|\alpha)dp^{\text{par}}=1$. The factor $\delta^{4}(p_{\text{initial}}^{\text{par}}-p_{\text{final}}^{\text{par}})$ ensures the conservation of energy and momentum in the process. If we know the final state four-momenta, we can use this delta function to infer information about the initial state momenta. More specifically, the delta function involving the $p_{\text{initial}}^{\text{par}}$ can be translated to a delta function on $x_{i}$’s, uniquely defining the $x_{i}$’s and collapsing the integrals:

[TABLE]

The events at the LHC, however, are defined according to reconstructed momenta ($p^{\rm rec}$) at the detector which may not equal the $p^{\rm par}$. Hence, we must modify the likelihood $(P)$ to be a function of detector level momenta:

[TABLE]

with an integration over the possible values of partonic momenta given the reconstructed momenta (represented by $\int d\Phi(p^{\rm par}_{\text{final}})dx_{{}_{1}}dx_{{}_{2}}$). When there is missing energy (e.g, neutrino) in the event, we integrate over all unconstrained momenta. However, if the particle is detected, we use the transfer functions $W(p^{\rm rec},p^{\rm par})$ to translate between the partonic momenta and reconstructed momenta. These functions are usually gaussian (or bi-gaussian), where the arguments are estimated according to Monte Carlo (MC) simulations222In the matrix element technique, the theoretical assumptions (the $|\mathcal{M}|^{2}$) and the assumptions about the experiment and detectors (the transfer functions) factorize and are independent of each other. Therefore, the improvement in any of these assumptions can easily be implemented in the matrix element method. Artoisenet:2008zz ; Mattelaer:2011ywa ; Mertens:2014iya . If the detected particle is a lepton or a photon, the transfer function are well approximated by a delta function $(\delta(p^{\text{par}}-p^{\text{rec}})$), while for colored objects, the reconstructed momenta may be significantly different from the partonic momenta.

After calculating the likelihood function ($P$) of the signal and the SM, we need to determine whether a given event is more likely to be due to the signal hypothesis or the background. Therefore, we study the likelihood ratio:

[TABLE]

where the $i$ subscript refers to the $i$-th event. If the value of $\alpha$ were known, we could plot the distribution of $\mathcal{L}$ for given events, just like any usual kinematic distribution that uses the reconstructed momenta. In the $\mathcal{L}$ distribution, larger values would indicate the signal hypothesis is favored and lower values meant the data is more consistent with the null hypothesis. Hence, with a cut on the larger values of $\mathcal{L}$, we could find the phase space that increases the signal fraction.

Having said that, in the case of new physics, we do not know the value of $\alpha$. Since we want to optimize our reach for the signal hypothesis, we want to choose $\alpha$ that maximizes $\mathcal{L}(p_{i}^{\rm rec};\alpha)$ (or equivalently $log\left[\mathcal{L}(\alpha)\right]$). Therefore, we plot $log\left[\mathcal{L}(\alpha)\right]$ with respect to $\alpha$ and look for maxima in the plot. This process can be done for each event. For multiple events, we simply sum over $i$:

[TABLE]

We denote $\alpha^{*}$ for the value of $\alpha$ that maximizes Eq. (5). As it will be important later on, we emphasize that Eq. (5) only yields one number for an entire set of events. Previous studies have shown that if $\alpha$ is the mass of a particle, the $\alpha^{*}$ returned by maximizing Eq. (5) is actually the same as the true value of $\alpha$; in fact, the current most precise measurement of the top quark mass is obtained with this process of maximizing the likelihood ratio with respect to top mass Fiedler:2010sg ; Bianchini:2016yrt . This procedure is used only to determine $\alpha^{*}$, and it gives no information on the prospect of discovering the signal. Once the optimal value $\alpha^{*}$ is determined, the sensitivity to the signal is determined by selecting a region of the $\mathcal{L}(\alpha^{*})$ (with $\alpha=\alpha^{*}$ in Eq. (4)) distribution that optimizes the signal over background ratio.

Despite the power of the MEM in discriminating signal against the null hypothesis, due to its computational difficulties, it is not commonly used in BSM searches. One reason may be the number of phase space integrations in the calculation of $P(p^{\rm rec}|\alpha)$. If the process of interest has multiple sources of missing energy or contains colored objects, the number of integrals can be high and may over-consume computational resources. Another reason is that the maximization of likelihood ratio with respect to $\alpha$ can be challenging depending on the nature of parameters. As we will see in Section 3.A, the maximization of log likelihood works if $\alpha$ is the mass of a new particle, but this approach breaks down if $\alpha$ is a coupling. That is because $log\left[\mathcal{L}(\alpha)\right]$ only increases with respect to $\alpha$ and has no local maxima. Although one naively might expect that larger values of $\alpha$ results in greater sensitivity to BSM hypothesis, in section IV.1, we will show that is not correct.

In this paper, we will study a simple model to give a working example of how some of the challenges in MEM for BSM can be overcome. The example we will consider is the gauged $U(1)_{L_{\mu}-L_{\tau}}$ symmetry, summarized in the next section. This model contains a massive $Z^{\prime}$ gauge boson that only interacts with the second and third generation of leptons. We will show that the LHC sensitivity can be improved by an order of magnitude compared with the cut-and-count method if we apply MEM.

III $L_{\mu}-L_{\tau}$ model

The difference between muon number and tau number $L_{\mu}-L_{\tau}$ is one of the accidental global symmetries present in the SM. This $U(1)_{L_{\mu}-L_{\tau}}$ symmetry is anomaly free and therefore can be gauged. However, from the oscillation of tau or muon neutrinos to electron neutrinos, we know this symmetry is not respected in nature and needs to be broken. The consequence of the breaking is a neutral, color singlet, massive $Z^{\prime}$ that couples only to muon number and tau number at tree level. The interactions of $Z^{\prime}$ are described by the Lagrangian below:

[TABLE]

where $Z^{\prime}_{\alpha\beta}=\partial_{\alpha}Z^{\prime}_{\beta}-\partial_{\beta}Z^{\prime}_{\alpha}$ is the field strength tensor, and $\ell_{2}=(\nu_{\mu},\mu)^{T}$, $\ell_{3}=(\nu_{\tau},\tau)^{T}$. As shown in the Lagrangian, the $Z^{\prime}$ has the same coupling to left handed and right-handed muon (tau), with a relative minus sign between the coupling of muons and taus He:1991qd . The new parameters in the model are $\alpha=(M_{Z^{\prime}},\epsilon)$, where $\epsilon g_{1}$ is the coupling of $Z^{\prime}$ to muons and taus.

Studying this model is important because some region of its parameter space can explain the long persisting discrepancy in the SM prediction and experimental measurement of muon anomalous magnetic moment $(g-2)_{\mu}$ Anastasi:2015oea ; Baek:2001kca ; Ma:2001md ; Gninenko:2001hx ; Pospelov:2008zw ; Heeck:2011wj ; Harigaya:2013twa ; Chen:2017cic . Some anomalies observed in B physics and flavor changing Higgs coupling Crivellin:2015lwa ; Crivellin:2015mga ; Heeck:2014qea can also be explained by gauged $L_{\mu}-L_{\tau}$, which further motivates studying this model. In particular, the anomalies recently observed in $R_{K}=\text{Br}(B\to K\mu^{+}\mu^{-})/\text{Br}(B\to Ke^{+}e^{-})\simeq 0.745$ with $2.6\,\sigma$ discrepancy between theoretical expectations, and in $R_{K^{*}}=\text{Br}(B\to K^{*}\mu^{+}\mu^{-})/\text{Br}(B\to K^{*}e^{+}e^{-})\simeq 0.7$ with $\sim 2.5\,\sigma$, by LHCb Bifani:2260258 can be explained by the $L_{\mu}-L_{\tau}$ model, assuming $M_{Z^{\prime}}\sim O(\text{TeV})$ and $g_{Z^{\prime}}\sim O(1)$ Crivellin:2015lwa ; Crivellin:2015mga ; Alonso:2017bff .

Because the $L_{\mu}-L_{\tau}\ Z^{\prime}$ only interacts with second and third generation leptons at tree level, it is not very constrained. One constraint comes from $Z-Z^{\prime}$ mixing that arises from loops of muons and taus (and their respective neutrinos), inducing a coupling of $O(10^{-3}\,\epsilon)$ between the $Z^{\prime}$ and all fermions. The factor of $10^{-3}$ is a rough estimate based on the loop suppression and the couplings of $Z^{\prime}$ and $Z$ with muon and tau. Precision measurements of $Z-$electron coupling at the BaBar and Belle II experiment Eigen:2015rea ; Curtin:2014cca ; Eigen:2015sea ; Essig:2013vha ; Wang:2015hdf ; Araki:2017wyg requires the coupling of $Z^{\prime}$ to electrons to be $\lesssim 10^{-3}$ for $10\ \text{MeV}<M_{Z^{\prime}}<10\ \text{GeV}$, which translates to $\epsilon\lesssim 1$ in the $L_{\mu}-L_{\tau}$ model.

The strongest current bound on $L_{\mu}-L_{\tau}$ for $M_{Z^{\prime}}<10\ \text{GeV}$ is from fixed target neutrino beam experiments. In particular, CHARMII Geiregat:1990gz and CCFR Mishra:1991ws ; Mishra:1991bv tightly constrain $L_{\mu}-L_{\tau}\ Z^{\prime}$ via the trident process: $N+\nu_{\mu}\rightarrow N+\nu_{\mu}+\mu\mu$ Altmannshofer:2014cfa . This process occurs through the exchange of off-shell $W^{\pm}/Z$ bosons in the SM. However, in the $L_{\mu}-L_{\tau}$ model, the exchange of $Z^{\prime}$ can significantly contribute to the rate of the process, especially if the $Z^{\prime}$ is produced on-shell  Geiregat:1990gz ; Mishra:1991ws ; Altmannshofer:2014cfa ; Kaneta:2016uyt . Neutrino trident experiment excludes a $L_{\mu}-L_{\tau}\ Z^{\prime}$ with $\epsilon\gtrsim 0.005$ and $M_{Z^{\prime}}\lesssim 1\ \text{GeV}$. The bounds loosen for heavier $Z^{\prime}$ to $\epsilon\sim 0.05$ for $M_{Z^{\prime}}=20\ \text{GeV}$.

LHC bounds on this model come from from recasting the $pp\to Z\to 4\mu$ searches by CMS and ATLAS CMS:2012bw ; Aad:2014wra . This bound surpasses the trident bound for $10\ \text{GeV}<M_{Z^{\prime}}\lesssim 45\,\text{GeV}$ Altmannshofer:2014cfa . In Ref.  Elahi:2015vzh , we discussed the potential LHC reach using a dedicated cut-and-count $Z^{\prime}$ analysis in the $pp\to Z\to 4\mu$ (for $M_{Z^{\prime}}>2m_{\mu}$) and $pp\to\mu^{+}\mu^{-}\not{E}_{T}$ (for $M_{Z^{\prime}}<2m_{\mu}$) channels. These channels were proposed for their cleanliness. Additionally, since the contribution of $Z^{\prime}$ is greatest when it is produced on-shell, the mass ranges of the channels were chosen such that an on-shell $Z^{\prime}$ can decay to muons ($pp\to 4\mu$) or neutrinos ($pp\to 2\mu\not{E}_{T}$). Although in Elahi:2015vzh a large region of parameter space could be uncovered with the cut-and-count method after the full $3\ \rm ab^{-1}$ of HL-LHC run, in this paper we will show that our sensitivity can be enhanced further if we use the matrix element technique.

IV LHC bounds on $L_{\mu}-L_{\tau}$ model using Matrix Element Method

IV.1 Looking for $Z^{\prime}$ with mass range $2m_{\mu}<M_{Z^{\prime}}<M_{Z}$ in $pp\rightarrow Z\rightarrow 4\mu$

The rare process $Z\to 4\mu$ occurs through $Z$ boson decay into two muons, one of which radiates a neutral boson $V$ that subsequently splits into a second pair of muons. In the SM, $V$ can be an off-shell $Z$ or photon, while in the $L_{\mu}-L_{\tau}$ model, the on-shell/off-shell $Z^{\prime}$ will also contribute.

Four muons reconstructing an on-shell $Z$ boson is a clean process that has been studied extensively. Therefore, it is an ideal channel for constraining $L_{\mu}-L_{\tau}$ model. Moreover, the related channels $pp\to Z\to 4e$ and $pp\to Z\to 2e\ 2\mu$ can be used as a background control sample to suppress the experimental uncertainties of this channel. Within the SM the cross section of $Z\to 4\mu$ is the same as the $Z\to 4e$ and $Z\to 2e\ 2\mu$ processes up to $O(m_{e}^{2}/m_{\mu}^{2})\sim 10^{-4}$, hence by measuring $Z\to 4e$ and $Z\to 2e\,2\mu$, we can obtain a precise prediction of $(pp\to Z\to 4\mu)_{SM}$. As such, in the following, we will assume that the systematic uncertainties of $(pp\to Z\to 4\mu)_{SM}$ are sub-percent.

To study the $L_{\mu}-L_{\tau}$ $Z^{\prime}$, we generated a Universal FeynRules Output (UFO) model Degrande:2011ua using Feynrules Alloul:2013bka . We then fed the model to MadGraph Alwall:2011uj to generate our events at leading-order. The background sample only contains the SM contributions to the four muon production $(pp\to 4\mu)_{SM}$, while the signal event sample include both the $Z^{\prime}$ as an intermediate state (where its width is calculated using MadGraph), and the SM gauge bosons $(pp\to 4\mu)_{SM+Z^{\prime}}$ to capture the interference among processes. The signal events were generated for various values of $\alpha^{\text{gen}}=(M_{Z^{\prime}}^{\text{gen}},\epsilon^{\text{gen}})$ between $2\,m_{\mu}\leq M_{Z^{\prime}}^{\text{gen}}\leq 40\,\text{GeV}$ and $0.001\leq\epsilon^{\text{gen}}\leq 0.1$.

Before applying MEM, we impose some preliminary cuts to ensure the events have been triggered upon. Specifically, we impose a di-lepton trigger used in LHC-13 Khachatryan:2016fll that selects events with $p_{T}\ (\mu_{1})>17\ \text{GeV}$ and $p_{T}\ (\mu_{2})>8\ \text{GeV}$, where $\mu_{1}$ is the leading muon and $\mu_{2}$ is the sub-leading one. We also require $p_{T}>4\ \text{GeV}$ for all muons and the separation $\Delta R_{\mu\mu}>0.05$  Khachatryan:2016fll ; CMS:2012bw ; Aad:2014wra . To ensure $\sqrt{\hat{s}}\sim M_{Z}$, we impose $76\ \text{GeV}<m_{4\mu}<106\ \text{GeV}$ and veto extra jets in the event, since they would hurt the cleanliness of the process.

As lepton momenta are accurately measured at the detector, we can be confident that the detected momenta very closely represents the partonic level momenta; stated in terms of transfer functions introduced in Eq. (3), we will assume $W(p^{\rm rec},p^{\rm par})=\delta(p^{\rm rec}-p^{\rm par})$. Furthermore, because the final states are fully reconstructable and we have vetoed jets333 we have done our analysis at parton level, but given $m_{4\mu}\sim M_{Z}$, it is reasonable to assume jet contamination is negligible. , we can find the initial state energies (or equivalently the $x_{i}$s in Eq. (1)) from conservation of energy and momentum, as shown in Eq. (2). As a result, all of the integrations collapse due to delta functions, and the likelihood ratio becomes

[TABLE]

Because the initial state quarks (with energies parameterized by $x_{{}_{1}}$ and $x_{{}_{2}}$) are the same for both the signal and the background at leading order, the PDF and $x_{i}$ factors cancel in the ratio. Hence, the likelihood ratio for a single event can be simplified to

[TABLE]

The analytical calculation of $|\mathcal{M}|^{2}$ can be simplified if we just consider the process $Z\to 4\mu$, where $Z$ is the vector sum of the four muons. Although the events could also be due to the processes $qq\to\gamma^{*}\to 4\mu$ or the interference of $\gamma^{*}$ and $Z$ mediators, requiring $76\ \text{GeV}<m_{4\ell}<106\ \text{GeV}$ assures us $|\mathcal{M}_{Z\to 4\mu}|^{2}$ is a reasonably accurate description of the events. For the analytic evaluation of $|\mathcal{M}|^{2}$, we use the form given in MCFM Campbell:1999ah ; Campbell:2011bn ; Boughezal:2016wmq , which is easily modified to include $Z^{\prime}$ intermediate states (for the case of the signal).

Even when all final states particles are reconstructed, there are complications in evaluating $|\mathcal{M}|^{2}$. Specifically, at the detector level, we do not know the ‘right’ pair of muons that reconstruct the $V=\gamma/Z/Z^{\prime}$, leading to a combinatorics problem. The most naive way we can account for this problem is by summing over all four possible $|\mathcal{M}|^{2}$ with different combinations of muon pairs originating from $V$. With a better algorithm, the sensitivity to signal may be further improved, but we do not attempt that here.

The likelihood ratio, as shown in Eq. (7) can be calculated for each event as a function of $M_{Z^{\prime}},$ and $\epsilon$. Before proceeding, we need to differentiate between the analysis $\alpha=(M_{Z^{\prime}},\epsilon)$ values used in the signal matrix element hypothesis ($\mathcal{M}(p^{\text{rec}},M_{Z^{\prime}},\epsilon$)) and the “truth” values – the values of $M_{Z^{\prime}}$ and $\epsilon$ used to generate the signal events (which, in actual data, would be unknown). To avoid confusion, we will use $\alpha^{a}=(M_{Z^{\prime}}^{a},\epsilon^{a}$) for the analysis values and $\alpha^{\text{gen}}=(M^{\text{gen}}_{Z^{\prime}},\epsilon^{\text{gen}})$ which were used for event generation.

Our goal is to determine the values of $M_{Z^{\prime}}^{a}$ and $\epsilon^{a}$ that maximize our sensitivity to the signal hypothesis. Keeping this task in mind, it is worth looking at how the $\alpha^{a}$ parameters enter into the matrix element:

[TABLE]

where $\Gamma_{Z^{\prime}}^{a}=\Gamma_{Z^{\prime}}(M_{Z^{\prime}}^{a},\epsilon^{a})$ and the ellipses represent the part of the matrix element that is independent of $M_{Z^{\prime}}^{a},\epsilon^{a}$. Knowing the dependence of $\mathcal{M}$ on $\alpha^{a}$ can help us better understand the behavior of the likelihood ratio for various $\alpha^{a}$ values. For example, the appearance of $M_{Z^{\prime}}^{a}$ only in the denominator suggests that the likelihood ratio is very sensitive to the value of $M_{Z^{\prime}}^{a}$ (e.g, the possibility of resonance). With $\epsilon^{a}$ in the numerator, we suspect the likelihood ratio to not be as sensitive to a particular value of $\epsilon^{a}$, because regardless of which testing hypotheses the events belongs to, their likelihood ratio will increase by increasing $\epsilon^{a}$. Although $\Gamma_{Z^{\prime}}^{a}$ in the denominator also depends on $\epsilon^{a}$, its contribution is suspected to be much smaller than the $\epsilon^{a}$ in the numerator.

To examine these conjectures and hunt for the $\alpha^{a}$ that maximize the likelihood ratio, we will proceed by plotting $\sum_{i}\log[\mathcal{L}(p_{i},\alpha^{a})]$ defined in Eq. (5) as a function of $\alpha^{a}$. Thereby, for each $\alpha^{a}$ we can compare the relative sizes of the summed log likelihood ratio of the signal sample events with that of background sample, and look for the values of $\alpha^{a}$ that maximize the likelihood of the signal with respect to the background. As we have two new parameters ($\alpha^{a}=M_{Z^{\prime}}^{a},\epsilon^{a}$), for simplicity we will fix one of the parameters and plot Eq. (5)) as a function of only one parameter.

First, we assume a non-zero value for $\epsilon^{a}$ and plot Eq. (5) for a range of $M_{Z^{\prime}}^{a}$ values. The value of $\epsilon^{a}$ can be any arbitrary non-zero value, and varying $\epsilon^{a}$ does not alter our results. For every event, the reconstructed invariant mass of two muons ($M_{\mu\mu}$) will produce a spike at $M_{Z^{\prime}}^{a}=M_{\mu\mu}$, which is not necessarily equal to the $M_{Z^{\prime}}^{\text{gen}}$. In fact, because we are summing over all four possible combinations of muon pairs, we have four spikes in $M_{Z^{\prime}}^{a}$ for each event. Therefore, summing over all events, the net distribution (in our case, 50000 events) has spikes at all values of $M_{Z^{\prime}}^{a}\lesssim M_{Z}$ as shown in the left plot of Fig. 1. In the signal sample, due to the higher number of events with $M_{\mu\mu}=M_{Z^{\prime}}^{\rm gen}$, the spike at $M_{Z^{\prime}}^{\rm gen}$ is more noticeable (see Fig. 1, right plots). The maximum likelihood ratio of the signal events is consistently at $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$ for any arbitrary value of $\epsilon^{a}\neq 0$, while the pure background events show no interesting behavior at $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$.

The height of the peak at $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$ relative to the other spikes depends on $\epsilon^{\rm gen}$ (regardless of $\epsilon^{a}$). The peak is more visible for larger values of $\epsilon^{\rm gen}$, but becomes less distinguishable from other spikes for smaller values of $\epsilon^{\rm gen}$. As one example, for $\epsilon^{\rm gen}\lesssim 0.01$ for $M_{Z^{\prime}}^{\rm gen}\gtrsim 10\,\text{GeV}$, we could not distinguish the peak at $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$ from other spikes). Given that our analysis highly depends on whether we are able to find $M_{Z^{\prime}}^{\rm gen}$ by plotting $\sum_{i}\log[\mathcal{L}(p_{i},M_{Z^{\prime}}^{a},\epsilon^{a}_{\text{ arbitrary}}\neq 0)]$ as a function of $M_{Z^{\prime}}^{a}$, our strategy needs to bifurcate depending on the value of $\epsilon^{\rm gen}$:

• •

large $\epsilon^{\text{gen}}$: we determine the value of $M_{Z^{\prime}}^{a}$ by observing the peak at $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$;

• •

small $\epsilon^{\text{gen}}$: we need to chose an arbitrary value of $M_{Z^{\prime}}^{a}$.

However, as the choice of the strategies depends on whether we are able to observe the peak at $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$, rather than trying to quantify what it takes to observe a peak, we will study both strategies for all values of $\epsilon^{\rm gen}$. More specifically, we will first assume that we can find $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$ even for small values of $\epsilon^{\rm gen}$, then we will fix an arbitrary value for $M_{Z^{\prime}}^{a}$ and study the signal assuming we cannot find $M_{Z^{\prime}}^{\rm gen}$ even for large values of $\epsilon^{\rm gen}$. The first method represents the best we can do with the MEM, while the second represents a more realistic reach for $\epsilon\lesssim 0.01$.

Now that we have discussed how to find the optimal value of $M_{Z^{\prime}}^{a}$, we follow the same procedure to determine $\epsilon^{a}$. For simplicity, we fix $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$ and let $\epsilon^{a}$ be a free parameter. The $\sum_{i}\log[\mathcal{L}(p_{i},M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen},\epsilon^{a})]$, where the sum is over $50000$ MC generated events for each sample, are plotted in Fig. 2 below as a function of $\epsilon^{a}$. Both the signal and background MC events increase as a function of $\epsilon^{a}$. Such behavior is anticipated, because as we can see in Eq. (8), larger $\epsilon$ results in larger $|\mathcal{M}(p^{\rm rec};M_{Z^{\prime}},\epsilon)|^{2}$ and thus larger likelihood ratio for any event. Stated another way, fixing $M_{Z^{\prime}}^{a}$ and studying the likelihood ratio with respect to $\epsilon^{a}$ is not useful because the simulated signal and background events behave the same way as a function of $\epsilon^{a}$. Thus, for the remainder of this section, we will explore an alternative different technique for estimating the optimal $\alpha^{a}$.

The problem with the approach of $\mathcal{L}$ maximization is that, for a given $\alpha^{a}$, it compares only two numbers: the summed log likelihood ratio of the simulated signal events vs. that of background events. With this approach, we are throwing away a lot of information about each event. We can get a better handle on the signal vs. the background, if we look at their distributions for a given $\alpha^{a}$ and how those distributions evolve as we vary $\alpha^{a}$. The distribution will illustrate to us the behavior of events for a given $\alpha^{a}$, rather than just their sum. We can use this distribution to find $\alpha^{a}$ that results in best discrimination of the signal from the background, even if the sum over events is similar. The signal usually resides in high values of the likelihood ratio, as expected according to the definition of the likelihood ratio (Eq. (4)). So, we distribute the events in the inverse likelihood ratio ($\mathcal{L}^{-1}$) to be able to see the excess of the signal spread over a narrower window $[0,1]$.444not all of the background will fall in this $[0,1]$ region, however we only care about how the background is distributed in the region where the signal resides.

As an example of $\alpha^{a}$ distribution method, in Fig. 3 below we show the distribution of the MC generated events for the signal and the background for different choices of $\epsilon^{a}$, assuming $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$. We can see that as we increase $\epsilon^{a}$, the resonance region – defined as events with $M_{\mu\mu}\sim M_{Z^{\prime}}^{\rm gen}$ – is more separated from the background. However for very large $\epsilon^{a}$ (for example $\epsilon^{a}=0.5$ in the Fig. 3), the background distribution also becomes more spread. This is again consistent with Eq. (8), because no matter whether the events are at resonance or not, increasing $\epsilon^{a}$ increases their likelihood ratio.

The optimal value of $\epsilon^{a}$ depends on the distribution of the signal with respect to the background. For a fixed sample size of signal and background simulated events, optimal is defined as maximization of

[TABLE]

After examining a few values555The benchmark points $\epsilon^{a}$ we studied are $\epsilon^{a}\in[0.01,0.5]$ with increments of $0.01$. With a larger sample of $\epsilon^{a}$, the best value of $\epsilon^{a}$ may slightly vary. of $\epsilon^{a}$, we find $\epsilon^{a}=0.05$ to give the best discrimination, regardless of the $\epsilon^{\rm gen}$. This result is consistent for all values of $M_{Z^{\prime}}^{\rm gen}$ that we studied.

Thus far, we have seen that with the assumption that we know $M_{Z^{\prime}}^{\rm gen}$, the optimal values for $\alpha^{a}$ are $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\text{gen}}$ and $\epsilon^{a}=0.05$. However, if we cannot find the value of $M_{Z^{\prime}}^{\rm gen}$ by maximizing $\sum_{i}\log[\mathcal{L}(p_{i},\alpha^{a})]$ with respect to $M_{Z^{\prime}}^{a}$, we must explore how our sensitivity changes if $M_{Z^{\prime}}^{a}$ is fixed to an arbitrary value $\neq M_{Z^{\prime}}^{\rm gen}$. Because the $Z^{\prime}$ mass for the generated events is no longer the same as the $Z^{\prime}$ mass used in the analysis matrix element, the events with $M_{\mu\mu}\sim M_{Z^{\prime}}^{\rm gen}$ are no longer at resonance in the matrix element $\mathcal{M}(p^{\rm rec};M_{Z^{\prime}},\epsilon)$ (Eq. (8)). Therefore, when we go to search for the optimal $\epsilon^{a}$ using distributions, the signal, and background are less separated than in the previous case. To increase the separation, we need to increase $\epsilon^{a}$. Consequently, the optimal $\epsilon^{a}$ is no longer fixed at $0.05$ and will depend on the difference between $M_{Z^{\prime}}^{a}$ and $M_{Z^{\prime}}^{\rm gen}$. For example, in Fig. 4, the distribution of the signal and the background MC generated events for $M_{Z^{\prime}}^{a}=25\ \text{GeV}$ and $M_{Z^{\prime}}^{\rm gen}=10\ \text{GeV}$, with $\epsilon^{a}=0.05,\,0.5,\,1$ are shown. In this example, we see $\epsilon^{a}=0.05$ does not give a good discrimination of the signal from the background, and we have to use larger values of $\epsilon^{a}$. Furthermore, as the plot illustrates, although large $\epsilon^{a}$ increases the separation between the signal region and the background, it will broaden the signal region. Because we are interested in distinguishing the region where the signal to background ratio is maximized, having a broad signal region is not ideal. Hence, we expect our reach for $M_{Z^{\prime}}^{a}\neq M_{Z^{\prime}}^{\rm gen}$ to be more limited that the previous case, even with the most optimal $\epsilon^{a}$.

Table 1 recapitulates the value $\epsilon^{a}$ depending on whether or not plotting $\sum_{i}\log[\mathcal{L}(p_{i},\alpha^{a})]$ as a function of $M_{Z^{\prime}}^{a}$ can tell us $M_{Z}^{\rm gen}$.

To quantify the reach of the strategy proposed here, let us define the following variables:

[TABLE]

The cuts on $\mathcal{L}^{-1}$ are imposed such that we get the maximum $S/\sqrt{B}$ with integrated luminosity of $300\ \rm fb^{-1}$. We also require $S>10$ to avoid confusion of signal events for statistical fluctuations. We find that, by using the MEM approach we are able to get $S/\sqrt{B}\geq 3$ up to $\epsilon^{\rm gen}\sim 0.002$ and $M_{Z^{\prime}}\lesssim 20\ \text{GeV}$, and up to $\epsilon^{\rm gen}\sim 0.005$ for $20\,\text{GeV}<M_{Z^{\prime}}^{\rm gen}\lesssim 40\,\text{GeV}$, provided that we can determine the true value of $M_{Z^{\prime}}^{\rm gen}$ by conventional MEM means. This best case scenario is indicated with the dashed red line in Fig. 5. The dashed-dotted gray line in Fig. 5 shows our reach assuming $M_{Z^{\prime}}^{a}=2\ \text{GeV}$ while optimizing $\epsilon^{a}$. Similarly, the dashed brown line is for $M_{Z^{\prime}}^{a}=10\ \text{GeV}$, and the dotted green line is for $M_{Z^{\prime}}^{a}=25\ \text{GeV}$. As expected, the lines with arbitrary $M_{Z^{\prime}}^{a}$ touch the dashed red line for $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$, and have a relatively good sensitivity when $M_{Z^{\prime}}^{a}\,\sim\,M_{Z^{\prime}}^{\rm gen}$, but their sensitivity declines as $M_{Z^{\prime}}^{a}$ moves away from $M_{Z^{\prime}}^{\rm gen}$. The current constraints from CCFR experiment and LHC-8 are also shown in Fig. 5 in black and purple solid lines, respectively. The expected exclusion bound (3 $\sigma$) from LHC-14, with $300\ \rm fb^{-1}$ integrated luminosity, using the cut-and-count method that is explored in Ref. Elahi:2015vzh is the dashed blue line in Fig. 5. With the assumption that we know $M_{Z^{\prime}}^{a}=M_{Z^{\prime}}^{\rm gen}$, our reach with MEM is about a factor of 10 greater than the canonical cut-and-count method, and even with $M_{Z^{\prime}}^{a}$ fixed to an arbitrary value, we have an improved sensitivity for some range of $M_{Z^{\prime}}^{\rm gen}$ compared with the cut-and-count method. By trying several values of $M_{Z^{\prime}}^{a}$ and choosing $\epsilon^{a}\sim 0.05$, one can achieve a sensitivity near the red curve shown in Fig. 5 at the LHC.

IV.2 Looking for $Z^{\prime}$ with mass range $M_{Z^{\prime}}<2m_{\mu}$ in $pp\rightarrow\mu^{+}\mu^{-}$ $E_{T}{/}\hskip 9.95863pt$

If $M_{Z^{\prime}}<2m_{\mu}$, the decay of on-shell $Z^{\prime}$ to muons is kinematically forbidden. In fact, for this mass range, on-shell $L_{\mu}-L_{\tau}\ Z^{\prime}$ can only go to neutrinos. While we could hunt for off-shell $Z^{\prime}$ in this mass range using the $4\mu$ final state (as in the previous section), off-shell $Z^{\prime}$ production is suppressed relative to on-shell production by two additional powers of $\epsilon$. Therefore, as proposed in Elahi:2015vzh , we will give up the benefits of a completely visible final state in favor of larger rate and hunt for $Z^{\prime}$ in $pp\to\mu^{+}\mu^{-}\not{E}_{T}$. The signal contribution to $pp\to\mu^{+}\mu^{-}+\not{E}_{T}$ comes from the production of a pair of muons, one of which radiates an on-shell $Z^{\prime}$ that decays to a pair of neutrinos 666The signal also captures the diagram where the neutrinos and muons are produced in opposite order: a pair of neutrinos are first produced, and then one of the neutrinos emits a $Z^{\prime}$ which splits to two muons. In this topology, the $Z^{\prime}$ must be off-shell, so the cross section is suppressed by more powers of $\epsilon$ and is, therefore, negligible in the parameter space of interest.: $pp\to\mu^{+}\mu^{-}Z^{\prime}\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell}$ .

Due to the presence of missing energy, we can no longer impose an invariant mass cut on the final state. As a consequence, some new, important backgrounds emerge:

[TABLE]

The last background arises as a result jet mis-measurements and pileup, and so peaks at low values of $E_{T}{/}\hskip 9.95863pt$. Di-tau production is the largest irreducible background, followed by diboson production. Technically, we include both resonant and non-resonant contributions in this category, as the latter can be non-negligible. Thus, the $VV$ background is more accurately described as $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM}$. Similarly, the ‘signal’ in this section is defined as $pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell}$ including $Z^{\prime}$ as a possible intermediate state: $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM+Z^{\prime}}$. As in the previous section, the signal is defined including SM contributions to incorporate interference.

To study this channel in more detail, we generated events for $pp\to\tau^{+}\tau^{-}\to\mu^{+}\mu^{-}+\not{E}_{T}$, $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM}$, and $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM+Z^{\prime}}$ via the MC chain MadGraph5-aMC@NLO Alwall:2011uj plus Pythia 6.4  Sjostrand:2006za , where the latter step is used here to decay the taus.777We restricted the possible $\tau$ decays to leptonic channels only within Pythia 6.4 to make event generation more efficient Before any MEM analyses, we require all events to pass the same dilepton trigger requirement used in Sec. IV.1 ($p_{T}\ (\mu_{1})>17\ \text{GeV}$ and $p_{T}\ (\mu_{2})>8\ \text{GeV}$, where $\mu_{1}$ is the leading muon and $\mu_{2}$), veto any jets with $p_{T}>20\,\text{GeV},\eta<2.5$, and impose a minimum missing energy cut of $\not{E}_{T}>20\,\text{GeV}$. The last cut is imposed to suppress the $pp\to\mu^{+}\mu^{-}+\text{jets}$ background.

After the initial set of cuts, muonic tau production $pp\to\tau^{+}\tau^{-}\to\mu^{+}\mu^{-}\not{E}_{T}$ is our main irreducible background. The cross section of the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ background is roughly two orders of magnitude larger than $\sigma(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM}$. Therefore, if we want to have any chance to be sensitive to a $Z^{\prime}$ signal, we need to first make the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ background more manageable. So, instead of trying to discriminate $Z^{\prime}$ signal against background, we will focus on distinguishing $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ from other processes.

The most efficient way to eliminate the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ background is by using a variable that is most faithful to $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ background and thus localizes its simulated events to a small region. Motivated by the benefits of the MEM discussed previously, we will use the $|\mathcal{M}_{\tau\tau}|^{2}$, the squared matrix element of the muonic di-tau production, as the discriminating variable. For each event we will calculate $|\mathcal{M}_{\tau\tau}|^{2}$ using the observed final state momenta, then search for and select out regions (using MC) of $|\mathcal{M}_{\tau\tau}|^{2}$ that $pp\to\tau^{+}\tau^{-}$ does not populate. In doing this, we are not following the traditional MEM in this section, because we are not using $|\mathcal{M}|^{2}$ of the signal to distinguish that from other processes. Rather, we are only using the $|\mathcal{M}|^{2}$ of (part of) the background, which makes our approach independent of the signal ($Z^{\prime}$) model.

As a further deviation from the traditional MEM, we will weight each event by only one squared matrix element, $|\mathcal{M}_{\tau\tau}|^{2}$, rather than two (a ‘signal’ hypothesis and a ‘background’ hypothesis). We may loose some discriminating power by not calculating the likelihood ratio as described in section II, but our approach is more time efficient as we do not have to deal with other squared matrix elements that contain missing energy.

Even after reducing our discriminant to the evaluation of a single $|\mathcal{M}|^{2}$ for each event, evaluating $|\mathcal{M}_{\tau\tau}|^{2}$ is still a difficult task. In $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ production, there are four sources of missing energy, which translates to 12 unknown momenta. Moreover, because we do not know the energies of the initial state, we have in total 14 unknowns. Energy-momentum conservation $\delta^{4}(p_{\text{initial}}-p_{\text{final}})$ reduces the number of unknowns to 10, but 10 integrations for each event is still extremely computationally cumbersome. Thankfully, the specific topology of the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ production can help us approximate the unknowns. The list of our assumptions to approximate the unknowns are the followings:

1. I.

We will assume that the $\tau$s were produced on-shell. We know that the invariant mass of an on-shell tau decay products is the tau mass. Therefore, we can determine two unknowns from this assumptions since there are two $\tau$s in each process. 2. II.

Instead of calculating $\tau\to\nu_{\tau}\bar{\nu}_{\mu}\mu$, we replace the two neutrinos with one massive scalar888Technically we should also consider vector massive neutrinos, though we do not expect this choice affects our results. For $m_{\nu_{s}}=0$ the vector results are identical to the scalar case, while there is a small shift in the matrix element if $m_{\nu_{s}}\neq 0$. neutrino ($\nu$s for the notation), and we calculate $\tau\to\mu\ \nu$s.

Three body decay has different kinematic distributions compared to two body decays, but in this approach we can reduce our number of unknowns by 4, leaving 4 remaining. Because we are no longer dealing with the actual $|\mathcal{M}|^{2}$ and we are calculating $pp\to\tau^{+}\tau^{-}\to\mu^{+}\mu^{-}\nu s\ \bar{\nu s}$, we refer to the matrix element we calculate as the “modified” $|\mathcal{M}_{\tau\tau}|^{2}$, or $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$. 3. III.

The tau pair can be produced from either a photon or a $Z$ boson. However, as a result of our basic cuts (di-lepton trigger and $\not{E}_{T}>20\ \text{GeV}$), we can be confident that the production of $\tau^{+}\tau^{-}$ is dominated by $Z$ exchange. Therefore, we will assume that the taus are produced from an on-shell $Z$, and so $\sqrt{\hat{s}}=M_{Z}$. This assumption eliminates another unknown. 4. IV.

Based on the previous assumptions, we expect the tau decay products to be nearly collinear. Consequently, the $\eta$ and $\phi$ of the $\nu s_{i}$ should be close to $\eta$ and $\phi$ of the corresponding $\mu_{i}$, with the subscript $i$ defined as the following: $\tau_{i}\to\mu_{i}\nu s_{i}$. For our analysis, we assume $\eta(\mu_{i})=\eta(\nu s_{i})$, and $\Delta\phi(\mu_{i},\nu s_{i})=\epsilon_{\phi}^{i}$ where $\epsilon_{\phi}\ll 1$. Therefore, we only work to first order in $\epsilon_{\phi}$. These assumptions specify two more unknowns and results in the relationship $p_{T}^{\tau_{i}}=p_{T}^{\nu s_{i}}+p_{T}^{\mu_{i}}.$ 5. V.

We can also assume $p_{T}^{\tau_{1}}=p_{T}^{\tau_{2}}$. This is the same as assuming there is no initial or final state radiation, which is reasonable given that $\sqrt{\hat{s}}\lesssim M_{Z}$, and we have vetoed jets in our events. This assumption leads to specification of one more unknown.

Making the above approximations 999 Although we have defined “modified” specifically for the second approximation, we generalize it definition to encompass all of the aforementioned approximations. , we can determine all of the unknown kinematic parameters and therefore calculate $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$ with no integrations, significantly reducing the time and computational power needed to do the analysis. The analytical expression of $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$ is given in Appendix A.

We emphasize that these assumptions are only reasonable in the context of the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ production. For the rest of the processes ($(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM}$ background and $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM+Z^{\prime}}$ signal), the approximations I-V are not faithful to the kinematics and we might get unphysical results, i.e. $|\mathcal{M}^{mod}_{\tau\tau}|^{2}<0$. To get an idea of how negative matrix element squared can arise, let us look at the conclusion of assumption IV: $p_{T}^{\tau}=p_{T}^{\nu s_{i}}+p_{T}^{\mu_{i}}$, which means $p_{T}^{\nu s}$ is calculated based on $p_{T}^{\tau}$. We can determine $p_{T}^{\tau}$ using conservation of energy and momentum and a combination of assumptions. We get

[TABLE]

Consequently, the $p_{T}$ of the vector sum of the two neutrinos coming from a tau can be deduced from these approximations: $p_{T}^{\nu s_{i}}=(p_{T}^{\tau})^{\rm approx}-p_{T}^{\mu_{i}}$. These approximations were reasonable in the framework of $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$. However, for $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})$ we can have $(p_{T}^{\tau})^{\rm approx}<p_{T}^{\mu_{i}}$ and therefore a negative (unphysical) magnitude for the transverse transverse momentum, which may lead to negative $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$. The weights (area normalized) of $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$ for MC generaed $(pp\to\tau^{+}\tau^{-}\to\mu^{+}\mu^{-}\not{E}_{T})$, $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM}$ and $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM+Z^{\prime}}$ are shown below in Fig. 6. As expected, the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ distribution is more localized and all of its events have $|\mathcal{M}^{mod}_{\tau\tau}|^{2}>0$. Therefore, if we restrict ourselves to events with $|\mathcal{M}^{mod}_{\tau\tau}|^{2}<0$, we can safely assume that $\tau\tau$ background is negligible. We have generated 10 million $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ events, and 1 million events for each of the $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM}$ and $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM+Z^{\prime}}$ processes, to make sure we have captured the tail of the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ distribution correctly.

Inspecting Fig. 6, we can see that the signal $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM+Z^{\prime}}$ and background $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM}$ have slightly different weights according to $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$. This is not completely surprising because there are contributions from more diagrams in the signal events. The difference is most significant at large, negative $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$; in particular, $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})$ including $Z^{\prime}$ intermediate states populates $|\mathcal{M}^{mod}_{\tau\tau}|^{2}<-4$ more than when the $Z^{\prime}$ is excluded. This $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$ region is shown in greater detail in Figure 7.

To comprehend why $(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM+Z^{\prime}}$ prefers to be in the region $|\mathcal{M}^{mod}_{\tau\tau}|^{2}<-4$, let us look once again at squared matrix element of the signal

[TABLE]

where $\mathcal{M}_{WW}$ refers to the portion of the SM matrix element that involves the $WW$ contribution. Obviously, the departure of the signal from the SM background is most ideal for discrimination when $Z^{\prime}$ is on-shell101010The greatest contribution to the signal (only considering the $Z^{\prime}$ contribution) is when on-shell $Z^{\prime}$ comes from an on-shell $Z$. On the other hand, because both $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ background and signal have $\sqrt{\hat{s}}\lesssim M_{Z}$, a significant fraction of the signal (only the $Z^{\prime}$ contribution) removed when removing $\tau\tau$ background. This is inevitable, and the same challenge was faced with the cut and count method described in  Elahi:2015vzh , when $M_{T}(\mu\mu,\not{E}_{T})<100\ \text{GeV}$ was imposed to remove $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$. Therefore, it is really important to use the interference to look for our signal. . That is because 1) the $|\mathcal{M}|^{2}$ is suppressed by only two powers of $\epsilon$ – only one powers of $\epsilon$ at the production vertex of $Z^{\prime}$ (amplitude level) and no $\epsilon$ suppression at the decay vertex, and 2) the topology of process with an on-shell $Z^{\prime}$ mediator is different from the SM background, and so with some careful cuts we can make the SM background small. The next most important contribution of the signal is in the $Z^{\prime}-SM$ interference, also suppressed by only two powers of $\epsilon$. The interference term is significant when either the SM piece is sizable or when the portion with $Z^{\prime}$ contribution is big. Each option requires different kinematics; large $Z^{\prime}$ contribution means the invariant mass of the neutrinos is small (or equivalently angular separation between the neutrinos is small), while large contribution of the SM could be when some of the intermediate states are produced on resonance. The only SM background that can have on-shell resonances and yet have other kinematics consistent with a (nearly) on-shell $Z^{\prime}$ is the t-channel $W^{+}W^{-}$ background. We suspect the excess in the signal in the region of $|\mathcal{M}^{mod}_{\tau\tau}|^{2}<-4$ is due to the interference of $Z^{\prime}$ piece with the $WW$ contribution.

To check this intuition, we study the distributions of the events in the invariant mass of a muon and the associated neutrino (i.e. $M_{\nu_{\ell}\,\mu^{+}}$ or $M_{\bar{\nu}_{\ell}\,\mu^{-}}$), and the separation between the two neutrinos $(\Delta R(\nu_{\ell},\bar{\nu}_{\ell}))$ after requiring $|\mathcal{M}^{mod}_{\tau\tau}|^{2}<-4$ in Fig. 8. These are not kinematic variables that we could actually measure, as neutrinos are not observed at the detector. However, studying them can help us understand the behavior of the simulated events for different regions of $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$.

We can see from Fig. 8 that in the region of $|\mathcal{M}_{\tau\tau}^{mod}|^{2}<-4$, all of the signal events have $M_{\mu\nu}\sim M_{W}$, and have a small separation between the two neutrinos. This is perfectly consistent with what we expected from the interference term; the effect of the interference is enhanced when $W$s are near resonance, and other kinematic distributions are more consistent with $Z^{\prime}$ signal. Thereby, we can be confident that the excess at $|\mathcal{M}_{\tau\tau}^{mod}|^{2}<-4$ is due to the interference between the signal and the $W^{+}W^{-}$ background.

With the basic cuts, $\not{E}_{T}>20\ \text{GeV}$, and modified $|\mathcal{M}_{\tau\tau}|^{2}<-4$, we get the cross section of the background $\sigma(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM}$ of $80\pm 3$ ab, while $\sigma(pp\to\mu^{+}\mu^{-}\nu_{\ell}\bar{\nu}_{\ell})_{SM+Z^{\prime}}$ is $133\pm 4$ ab. The uncertainties on the cross sections are derived based on the number of events in our simulation that survive the imposed cuts. Let us define:

[TABLE]

If we use the significance $S/\sqrt{B}$ as a test statistic and assume that $pp\to\mu^{+}\mu^{-}+\text{jets}$ (which was our reducible background) is zero 111111This assumption is backed up by a MC study of 500K $pp\to\mu^{+}\mu^{-}$ events of generated at the detector level (PGS PGS ) with the default smearing algorithm. Requiring events pass the dilepton trigger and contain no jets with $p_{T}^{j}>20\,\text{GeV}$ and photons with $p_{T}^{\gamma}>10\,\text{GeV}$, the cross section was 111 pb. After imposing $\not{E}_{T}=p_{T}(\text{dimuon})>20\,\text{GeV}$, the cross section drops to $0.22\,\text{pb}$. With the further requirement of $|\mathcal{M}_{\tau\tau}^{mod}|^{2}<-4$, we find the $pp\to\mu^{+}\mu^{-}$ background can be removed completely. after imposing the $E_{T}{/}\hskip 9.95863pt$cut and $|\mathcal{M}_{\tau\tau}^{mod}|^{2}<-4$, we find that after $300\ \rm fb^{-1}$ luminosity, we get $S/\sqrt{B}\geq 3$.

Our reach in the region of the parameter space using MEM in contrast with the cut-and-count method and the $(g-2)_{\mu}$ band is shown in Fig. 9. The MEM and cut-and-count method bounds are based on the benchmark point of $M_{Z^{\prime}}^{\text{gen}}=0.05\,\text{GeV}$ and $\epsilon^{\text{gen}}=0.001$. We see that with MEM, our sensitivity improves by a factor of 5-10 compared to the cut-and-count method, and we can explore a greater region of parameter space including the $(g-2)_{\mu}$ band.

V Discusssion

In this paper, we applied the Matrix Element method (MEM) to the $L_{\mu}-L_{\tau}$ model and concluded that our sensitivity improves by up to an order of magnitude compared with the cut-and-count method. The MEM uses the matrix element of a process to distinguish signal from background, and it has proven to be a powerful tool in several Standard Model (SM) measurements (e.g, top quark mass) and searches (e.g, Higgs to four lepton and electroweak single top production). However, it has not been extensively applied to beyond the SM (BSM) searches. The two main hurdles are that the MEM needs the physical parameters of the model and the four momenta of all initial and final states as inputs. In the quest of finding BSM signatures, we do not know the value of the new model parameters and the processes often contain unknown momenta in the form of missing energy.

To investigate how we can combat these difficulties, as a first small step, we applied the MEM to the $L_{\mu}-L_{\tau}$ model. This model, being one of the simplest extensions of the SM, is already very well motivated because it can explain some of the current observational anomalies in the $(g-2)_{\mu}$ and B decays. In this model, there is a $Z^{\prime}$ that couples to only second and third generation leptons at tree level. As a result, any tree-level process at the LHC involving $Z^{\prime}$ has to include four leptons of second or third generation. We considered two cases, $4\mu$ and $2\mu+\not{E}_{T}$. Both of these processes have a large number of kinematic observables, making them ideal test grounds for the MEM approach. The new parameters introduced by the $L_{\mu}-L_{\tau}$ model are the $Z^{\prime}$ mass $M_{Z^{\prime}}$ and coupling (parameterized by $\epsilon$: $g_{Z^{\prime}}=\epsilon g^{\prime}$).

In the mass range $2m_{\mu}<M_{Z^{\prime}}<M_{Z}$, we looked at the process $pp\to Z\to 4\mu$. This channel is clean and well-understood, and the presence of $Z\to 4e$ and $Z\to 2e\ 2\mu$ control samples can be used to mitigate systematic uncertainties. Therefore, we can be sensitive to percent-level deviations. As the matrix element of the signal depends on $M_{Z^{\prime}}$ and $\epsilon$, we first discussed how we can find values of these parameters that best separate signal from background. The optimal value of $M_{Z^{\prime}}$ can be found by maximizing likelihood ratio with respect to $M_{Z^{\prime}}$, depending on the strength of the $Z^{\prime}$ coupling. However, the likelihood ratio function increases monotonically as a function of $\epsilon$ regardless of whether events belong to the signal sample or the background sample. Hence, we had to deviate from the conventional MEM and look for the most optimal analysis value of $\epsilon$ by studying the distribution of the signal and background MC generated events as a function of likelihood ratio for various fixed values of $\epsilon$. In the best case scenario, $L_{\mu}-L_{\tau}$ model can be explored up to 3$\sigma$ for $\epsilon\gtrsim 0.002$ for $2m_{\mu}<M_{Z^{\prime}}<20\ \text{GeV}$ and $\epsilon\gtrsim 0.005$ for $20<M_{Z^{\prime}}<40\ \text{GeV}$ with luminosity of $300\ \rm fb^{-1}$ at the LHC, which is about an order of magnitude improvement compared with the cut-and-count method.

For lighter $Z^{\prime}$, we studied the process $pp\to 2\mu\not{E}_{T}$. Due to the presence of missing energy, this channel is not as clean as the all muonic final state and is afflicted by several backgrounds. One significant background is $pp\to\tau^{+}\tau^{-}$ with the taus decaying to muons, which has a cross section that is orders of magnitude greater than signal. To be sensitive to the signal, we first work towards eliminating the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$. We use the squared matrix element of $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ for this task. This is a departure from the canonical MEM, as we weight the events by the squared matrix element of only one process, whereas in MEM we usually weight the events by the ratio of the squared matrix element of the signal processes over the background ones. This alternative approach has two main benefits: 1) given that we have missing energy in the process, calculating the squared matrix elements is challenging; Hence, focusing our attention to only one can save us time and computational power. 2) Furthermore, this approach is independent of the model parameters and can be used for any BSM physics with this signature at the LHC.

Due to the presence of 4 sources of missing energy in the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ process, we have 10 unconstrained momenta in this channel. Traditionally, one would proceed by integrating over the unknown momenta. In this paper, we instead showed how the unknown momenta could be estimated by exploiting the topology of the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$. Without the need for any integrations, we calculate the (modified) squared matrix element of $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ ($|\mathcal{M}_{\tau\tau}^{mod}|^{2}$) and use the resulting weight as a discriminant. With a judicious cut on $|\mathcal{M}_{\tau\tau}^{mod}|^{2}$, we find the $\left.\tau^{+}\tau^{-}\right|_{\text{dimuon}}$ background can be completely eliminated. We then observe that the signal reacts differently to the $|\mathcal{M}_{\tau\tau}^{mod}|^{2}$ compared to other (non-$tau$) SM backgrounds, and trace the difference to interference between $Z^{\prime}$ contributions to the amplitude and contributions containing two on-shell $W$’s. As a result, we can differentiate the signal from all SM backgrounds using $|\mathcal{M}_{\tau\tau}^{mod}|^{2}$ alone. With this method, we find we can reach to 3$\sigma$ up to $\epsilon\sim 0.001$ for $M_{Z^{\prime}}<2m_{\mu}$ assuming an integrated luminosity of $300\ \rm fb^{-1}$, covering the $(g-2)_{\mu}$ band. This result may be improved if we relax some of the kinematic assumptions and instead integrate over a subset of the unconstrained momentum, something which may be worth investigating in the future.

Because we did not use the squared matrix element of the signal in our analysis of $pp\to 2\mu\not{E}_{T}$, our procedure can be applied to any BSM searches with leptons and missing energy in the final states. Scenarios with leptons and missing energy are particularly well-motivated in many dark matter and dark photon searches at the LHC Gramstad:2013loc ; Aad:2012pxa ; Aad:2014vma ; Baek:2013fsa ; Kruker:2013joa ; Wittkowski:2013fua ; Cheng:2013hna ; Teroerde:2016dtv ; Altmannshofer:2016jzy . More generally, we argue that processes with $E_{T}{/}\hskip 9.95863pt$that have a specific topology can benefit from MEM, while not suffering from its computational challenges.

In conclusion, in this paper we provide a working example where, after approximating unknown momenta using the topology of the process, the MEM yields a superior sensitivity compared to the cut-and-count method without having to integrate. Consequently, even if applying the canonical MEM to BSM searches has obstacles and appears to be computationally challenging, we may be able to modify MEM with reasonable assumptions to ease the computational difficulty and yet gain a better sensitivity than the cut-and-count method.

Acknowledgments

We thank Nirmal Raj for his valuable suggestions on the draft and Joe Bramante, Rodolfo Capdevilla, Carlos Alvarado, and Antonio Delgado for useful discussions. This work was partially supported by the National Science Foundation under Grants No. PHY-1417118 and No. PHY-1520966.

Appendix A The Modified Squared Matrix Element of $\left.\tau^{+}\tau^{-}\right|_{\rm dimuon}$ Process ($|\mathcal{M}^{mod}_{\tau\tau}|^{2}$)

In this appendix, we detail the calculation of $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$. The first step is to define the four-vector of $\nu s$ (the vector sum of the two neutrinos coming from each tau) using the momenta of muons in the framework of assumptions discussed in Section IV.2, keeping in mind the $\nu s$ are not massless:

[TABLE]

where the subscript $i$ is defined such that $\tau_{i}\to\mu_{i}\nu s_{i}$. In the following, we will define $p_{T}^{\nu s_{i}},\,m_{\nu s_{i}},\,\eta_{\nu s_{i}},$ and $\phi_{\nu s_{i}}$ in terms of known or measurable parameters:

[TABLE]

where we have ignored the muon mass $(m_{\mu}=0)$, and $\epsilon_{\phi}$s are calculated from the conservation of momenta in the transverse plane.

Furthermore, one of the assumptions (III) in Section IV.2 is that the taus are produced from an on-shell $Z$. Therefore, in calculating $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$, we will also assume the process is $Z\to\tau^{+}\tau^{-}\to\mu_{1}\mu_{2}\,\nu s_{1}\nu s_{2}$, where $Z$ is simply the vector sum of the final state products, shown by $p^{Z}$. The modified squared matrix element of $\left.\tau^{+}\tau^{-}\right|_{\rm dimuon}$ is

[TABLE]

The numerical coefficient in front of $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$ is irrelevant and thus can be ignored. This is because it does not matter whether the plots are with respect to $|\mathcal{M}^{mod}_{\tau\tau}|^{2}$ or $16\pi|\mathcal{M}^{mod}_{\tau\tau}|^{2}$.