Interference effects in dilepton resonance searches for Z′ bosons and dark matter mediators

New Z′ gauge bosons arise in many extensions of the Standard Model and predict resonances in the dilepton invariant mass spectrum. Searches for such resonances therefore provide important constraints on many models of new physics, but the resulting bounds are often calculated without interference effects. In this work we show that the effect of interference is significant and cannot be neglected whenever the Z′ width is large (for example because of an invisible contribution). To illustrate this point, we implement and validate the most recent 139 fb−1 dilepton search from ATLAS and obtain exclusion limits on general Z′ models as well as on simplified dark matter models with spin-1 mediators. We find that interference can substantially strengthen the bound on the Z′ couplings and push exclusion limits for dark matter simplified models to higher values of the Z′ mass. Together with this study we release the open-source code ZPEED, which provides fast likelihoods and exclusion bounds for general Z′ models.


Introduction
The dijet and dilepton final states are amongst the simplest channels currently considered by the LHC collaborations. While dijet resonance searches [1,2] have the advantage that any new particle produced from qq annihilation necessarily can decay back into a pair of quarks, searching for an excess can be difficult due to the large QCD background. Instead, dilepton searches look for a similar bump-like feature above a much smaller electroweak background and achieve great sensitivity to any new particle that couples to Standard Model (SM) leptons [3][4][5][6].
In particular, these searches for exotic resonances offer us a powerful way to probe theories with a new spin-one mediator Z . Such Z bosons generically appear in many extensions of the SM [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] and are an essential part of Grand Unified Theories [22,23]. On a more phenomenological level, they have also received substantial attention as the mediator of spin-one simplified models of Dark Matter (DM) . These simplified models have been advertised by the LHC DM working group [45,46] in order to explore the complementarity between different LHC analyses and across different DM experiments, which include direct and indirect detection, as well as observations of the DM relic density.
While originally these simplified models focused exclusively on the interactions between DM and quarks, it was soon pointed out that lepton couplings cannot be neglected.

JHEP03(2020)104
In models where the Z couples differently to left-and right-handed quarks, the presence of lepton couplings is imposed both by considerations of gauge invariance and by the requirement that there are no gauge anomalies (assuming no exotic SU(2) fermions or additional Higgs doublets). But even in models with vector-like couplings to quarks, lepton couplings generally arise through loop-induced kinetic mixing. Searches for dilepton resonances therefore often place the strongest constraints on simplified DM models and in many cases exclude the most interesting regions of parameter space [47][48][49][50][51][52][53].
In this work, we point out that existing bounds on simplified DM models from dilepton resonance searches are inaccurate, because they neglect the effect of interference between the Z signal and the SM Drell-Yan background pp → Z * /γ * → + − [9,13] (see also ref. [54]). It is commonly assumed that the impact of this interference is negligible, which is typically a good approximation for narrow resonances. However, in the context of simplified DM models, this assumption is not justified because decays of the Z into DM particles can give a large additional contribution to the width of the Z , called the invisible width. If the size of the DM coupling is larger than the SM couplings, the width of the Z will significantly increase as the phase space for the invisible decay opens up.
We demonstrate that the effect of interference can be large, in particular if the signal is smaller than the background and spread out across several bins. In particular for small Z masses (m Z < 2 TeV) and large widths (Γ Z /m Z > 3 %), upper bounds on the couplings can improve by up to a factor of 1.5. The code used to obtain these results is publicly available and can be downloaded from https://github.com/kahlhoefer/ZPEED. In the interest of computational speed, the code makes essentially no use of Monte Carlo event generators, relying instead on analytical cross section calculations and exploiting that the effects of parton distribution functions (PDFs) and analysis cuts are essentially modelindependent.
Our paper structure is then as follows. Section 2 introduces our calculation of cross sections for dilepton processes and shows the effect of interference on signal shapes. In section 3 we describe our implementation of an ATLAS search for dilepton resonances with 139 fb −1 of data [5], including the modeling of detector effects, the statistical method, and a validation via comparison to published bounds. In section 4 we then present our results on the importance of interference effects for various Z models, with a special focus on a DM simplified model, with benchmark couplings recently proposed by the LHC DM working group [46]. Finally, our conclusions are presented in section 5.

Interference effects for vector resonances
In this section we describe the calculation of the cross section for pp → Z → + − (see figure 1) at leading order, including the effect of interference with the SM background processes mediated by the Z boson and the photon. This issue has previously been studied in the context of specific Z models in ref. [9]. For this purpose we introduce a generic Z model with the interaction Lagrangian where Z is the spin-one mediator (with mass m Z ), f is a SM fermion and g V /A are vectorial/axial couplings. Since we wish to remain agnostic about the possible existence of additional contributions to the total width, we treat the decay width of our Z , denoted by Γ Z , as a free parameter in this section. The cross section for the full hadronic process can be related to the partonic one for the hard process as where the sum is performed over all quark and anti-quark flavours. Here, the x i denote the momentum fractions of the individual partons and f q and fq are the MSTW PDFs [56], which we evaluate setting the factorisation scale to µ = m . It is straight-forward from this expression to calculate the differential cross section with respect to the dilepton invariant mass dσ/dm (see appendix A). This cross section can be split into several parts: In dilepton resonance searches the SM background is typically large (at least for m 2 TeV) but known with a high level of precision. These searches are therefore potentially sensitive to exotic resonances even if in any given bin σ signal σ background . 1 For many Z models the width Γ Z is small compared to the bin size. In this case the signal will only be observable if dσ signal /dm dσ background /dm for m ≈ m Z . Since it follows that dσ signal /dm dσ interference /dm , so that interference effects are typically not important. If on the other hand Γ Z is comparable to the bin size (for example because of an invisible decay mode), dilepton resonance searches are potentially sensitive to signals with dσ signal /dm dσ background /dm for all values of m . For such small signal cross sections, interference effects can potentially be very important. 1 Here we define σ = b a (dσ/dm )dm for a bin given by m ∈ [a, b]. Typical bin sizes are comparable to the detector resolution, which is approximately 1-2% in the electron channel and 5-10% in the muon channel. Note that for the purpose of this section we neglect detector effects, which will be discussed in detail in section 3.1. We illustrate the effect of interference in figure 2 for a narrow Z signal with width 2.5 GeV (left panel) and a broad signal with width 15 GeV (right panel), keeping the couplings and resonance mass fixed. 2 The differential cross section as a function of m is shown for the naive sum of signal and background without interference and with interference included.
As expected, we find that interference becomes more important for larger widths, because of the suppression of the pure signal term compared to the interference term. In detail the effect of interference depends on the sign of the Z couplings. For g V q g V > 0 interference is constructive for m < m Z and destructive for m > m Z . Since the background is monotonically falling, this leads to an increase in the height of the peak and a shift of its location to smaller values of m . For the opposite case (g V q g V < 0) the height of the peak still increases, but the peak is now shifted to larger values of m . In the following, we will focus on the case that g V q g V > 0. Results for the opposite case are summarized in appendix B.
To conclude this section, we note that interference effects are more relevant for Z mediators with vector couplings than for those with axial couplings. The reason is that axial mediators do not interfere with the photon, which gives the dominant contribution to interference for vector mediators. We will therefore restrict ourselves to vector mediators in the following.

Analysis set-up
In this section we describe how to translate the theoretical cross section from above into realistic predictions of the expected number of events in a given set of bins of the dilepton JHEP03(2020)104 invariant mass m , including analysis cuts, detector efficiencies, energy resolution and higher-order effects. We then give a brief summary of the statistical method that we employ to test whether or not the resulting signal prediction is compatible with data at a given confidence level. Finally, we perform a validation of our analysis set-up by reproducing the published bounds on the production cross section of Z bosons with given width from the ATLAS collaboration [5].

Signal prediction
For a given bin i covering some range of m , the prediction for the number of detected electron or muon pairs ( = e, µ) is written as where L is the luminosity, dσ dm is the differential signal cross section including interference, W i (m ) denotes a window function reflecting the finite detector resolution, and ξ (m ) is a rescaling factor taking into account higher-order corrections and detector efficiencies. The different ingredients of the predictions will be explained in the following.
We perform a fully differential leading-order (LO) computation for the Drell-Yan cross section including a Z mediator. The computation is implemented in a fast and efficient computer code as further detailed in section 3.4 and appendix A. The fiducial phase-space volume of the ATLAS analysis is defined by p T > 30 GeV for electrons as well as muons. Concerning rapidity, we accept electrons with |η| < 1.37 or 1.52 < |η| < 2.47 and muons with |η| < 2.5. Integrating over the fiducial volume for fixed invariant dilepton mass m , we obtain dσ dm . We also calculate the SM Drell-Yan background dσ SM dm , i.e. the first three terms in eq. (2.3), in complete analogy to the signal.
The limited detector resolution is reflected in our analysis using a simple Gaussian kernel which smears the calculated invariant mass spectrum. For a bin defined by m ∈ [a i , b i ] the Gaussian smearing is implemented using the window function where the detector resolution s(m ) is taken from the auxiliary figures of the ATLAS analysis in ref. [5]. Unfortunately, detector efficiencies cannot be included at the fully differential level since we lack the full experimental information. In particular, quality requirements for the muon or electron identification cannot be approximated by a simple detector simulation like DELPHES [57]. 3 However, we can make use of the published predictions for the SM Drell-Yan background in order to estimate detector efficiencies as a function of m and then improve our LO prediction dσ SM dm by appropriate rescaling factors ξ (m ). In addition, the rescaling also approximately captures higher-order corrections beyond LO in perturbation theory as discussed at the end of the section.  The rescaling factors ξ (m ) are derived as follows. Tables 3 and 4 in ref. [3] list the expected event yields s exp ,i for the Drell-Yan background in wide bins of m . 4 We calculate the corresponding event yields based on our LO calculation including detector resolution. We then define ξ The rescaling factors obtained in this way are stated in table 1. The function ξ (m ) is then obtained by linear interpolation. Since the resulting functions ξ (m ) depend only weakly on m , the simple linear interpolation turns out to be a sufficient approximation.
As an alternative approach, we have first used DELPHES on a fully differential level to include those detector effects that are implemented. Additional detector effects not included in DELPHES are then again included by our rescaling approach. The differences between the two approaches are negligible. Hence, for simplicity, we do not use any detector simulation by DELPHES for the results presented in the following.
As noted above, our LO cross section is not only modified by detector effects but also by higher-order corrections. The dominant higher-order corrections are approximately included in our rescaling procedure as well because they are included in the expected event yields s exp ,i . Like the detector efficiency, the corrections are not included at the fully differential level but they are effectively treated as m -dependent K-factors along with the detector effects. Here, we assume that the higher-order corrections affect the SM background in the same way as the differential signal cross section including interference. This is certainly true for the QCD corrections to the Drell-Yan process, which only concern the initial state.

Statistical method
Having calculated the predicted signal s i in each bin, we can construct the likelihood may depend on additional nuisance parameters, in which case −2 log L(µ) denotes the profile likelihood (where all nuisance parameters have been set to the values that maximise the likelihood for given µ). The contribution from interference between signal and background is included in the predicted signal s i . Since signal and interference depend differently on the parameters of the underlying model, the term µs i is unphysical for general values of µ in the sense that it does not correspond to any parameter combination. Nevertheless, introducing µ is a useful construction to interpolate between the signal+background hypothesis (µ = 1) and the background-only hypothesis (µ = 0) without changing the shape of the signal. The value of µ that maximises the likelihood is calledμ. Having found this value, we calculate the test statistic which is expected to follow a χ 2 distribution with 1 degree of freedom. Rather than calculating exclusion bounds directly from q µ , we employ the CL s method [58]. In the asymptotic regime (b i , o i 1), the modified p-value of the signal+background hypothesis is given by 5 Here Φ is the cumulative distribution function of the normal distribution and q A,µ is the value of the test statistic q µ for the Asimov data set [59], in which all observations exactly match the background expectation The signal+background hypothesis (µ = 1) can now be rejected with (at least) 95 % confidence level if CL s ≤ 0.05. It is common practice to solve CL s = 0.05 for µ in order to find the smallest value of µ that is excluded. However, as discussed above only µ = 0 and µ = 1 represent actual physical models. In the following, we will therefore not quote bounds on µ but instead apply the CL s method to every point in parameter space in order to identify those parameter regions where µ = 1 is excluded.
At present only ATLAS provides publicly available data for dilepton resonance searches based on the entire data from Run 2 [5], and we focus on their analysis here. 6 In contrast to previous dilepton resonance searches, ATLAS does not rely on Monte Carlo simulations to estimate backgrounds, but instead obtains the background estimates by fitting a smooth function to the observed data. In principle, the uncertainties on the fit parameters obtained in this way should be included as nuisance parameters. However, given that the background is fitted across many different bins, while the signal is more localised, the uncertainties in the nuisance parameters have a negligible impact on the profile likelihood. For our implementation we therefore simply take b i to be the central value of the background prediction. 5 For large values of m Z the assumption of asymptotics leads to exclusion limits that are too strong by a factor of 2 or more. The main focus of the present work is however on m Z 2 TeV, where the asymptotic expression for CLs provides a very good approximation. 6 We have checked that including the publicly available data from CMS based on an integrated luminosity of 36 fb −1 [4] does not substantially change any of the results that we present. To reproduce the ATLAS analysis as closely as possible, we exclude the contribution from off-shell Z bosons at small m . Specifically, we limit ourselves to the signal region is the detector resolution in the e + e − channel. In general m ,min will not coincide with the boundary of any bin. The bin [a, b] that satisfies a < m ,min < b is included in the likelihood, but its contribution is multiplied with the weighting factor where s min is the number of signal events in the interval [m ,min , b]. This approach ensures that the likelihood is a continuous function of m Z . We impose no upper bound on m other than the requirement m < 6253 GeV implied by the ATLAS data.

Validation
In order to check our implementation of the detector efficiencies, smearing and rescaling functions, as well as our statistical analysis, in this subsection we validate our results by comparing to bounds published by the ATLAS collaboration in ref. [5]. We show in figure 3 our bound on the cross section as a function of the Z mass, compared to that published by ATLAS [5]. Note that in order to reproduce the approach taken by the experimental analysis, these bounds are calculated ignoring the effect of interference with SM Drell-Yan processes. The bounds are for 95 % C.L. and show good agreement for a Z width of 0.5 % (left) and 3 % (right). We have also checked the bounds for larger Z widths and find good agreement up to Γ Z /m Z ≈ 6%. For even larger widths correlated background uncertainties, which cannot be properly included with publicly available information, become important and our approach yields bounds that are slightly stronger than the published ones. For signal widths smaller than 0.5 %, on the other hand, bounds will be dominated by detector resolution and will be very similar to the case shown in the left JHEP03(2020)104 panel of figure 3. We hence conclude that our implementation is reliable for any Z signals with Γ Z /m Z 0.06.

ZPEED
To obtain these results, we have developed a highly efficient numerical code called ZPEED (Z Exclusions from Experimental Data), which is capable of calculating the likelihood and CL s value for a given Z parameter point within less than a second on a single CPU. The code implements the approach outlined in appendix A, i.e. it uses analytical expressions for the differential cross sections of signal and interference terms together with tabulated values of the function T q,2 (m ) as defined in eq. (A.5), which accounts for PDFs and phase space cuts. The differential cross sections are then multiplied with the rescaling factors ξ (m ) and the window functions W i (m ) introduced in eq. (3.1). Indeed, the integration over m in eq. (3.1), which needs to be performed at runtime, is the only computationally expensive step. Once the predictions s i have been calculated, it is straightforward to calculate the likelihood defined in eq. (3.4) as a function of the signal strength µ, determineμ and obtain the CL s value. At present only the ATLAS analysis based on 139 fb −1 has been implemented, but future updates will be provided whenever new data becomes publicly available. The code is open source and can be downloaded from https://github.com/kahlhoefer/ZPEED.

Results
In this section we illustrate the importance of interference effects by showing how they impact bounds derived from experimental data. We will first do this in a model-independent way by treating couplings and width as independent parameters and then focus on a specific simplified model, in which the width of the Z is calculated self-consistently as a function of the underlying parameters.

Model-independent bounds
We first consider a general Z model with vector couplings and define the effective coupling g ≡ (g V q g V ) 1/2 . For fixed total width Γ Z the Z production cross section is proportional to g 4 , while interference effects scale as g 2 . We can therefore use the analysis chain presented in section 3 to calculate bounds on g with and without interference for different values of Γ Z .
The resulting exclusion bounds are shown in figure 4 as a function of the mediator mass m Z for Γ Z /m Z = 0.01, 0.03 and 0.06. As expected, interference effects are negligible when the relative width is small (top panel) and become increasingly important as Γ Z /m Z increases. Interference effects are largest for small values of m Z , which is a consequence of the steeply falling SM background. In the bottom panel, which assumes a 6 % relative width, interference effects lead to a strong enough distortion of the input signal such that the exclusion limits are changed significantly. For instance, for m Z ≈ 500 GeV, the exclusion limit on g obtained from the pure Z signal is about a factor of 1.  . Upper bound on the effective coupling g = (g V q g V ) 1/2 at 95 % confidence level, with and without interference effects. We consider Z bosons with vanishing axial couplings and different relative widths Γ Z /m Z = 1%, 3%, 6%.

JHEP03(2020)104
in the dilepton invariant mass spectrum to smaller values (see figure 2), which results in a shift of the exclusion bound to larger masses. For example, the dip around m Z ≈ 1.9 TeV in the bottom panel is shifted to about m Z ≈ 2 TeV once interference effects are included.
Instead of assuming equal couplings to electrons and muons, one can also calculate constraints on the effective coupling g to each lepton family separately. The resulting upper bounds are provided in appendix B.
We emphasize that for large relative widths the impact of interference effects is at least as important as the impact of higher-order QCD corrections. In particular, the former can significantly change the shape of the signal, while the latter only result in an effective rescaling of the cross section that can be applied after signal events have been generated. Interference effects, on the other hand, need to be included during signal generation and depend in a more complicated way on the underlying parameters. It is essential to include these effects in order to obtain accurate bounds on the parameter space of a given Z model. In most cases including interference effects leads to stronger exclusion limits, which further enhances the potential of dilepton resonance searches to constrain models of BSM physics.

Bounds on dark matter simplified models
As we have seen above, interference effects are most important for large relative widths. Such large widths typically cannot be obtained from decays into SM particles (as the required couplings would violate experimental constraints), but they are a generic prediction in models with additional contributions to the Z width arising from decays into new invisible light degrees of freedom. As a specific example of such a model, we consider a spin-one simplified DM model [25], which has been employed by the LHC collaborations [45,46] to create benchmark points in theory space that allow for different LHC DM searches to be compared to each other and to non-collider experiments.
We extend eq. (2.1) to include a coupling to a SM singlet Dirac fermion χ with mass m χ as a DM candidate. The corresponding interaction Lagrangian reads Then each partial width of the Z is where N c is the number of colours. It has been shown that (for a minimal Higgs sector) gauge invariance requires g A = g A q [47], which typically leads to overwhelmingly strong constraints from dilepton resonance searches in models with non-zero axial couplings. We therefore focus on the case g A q/ /χ = 0, while the three remaining couplings g V q , g V and g V χ are treated as independent parameters. A particularly well-motivated possibility is that g V vanishes at high scales and is only introduced at low scales through kinetic mixing [48]. In this case one naturally finds g V q g V > 0, such that bounds from dilepton resonance searches are suppressed but still relevant. In the simplified DM model introduced above, Γ Z depends decisively on the mass hierarchy. For m χ > m Z /2, invisible decays are kinematically forbidden and the relative width is very small. In the opposite case, the partial width Γ (Z −→ χχ) may contribute significantly to the total width, in particular if g V χ g V q . Following the recommendations of the LHC DM working group [46], we therefore consider the benchmark choice g V χ = 1.0 and g V q = 0.1, such that Γ Z /m Z ≈ 0.5 % for m χ > m Z /2 and Γ Z /m Z ≈ 3 % for m χ m Z /2. We consider the two choices g V = 0.01, 0.02, corresponding to an effective coupling g = (g V q g V ) 1/2 = 0.032 and g = 0.045, respectively. Figure 5 displays the resulting exclusion limits in the m Z -m χ -plane both with and without the inclusion of interference effects. 7 We emphasize again that in these plots the 7 We note that the exclusion limit obtained in the absence of interference effects is slightly stronger than JHEP03(2020)104 decay width is computed following eq. (4.2). As expected, interference effects are most important for m χ < m Z /2, corresponding to larger relative width of the Z , and for small g V . When interference effects are neglected, the parameter region with small m χ is essentially unconstrained for m Z 850 GeV (m Z 1650 GeV) in the case that g V = 0.01 (g V = 0.02). Including interference effects, the parameter region probed by dilepton resonance searches is extended to m Z 1200 GeV (m Z 2050 GeV). Although the precise parameter regions excluded by the ATLAS analysis depend sensitively on fluctuations in the data, the general trend is clear: interference effects lead to stronger bounds on the simplified DM model.
To conclude this discussion, we note that χ should, as a DM candidate, also satisfy bounds coming from the relic density of DM and from direct and indirect detection experiments in addition to collider bounds. A number of works have investigated in detail the complementarity of these different constraints (see e.g. refs. [24-35, 37-42, 44-49, 51-53]). Here we focus on the contribution of χ to the total decay width of the Z and the resulting interference effects. Therefore, we do not make any assumptions on the cosmological history and the relic abundance of χ. In fact, all of the results presented in this work remain valid even if χ is unstable and decays into either SM particles or other BSM states.

Conclusions
In this work we have investigated the sensitivity of the LHC to new Z bosons with a focus on the effect of interference between the Z signal and SM Drell-Yan background in the dilepton channel. Interference is enhanced for Z bosons with large width (compared to the detector resolution), arising for example from invisible decay modes into new light degrees of freedom, and results in an asymmetric signal with modified peak amplitude and position (see figure 2). Details of our calculations and of the fast numerical implementation can be found in appendix A.
In order to quantify the impact of interference on bounds derived from experimental data, we have implemented an existing ATLAS search for dilepton resonances. We use smearing functions to model energy resolution based on experimental data and estimate experimental efficiencies and higher-order corrections by rescaling our predicted Drell-Yan background to published background estimates. We have calculated exclusion bounds on the fiducial cross section neglecting interference with the CL s method and found excellent agreement with published limits (see figure 3). We have made the code used to obtain these results publicly available. 8 We then applied this analysis to the case of a Z with purely vectorial couplings in order to obtain bounds on the effective coupling g = (g q g ) 1/2 as a function of m Z for different values of Γ Z (see figure 4). As expected, interference effects are most important for large widths and can substantially strengthen the bounds on the effective coupling g. For example, for a Z with 6 % relative width the bound on the couplings improves by up to a factor of 1.5 once interference is included. the one provided by the ATLAS collaboration. In the absence of a detailed documentation it is difficult to identify the origin of this discrepancy. 8 ZPEED -Z Exclusions from Experimental Data: https://github.com/kahlhoefer/ZPEED.

JHEP03(2020)104
We also considered a specific example for a model where the Z width can be large in spite of small couplings to quarks and leptons, namely a simplified model of DM with a spin-1 mediator. Assuming the DM coupling g χ is large compared to g q and g , decays of the mediator into DM particles give rise to a large invisible width and therefore a substantial increase of the total width, whenever decays into DM are kinematically allowed. In this model the Z width can easily be large compared to the detector resolution and therefore large enough for interference effects to be relevant. As a specific benchmark we considered g χ = 1, g q = 0.1, and g = 0.01 (a choice recommended by the LHC Dark Matter Working Group and used by both ATLAS and CMS to present exclusion limits [46]), as well as an additional example with g = 0.02. Deriving bounds on this model as a function of DM mass and mediator mass, we demonstrated that interference effects lead to substantially stronger constraints on the parameter space of this model (see figure 5).
The LHC is entering the phase where precise signal predictions are essential in order to fully exploit the benefits of high statistics. We argue that in order to accurately calculate constraints on Z bosons with large widths from dilepton resonances, interference effects must be included. This is particularly true in DM models where the Z acts as the mediator between DM and SM fermions and obtains a large invisible width. We encourage the experimental collaborations to incorporate the modified signal shapes and look forward to the inclusion of interference effects in bounds coming from existing and new LHC data.

A Cross section calculations
In this appendix, we provide more details on the calculation of the partonic and hadronic Drell-Yan cross section in the Z model under consideration. The interaction Lagrangian of the Z has been introduced in eq. (2.1). The differential LO result for the partonic signal cross sectionσ Z Z reads whereŝ andt are the usual Mandelstam variables, N c = 3 for QCD, m Z is the mass of the Z and Γ Z its total width. The coupling coefficients read (A.2)

JHEP03(2020)104
Convolving the partonic cross section with parton-distribution functions of the quarks f q and anti-quarks fq, the fully differential hadronic cross section is given by where m = √ŝ is the dilepton invariant mass, η ± are the rapidities of the positively and negatively charged leptons in the lab frame, x i are the momentum fractions of the partons, and y = 1 2 (η + − η − ). The sum runs over all light quark and anti-quark flavours. To obtain this result, we have made use of the following relations between the different kinematic variables:t = − m 2 2 cosh y e −y , where Y = 1 2 (η + + η − ). Hence, we can define where i = 0, 1, 2 and it is understood that we only integrate over the fiducial region, i.e. the cuts on the rapidities and the lepton transverse momenta p T = m /(2 cosh y) are respected. With this definition, we write the differential cross section as a product of model-independent but cut-dependent function T q,i (m ) and simple modeldependent factors consisting of couplings and propagators. Employing MSTW parton distribution functions [56], the T q,i -functions can be evaluated once on a fine discrete m -grid and linearly interpolated, such that no numerical integrations have to be performed when the cross section is evaluated for different model parameters. Hence, eq. (A.6) is a particularly efficient implementation for parameter scans. Note that this separation only works for s-channel mediated interactions like the Drell-Yan like process under consideration, since only in this case the propagator does not depend on the rapidities. The evaluation of the hadronic cross section further simplifies, since for our (symmetric) fiducial volume one has T q,0 + 2T q,1 = 0 (A.7) and c q 1 = 2c q 0 implies that T q,0 and T q,1 do not contribute to the cross section. So far, we have only discussed the signal cross section dσ Z Z /dm without interference. However, all considerations apply with trivial modifications to the interference terms and the SM Drell-Yan background as well. Hence, as our final result, the cross section in eq. (2.3) can be calculated from is given in terms of the couplings of the vector bosons to fermions defined for the photon and the Z boson in analogy to eq. (2.1).

B Further exclusion limits
In figure 6 we show separate constraints on the effective coupling g to electrons and muons. We find constraints on the former to be slightly stronger than on the latter, which is a JHEP03(2020)104 direct consequence of the better detector resolution for the electron final state. Although fluctuations are more pronounced in the electron channel, the effect of interference is similar in both cases.
In figures 7, 8 and 9 we present our results for the case that the product of quark and lepton coupling are negative (g V q g V < 0). As can be seen from figure 2, this changes the shape of the expected signal substantially. Crucially, interference still leads to an increased height of the peak and therefore including interference effects typically leads to stronger exclusion bounds.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.  Figure 9. Same as figure 6 but for the case that g V q g V < 0.