Topping-up multilepton plus b-jets anomalies at the LHC with a $Z'$ boson

During the last years ATLAS and CMS have reported a number of slight to mild discrepancies in signatures of multileptons plus $b$-jets in analyses such as $t\bar t H$, $t\bar t W^\pm$, $t\bar t Z$ and $t\bar t t\bar t$. Among them, a recent ATLAS result on $t\bar t H$ production has also reported an excess in the charge asymmetry in the same-sign dilepton channel with two or more $b$-tagged jets. Motivated by these tantalizing discrepancies, we study a phenomenological New Physics model consisting of a $Z'$ boson that couples to up-type quarks via right-handed currents: $t_R\gamma^\mu \bar t_R$, $t_R\gamma^\mu \bar c_R$, and $t_R \gamma^\mu \bar u_R$. The latter vertex allows to translate the charge asymmetry at the LHC initial state protons to a final state with top quarks which, decaying to a positive lepton and a $b$-jet, provides a crucial contribution to some of the observed discrepancies. Through an analysis at a detector level, we select the region in parameter space of our model that best reproduces the data in the aforementioned $t\bar t H$ study, and in a recent ATLAS $t\bar t t \bar t$ search. We find that our model provides a better fit to the experimental data than the Standard Model for a New Physics scale of approximately $\sim$500 GeV, and with a hierarchical coupling of the $Z'$ boson that favours the top quark and the presence of FCNC currents. In order to estimate the LHC sensitivity to this signal, we design a broadband search featuring many kinematic regions with different signal-to-background ratio, and perform a global analysis. We also define signal-enhanced regions and study observables that could further distinguish signal from background. We find that the region in parameter space of our model that best fits the analysed data could be probed with a significance exceeding 3 standard deviations with just the full Run-2 dataset.


Introduction
Since the beginning of operation of the Large Hadron Collider (LHC), the ATLAS and CMS collaborations have carried out a broad program of precision measurements of Standard Model (SM) parameters and processes, as well as sensitive searches for new phenomena beyond the SM (BSM), all of which have furthered our understanding of the fundamental interactions of Nature. However, we know the SM to be an incomplete theory and we still hope to find New Physics at the LHC. Being a proton-proton collider, the LHC is a discovery machine with an energy range that would ideally produce BSM resonances at the TeV scale. These resonances have been searched for in many different channels, motivated by a variety of BSM models. As these resonance searches have excluded a larger portion of the parameter space for the simplest models, the attention has turned to precision measurements where deviations from SM predictions could provide hints to BSM signatures. Such BSM could still be produced in a resonant way but, without specific dedicated searches, this resonant feature could go unnoticed.
Among the most recent SM benchmarks being explored, the LHC program has made impressive progress in the measurement of the cross-sections of ttW ± [1,2], ttH [3,4] and four-top-quarks (tttt) [5,6] production. To measure these cross-sections, a careful choice of the final state must be made. The same-sign dilepton and multilepton final states with a high b-jet multiplicity present a good balance between low SM backgrounds and high signal yield. In these final states, a persistent discrepancy has been present in the ttW ± normalization, as the observed yields are consistently larger than those expected from state-of-art theoretical predictions. This discrepancy has been found both in dedicated ttW ± cross-section measurements [1,2] and in several searches that consider the ttW ± process as a background [3][4][5][6], as summarised in Table 1. We see that the ATLAS and CMS ttH measurements [3,4] analyse a larger dataset than the ttW ± cross-section measurements [1,2]. If we consider the same reference cross-section for all searches, using for example the reported cross-section of Ref. [7], we see that the ttH measurements show a larger tension than the ttW ± measurements. Because of this two facts, we focus on the ttH searches. Despite having similar strategies to suppress non-prompt lepton background, there are several differences in the ttH analysis approach between ATLAS and CMS, such as different fake estimation techniques (simultaneous profile likelihood template fit in ATLAS vs misidentification probability method in CMS) and different analysis strategies (jet multiplicity, total lepton charge, and Boosted Decision Trees (BDT) categorisation in ATLAS vs Deep Neural Networks in CMS). A summary of the event selection used in the ATLAS and CMS ttH analyses can be found in Table 2. We focus on the two same-sign dilepton and trilepton channels where hadronically decaying τ s are vetoed.
We choose the ATLAS ttH measurement for the re-interpretation in this work since it provides complete information on the multilepton discrepancies observed as a function of the total lepton charge and the b-jet multiplicity of the event. A recent CMS search for new physics within the effective field theory (EFT) framework [8] using 41.5 fb −1 of integrated luminosity exploits a similar categorisation. However, since there is no information on the two same-sign dilepton selection with exactly 1 b-jet and given the lower data statistics, we do not include the CMS EFT results in the current work.
In addition, a four-top-quarks search by the ATLAS Collaboration in the same final states has measured a four-top-quarks cross-section a factor of two higher than the SM prediction [5]. A summary of the signal region event selection used in the ATLAS and CMS four-top-quarks analyses can be found in Table 3.
It should be pointed out that these measurements are individually consistent with the SM and thus the observed discrepancies could be merely the result of statistical fluctuations, and/or unaccounted experimental or theoretical uncertainties. However, when combined, they paint an interesting picture that is worth exploring, as it may open the door to new exciting discoveries.  Table 1: This table lists the ttW ± reference cross-sections and the corresponding signal strengths for different searches. The last column is the signal strength corresponding to the reference cross-section listed in Ref. [7], 600.8 fb. The reference cross-section was not listed in Ref. [6] and so it was taken from Ref. [9].  > 500 > 300 |m e ± e ± | (2LSS) or > 15 -> 12 |m OSSF | (3L) [GeV] |m e ± e ± − m Z | (2LSS) or > 10 -> 15 |m OSSF −m Z | (3L) [GeV] Other -Missing transverse momentum cuts Table 3: Comparison of event selections between the ATLAS [5] and CMS [6] four-top-quarks analyses. H T is the scalar p T sum of jets, leptons and b-jets.
In this work, we propose a purely phenomenological BSM model, where a Z spin-1 boson with mass below 1 TeV is responsible for the reported discrepancies in ttH and four-topquarks searches. We choose this model with a particular coupling structure as it is able to provide both charge asymmetry and high b-jet multiplicity, key ingredients to accommodate the different experimental signatures. This paper is organized as follows. We present our phenomenological model in Sect. 2. We then study its experimental imprints in Sect. 3, including estimating the region of parameter space that best fits the ttH and four-top-quarks ATLAS data, and possible constraints from other observables on this parameter space. In Sect. 4 we propose a global search strategy to discriminate signal from background. Finally, we present our conclusions in Sect. 5.

The model
We aim to construct a model that can account for several slight excesses in multileptons plus b-jets final states at the LHC. More precisely, we require the model to be capable of producing more positive than negative leptons to better fit the imbalance suggested by the results in Ref. [3]. We motivate the model from a phenomenological point of view.
Since the imbalance in the final states reported in Ref. [3] has more positive than negative leptonic charge, the model needs to capture the excess in positive charge present in the colliding protons. Thus the BSM should couple to up quarks, while being safe to low-energy physics observables. If the new particle coupling to up quarks would be charged, then it would also couple to bottom quarks (a W boson) or to leptons (a Leptoquark). If the new particle were a W , then it would be difficult to produce an excess in positive multileptons with b-quarks in the final state. Were the new particle a Leptoquark then, being charged under SU (3) C , the bounds on its mass from pair production would make it more difficult to reproduce the observed excesses. We are then left with neutral particles coupling to up quarks (u, c and t) and with non-negligible Flavour Changing Neutral Currents (FCNC) to account for the deviations. Among the three usual spins 0, 1 or 2, we find suitable to study spin-1 since, if color neutral, its gluon fusion production is protected through the Landau-Yang theorem [10][11][12], which also guarantees that the restrictive bounds in di-photon do not apply [13][14][15][16][17]. Examples of studies of a FCNC spin 0 scalar boson phenomenology at the LHC, can be seen in Refs. [18,19].
A color neutral spin-1 particle with non-diagonal couplings is known as a FCNC Z . Stringent constraints from the LEP II and Drell-Yan experiments require tiny or null coupling to leptons, thus for the purposes of this work we restrict to a leptophobic Z . Moreover, in light of limits set by low energy physics experiments and precision tests [20], we restrict Z to only couple to right quark flavour changing ut and ct and flavour conserving tt currents, In addition, one should consider a kinetic Lagrangian and an eventually negligible Z-Z mixing because of the restrictive bounds imposed by LEP [21]. For the sake of simplifying notation, in the following of this article -except in Appendix B-, we drop the R supraindex from the couplings.
In Section 3.2 we study bounds to the parameter space of this BSM Lagrangian coming from low energy physics and collider phenomenology. Effective theories similar to the one described by the interaction in Eq. 1 have been studied in different contexts, as for instance in Refs. [22][23][24][25][26]. In particular, a very similar set-up has been implemented in Ref. [27], albeit in a different mass range and assuming a hidden sector that increases the width-to-mass ratio, and in Ref. [28] where a similar model was implemented to study the forward-backward tt asymmetry at the Tevatron.

Phenomenology
To account for the observed data (see Sect. 3.1 for more details), we need to produce same electric charge dilepton (denoted 2LSS, with SS standing for same-sign) and multilepton (at least three leptons, denoted 3L) final states with charge asymmetry and high b-jet multiplicity. To accomplish this, we consider two relevant Z -induced processes, tZ + tZ (denoted tZ in the following) and ttZ , with a hierarchy between the relevant couplings to enforce a high probability of three-and four-top-quarks final states. We show two examples of the Feynman diagrams in Fig. 1, and the relevant cross-sections and branching ratios in Fig. 2 as a function of M Z for a given set of couplings. The cross-sections for a different set of couplings can be obtained by simple re-scaling. The branching ratios for a different set of couplings can be obtained by using Γ(Z → tu)/(g ut ) 2 = Γ(Z → tc)/(g ct ) 2 and re-scaling.  Other processes that could be tested in these channels (see Sect. 3.2) are either numerically irrelevant such as same-sign top-quark pair production (tt, t t) [27], or chirality suppressed such as radiative Z production (tZ j, tZ j) [27,29]. Note that these results can be fairly specific to our model, see e.g. Ref. [30] for a different Z model where cg → tZ is the main production channel. In particular, non-resonant effects are suppressed both by the small Z width we consider by neglecting hidden sector decays, and by our particular choice of coupling structure.
From Fig. 1, we see that contributions to tZ are proportional to (g ut ) 2 and to (g ct ) 2 while ttZ is proportional to (g tt ) 2 . From Fig. 2a, we see that even for g tt g ut , g ct , the cross-section for tZ production is larger than for ttZ , mainly due to the kinematic requirements that must be met to produce each of the final-state particles on-shell. When comparing the four possible tZ processes, the largest cross-section is ug → tZ , as expected from the model motivations. This is due to u-quark abundance in the proton, which ensures that the g ut -induced processes yield a considerable charge asymmetry that is not present in the other production processes. We are interested in this charge asymmetry and how it is reflected in current experimental searches.
If we take the possible relevant decays obtained from the Lagrangian in Eq. 1 into account, and assuming they are all kinematically accessible, we see that tZ can yield the following final states: ttj, t tj, ttt and ttt. We are discarding the ttj final state because of the 2LSS and 3L selection criteria. On the other hand, ttZ can yield tttj, tttj and tttt final states.
If we want events enriched with leptonic charge asymmetry and b-jets, we need tZ to decay mostly to three top quarks while being as charge asymmetric as possible. That is, we need BR(Z → tt) > BR(Z → tj + tj). As we consider relatively low M Z masses, we need g tt g ut to avoid phase-space suppression. If M Z < 2m t , we consider the three-body decay Z → tW − b, tW + b. From Fig. 2 we see that for the benchmarks points we display, when M Z is large enough the tt decay mode dominates. When combined with tZ and ttZ production, these decays produce three-and four-top-quark final states.
After this overview of the basic phenomenology of the proposed model, we turn to studying its effect on the relevant observables, and how these observables determine the region in its parameter space most compatible with the experimental results.

Experimental imprints and model tuning
In this section we study how the Z model detailed in Sect. 2 is probed by different existing experimental results. A Z as in Eq. 1 affects both high-and low-energy observables. In light of recent experimental results from ttH and tttt searches in ATLAS and CMS, we are specially interested in multilepton-plus-b-jets final states. In particular, multilepton final states with non-zero total leptonic charge are highly sensitive to Z . Taking this into account, we detail in Sect. 3.1 how on-shell Z production in association either with a single top quark (tZ or tZ ) or with a top quark pair (ttZ ) could explain the need for ttW ± re-scaling to account for tensions in data in recent ttH results [3], while yielding interesting signatures in four-topquarks. In Sect. 3.2 we test whether the parameter space indicated by multilepton and high b-jet multiplicity is safe to other observables that could be affected by our model. The more relevant observables we consider are D 0 −D 0 -meson mixing, top-quark rare decays through flavour-changing neutral currents (FCNC), top-quark pair production, same-sign top-quark pair production, resonant tj production in tt+jets events, and Z radiative production tZ j. Other processes involving resonant Z production, such as single Z , Z j and single top quark production tj + tj, are absent due to having approximated as null the g uu and g cc couplings.

Fits to experimental data
We study how our phenomenological Z model detailed in Sect. 2 accommodates recent results reported for ttH production while avoiding constraints and yielding potentially interesting results in tttt searches. Regarding ttH, we focus on the ATLAS preliminary results reported in Ref. [3], although results from CMS [4] are consistent with our conclusions. Regarding four-top-quarks we consider the results reported by ATLAS in Ref. [5], having corroborated that the results obtained are compatible with the combination of the ATLAS and CMS [6] results. Both searches target same-sign dilepton and multilepton processes but their channel definitions are not the same and the reconstructed objects, both leptons and jets, have different kinematic cuts and tagging efficiencies, which we take into account in our study.
For the case of ttH, the results reported in Ref. [3] are particularly interesting because of the difficulties reported when dealing with the irreducible ttW ± background. Figure 2 of Ref. [3] highlights the need for a missing charge asymmetric contribution to match the data, which yields a normalization factor for ttW ± larger than one. This is consistent with other measurements, for example Refs. [1,2], and has motivated a push for higher theoretical accuracy in the ttW ± calculations [31][32][33][34][35][36][37][38], which nevertheless have not fully explained the discrepancy between the expected and the observed ttW ± event yields.
If we treat the ttW ± background as well modelled, and thus constrained to have a normalization factor consistent with unity within the uncertainty in its theoretical cross-section, we are faced with charge asymmetric anomalous events with high b-jet multiplicity. Our Z model is designed to accommodate these two features. A recent example of BSM effects in ttW can be found in Ref. [9], where, in contrast to a resonant BSM physics model, they study an Effective Field Theory in the top-quark sector.
As detailed in Sect. 2.1, we study Z production in association either with a single top quark or with a top-quark pair. These processes, tZ and ttZ , along with a considerable Z → tt branching ratio achieved with a suitable coupling hierarchy, provide the necessary same-sign dilepton and multilepton signatures, with both charge asymmetry and high b-jet multiplicity.
As we consider three-top-quarks production and four-top-quarks production, we need to make sure that our results are compatible with four-top-quarks limits. Even if four-top-quarks production will be mostly sensitive to (g tt ) 4 , three-top-quarks production is proportional to (g ut ) 2 (g tt ) 2 and will introduce a charge asymmetry.
To see how well our model can agree with the data, we obtain the event yields expected from tZ and ttZ in the different reported bins in Fig. 2 of Ref. [3], and the expected events in four-top-quarks searches with the selection criteria of Ref. [5] as a function of (g ut , g ct , g tt ) for different values of M Z . To do this, we have implemented the Z model in Eq. 1 using Feynrules [39] and then simulated a fixed set of points through Madgraph 5 [40] for production and decay of tZ and ttZ . The signal events have been generated using a leading-order matrix element and the NN23LO1 PDF set [41], and have been processed through Pythia 8 [42] for the modelling of parton showering and hadronization, as well as through a simulation of the detector response as implemented in Delphes [43]. We use the Monash tune [44] of Pythia 8 with a few changes aimed to reproduce the ttW ± N j distribution as faithfully as possible. These changes are detailed in Table 4 in Appendix A. We also modify the default Delphes card, producing two new cards that allow us to reproduce the reported expected events in Refs. [3] and [5]. These changes are detailed in Appendix A. The only difference between the two Delphes cards is the b-tagging efficiency. The four-top-quarks search [5] uses a higher b-tagging efficiency working point (77% average b-tagging efficiency) than the one used by the ttH search [3] (70% average b-tagging efficiency).
After simulating the events, we implement the event selection cuts and obtain the event yields. Each search has a different event selection for each channel as detailed in Tables 2  and 3. Additionally, we require at least 4 jets in the 2LSS ttH selection and at least 3 b-jets in the 2LSS and 3L four-top-quarks selection. The former modification is needed to compare our results to Fig. 2 in Ref. [3] and the latter is needed to obtain a signal-enhanced selection similar to the one defined by the use of the BDT in Ref. [5].
We also incorporate specific trigger selection efficiencies for each leptonic channel to the ROOT [45] code used to analyse the Delphes output. The simulated signal samples have been normalized using k-factors obtained from simulating similar events to NLO with the same set-up, and which are consistent with those in the literature [46].
After simulating a fixed set of events, we observe that the expected number of events from a given process in an analysis channel, N channel process , can be parametrized as follows 1 : where all the coefficients A ch i , B ch i are functions of M Z and absorb the acceptance of the channel and the cross-section for the specific process normalized to the corresponding coupling set to unity. The channels we consider are the different bins of Fig. 2 of Ref. [3] and the fourtop-quarks event yield with Q > 0 and with Q < 0 with the selection criteria of Ref. [5]. We obtain the A ch i and B ch i with Weighted Least Squares, where the uncertainty of each measurement is due to the Monte Carlo finite sampling, and with them we generate arbitrary points in the parameter space.
The likelihood fit to the eight bins of Fig. 2 of Ref. [3] is performed using the HistFitter package [47], which relies on RooFit [48] and the minimization algorithms from MINUIT [49]. Additionally, systematic uncertainties affecting the overall normalization of the SM processes are included in the fit as nuisance parameters (NP) with Gaussian constraints: 20% uncertainty is assigned to ttH, ttW , ttZ, and signal, and 50% uncertainty to diboson, while the various fake lepton components have uncertainties assigned corresponding to the normalisation factor precision reported in the ATLAS result.
In Figs. 3-5 we compute the impact of different points in parameter space for different M Z masses on the two experimental analysis. In the left column of these plots, we obtain the point that minimizes the Negative Log-Likelihood (NLL) [20] for the reported observed events in Fig. 2 of the ttH search [3], considering the pre-fit SM background contributions. After finding this minima, we plot the 1 and 2 standard deviations (s.d.) regions and the goodnessof-fit contour lines [50]. We include in this same left column the contour lines indicating the four-top-quark BSM and SM to SM only event ratio, N four-tops SM +BSM /N four-tops SM . These events need to pass the four-top-quarks-like selection cuts described in Table 3 and include SM four-top-quarks and BSM ttt + ttt, tttj + tttj and four-top-quarks processes. We can see the interplay between a good ttH fit, which requires asymmetry (non-negligible g ut ) but also a high b-jet multiplicity (large g tt ), and the four-top-quarks fit (not so large g tt ). We combine both measurements in the right column where we minimize the ttH+ four-top-quarks data For four-top-quarks we consider two N four-tops obs /N four-tops In all cases we find that the best fit points are compatible and we chose to focus on the AT-LAS result to tune our simulation and obtain further information. We plot the 1 s.d. and 2 s.d. regions for the ATLAS value and we showcase how the Z introduces an imbalance in the ratio of yields of ≥ 3 top-quarks events with positive and negative total leptonic charge, in the four-top-quarks search. Observe that the ATLAS four-top-quarks analysis [5] does not distinguish three-and four-top-quarks and therefore this imbalance is also induced by diagrams as in Fig. 1a because of an up-quark in the initial state and, as expected, grows with g ut . We plot the contour levels for r(4t) in the figures since it provides a qualitative insight on the behaviour of the asymmetry coming from the g ut -mediated BSM contributions in a four-top-quarks-like selection. We also observe that this parameter is not currently reported by the experimental collaborations, and could provide clues of BSM contributions.
To better identify interesting features, we plot in Fig 3 the region M Z < 2m t and g ct = 0; in Fig. 4 we study M Z > 2m t and g ct = 0; and in Fig. 5 we explore the g ct = 0 region while fixing g tt = 0.2 and 0.4. In the following paragraphs we discuss each one of these figures in detail. Although we find many points in parameter space compatible with the experimental results, we observe that M Z = 400 GeV has a slightly better accordance with the data. For M Z < 2m t (Fig. 3), the 1 s.d. regions in all plots include g ut ≈ 0. This corresponds to no charge asymmetry and is indicative of the fact that none of the points in the parameter space are a particularly good fit for the data in this region, as quantified by the poor goodnessof-fit, not being significantly better than that of the SM hypothesis (g ut = g ct = g tt = 0). This is because Z → tt is suppressed and thus the events providing charge asymmetry with large b-jet multiplicity are suppressed as well. This leaves Z → tu + tu as the dominant decay mode, worsening the fit. When comparing the right column (ttH and tttt) to the  Figure 4: Idem as in Fig. 3, but for M Z > 2m t . In this case, the opening of the decay Z → tt enhances the channel ug → tZ → ttt which induces a charge imbalance that favours considerably the fit in ttH. This can be appreciated in the goodness-of-fit (green dotted), which is considerably better in this case than for those in Fig. 3.
left column (only ttH), one should keep in mind that yields in four-top-quarks searches are obtained mainly through the three-top-quarks production diagrams, which are proportional to (g tt ) 2 BR(Z → tu + tu) and (g ut ) 2 BR(Z → tt * + t * t). For M Z = 200 GeV, the second process is negligible and therefore the effect of incorporating four-top-quarks to the NLL is to constrain g tt . Whereas for M Z = 300 GeV, there is a slight opening of Z → tt * + t * t and thus the 1 s.d. allowed region enlarges to the medium g ut region from Fig. 3c to 3d.  Figure 5: Idem as in Fig. 4, but for different points in the (g ut , g ct ) plane. The blue regions correspond to regions in parameter space disfavoured by D 0 ↔D 0 mixing as detailed in Sect. 2.
The goodness-of-fit to ttH data increases when we consider M Z > 2m t (Fig. 4). In fact, for these M Z values the Z → tt channel opens up and diagrams as in Fig. 1a provide multilepton events with charge asymmetry and large b-jet multiplicity, which are crucial to improve the fit. This important improvement can be explicitly seen by comparing the goodness-of-fit (green dotted contours) between the M Z < 2m t (Fig. 3) and M Z > 2m t (Fig. 4) results in the left plots. We can see in the left plots that as M Z increases from 400 GeV to 600 GeV, the tZ and ttZ production cross-sections decrease and therefore larger couplings are needed for the best-fit regions. We see however in these plots that the best-fit regions are in potential tension with four-tops-quark production cross-section (black dashed contours), which indicate a preference for lower values for g tt . In the right plots of Fig. 4 we include four-top-quarks data to the NLL and we observe that the best-fit point has a noticeable lower g tt and similar g ut . This is because four-top-quarks searches are more sensitive to g tt than to g ut . We do not consider masses above 600 GeV because we find that larger masses would yield H T distributions that enter into conflict with those reported in four-top-quarks analyses, as discussed below in Fig. 9.
We plot the g ct = 0 cases in Fig. 5 for M Z = 400 GeV and two values of g tt . In all cases we see that the best-fit point is in g ct = 0, which indicates that the main handles to accommodate the data are g ut and g tt . Nevertheless, we observe that g ct = 0 is allowed at the 1 s.d. level in a large part of parameter space. We show in blue the region disfavoured by D 0 ↔D 0 mixing which limits the g ut and g ct couplings as detailed in Sect. 3.2.
From these fits to the data we can conclude that the Z is not only compatible with the ttH and the four-top-quarks data, but also in many regions is more compatible that the SM. The tZ process can provide the necessary charge asymmetry and b-jet multiplicity while still being hidden in four-top-quarks production. To summarize how selecting a good benchmark point in parameter space reduces the tension in the reported results in ttH (Fig. 2 in Ref. [3]), we show in Fig. 6 how these results are modified if the Z BSM is added to the SM yields. We consider the best-fit point for ttH data only for M Z = 400 GeV: g ut = 0.04, g ct = 0.0 and g tt = 0.4. For comparison, we show the equivalent post-fit plot in the case where only SM processes are considered in the fit. In both fits, the same systematic model is used for the SM processes.
Although all masses above 2m t provide good fits, we consider the best-fit point for ttH and four-top-quarks data corresponding to M Z = 400 GeV and couplings g ut = 0.04, g ct = 0.0 and g tt = 0.2 as our benchmark point for further studies. As we show in Sect. 4, this mass has the interesting feature of having an H T distribution similar to, although slightly softer than, that of SM four-top-quark production.
Finally, it is interesting to observe that the preferred region in parameter space by the fit indicates that Z associated production either with a single top-quark or with a top-quark pair can mimic ttW ± in the ttH search and four-top-quarks in the four-top-quarks search. In Sect. 4, we study kinematic distributions and properties of the new signals compared to the main SM background processes in order to find ways to break the degeneracy.

Constraints from other observables
In this section we analyse how other observables besides ttH and tttt constrain the previously studied parameter space, namely 0 ≤ g ut ≤ 0.1, 0 ≤ g ct ≤ 0.1 and 0 ≤ g tt ≤ 1.0; and 200 ≤ M Z ≤ 600 GeV.
As a low-energy physics phenomena, D-meson mixing is particularly sensitive to any BSM effects. In particular, since virtual Z and top quarks can contribute to D 0 ↔D 0 mixing at loop level, it is a sensitive observable to our model. Since we approximate g uc ≈ 0, there is no tree level contribution to D-meson mixing. We find that the Lorentz structure in Eq. 1 provides a contribution to D-meson mixing that is translated to the D-meson mass difference as: where x = (M Z /m t ) 2 . The deduction of this expression with the explicit expression for f (x) and the numerical factors involved are detailed in Appendix B.
Although in principle one should compare the experimental value to the predicted value due to SM and BSM contributions, a detailed knowledge of ∆M SM D is currently lacking [20,51,52]. Therefore, as a naive estimation, we require that the BSM contribution to ∆M D is smaller than the uncertainty in its measurement: Using current available data [20] we obtain that the above constraints are translated into the model parameters as g ut g ct < 2.0 × 10 −3 to 4.5 × 10 −3 for M Z ranging from 200 GeV to 600 GeV, respectively. These bounds are plotted as blue regions in Fig. 5.
High-energy collider physics is sensitive both to inclusive on-shell production of the Z + X and to non-resonant behaviour, such as the re-scaling of different top-physics crosssections [22][23][24]27]. As mentioned before, taking g uu , g cc ≈ 0 makes our model insensitive to many of them. Other effects are greatly diminished.
Measuring top-quark rare decays opens a window to BSM effects. Z -induced top-quark rare decays have been studied in Ref. [53] and can be confronted to the experimental bounds set in Refs. [54][55][56][57]. Using the same matrix elements we can produce the relevant decay branching ratios for the model described in Eq. 1. After adding the smaller SM contributions, we find that rare decays branching ratios are below current limits, which are of O(10 −5 ). The main difference between our model and those considered in Ref. [53] is the absence of g uu and g cc couplings and of left-handed couplings, which effectively reduces the Z -induced rare decays and allows us to go to smaller masses.
As mentioned above, in high-energy collider searches our model yields no pp → Z production. In tt production, it can only contribute through the t-channel and cannot be probed by tt resonance searches. The only constraints coming from tt are then those on the inclusive crosssection, where we find that Z effects are well below the current experimental bounds [58, 59] for the parameter space we are exploring, being at most 45 fb for g ut = 0.1, g ct = 0.1, g tt = 1.0 and M Z = 200 GeV.
Our Z model is also sensitive to same-sign top-quark pair searches, both to prompt samesign top-quark pair production mediated by a Z in the t-channel and to same-sign top-quark production with an associated jet. The former is sensitive to (g ut ) 4 , (g ct ) 4 , and (g ut ) 2 (g ct ) 2 , the dominant process being uu → tt due to the PDF imbalance between the up quark and the others. The latter corresponds to the same tZ production mechanism we discuss in Sect. 2.1, but we now select events where the Z decays to tū + tc. Current experimental limits from prompt same-sign top-quark pair production are σ tt ≤ 1.2 pb [60] and σ tt ≤ 89 fb [61]. However, the latter limit is not fully model-independent, as the signal regions are optimized for higher mass vector mediators (with M V ≥ 1 TeV). As for our Z model, the highest possible σ tt we can produce is 76 fb, for M Z = 200 GeV and g ut = 0.1, which is just below the most stringent current experimental limits. M Z is so low that we should expect the acceptance to be significantly affected by the experimental cuts. In particular, the H T ≥ 750 GeV cut imposed in Ref.
[61] could be too high for such a low M Z . The cross-section σ tt is proportional to (g ut ) 4 and decreases as we increase the mass, e.g. σ tt = 15 fb for M Z = 400 GeV and g ut = 0.1. All things considered, we conclude that prompt same-sign top-quark pair production does not place relevant constraints on our parameter space. Ref.
[61] also casts experimental limits on ttj, although they consider a heavier Z . Taking into account the re-scaling due to the BR(Z → tj) and the same caveats about the acceptances, we find that these limits are also avoided by our model. A dedicated search for a light Z in this channel would be interesting, as one could potentially reconstruct the Z mass. Relaxing the H T cut and selecting events with a hard light jet would potentially augment the signal acceptance, while the background would still be relatively small. Lowmass tj resonances have been searched in the tt+ jets channel [62, 63], albeit with very low sensitivity.
While Z -mediated single top-quark production tj + tj is absent due to g uu , g cc ≈ 0, there is the possibility of Z radiative production tZ j, where Z is emitted from either the t/t or the jet (which can then be either u, u, c or c). However, in our model this channel is also severely suppressed due to the chiral nature of our coupling choice. All radiative tZ j production requires a W ± , either in the s-channel or the t-channel, with the latter dominating 2 . To interact with this W ± , we need left-handed quarks. However, in our model we require righthanded quarks to radiate a Z , which yields a considerable suppression. In particular, due to PDF imbalance the most relevant diagram with W ± for tZ j is ud → du → dtZ . But this diagram requires a left-handed up quark to interact with the W and a right-handed up quark to interact with the Z . As the up quark is essentially massless, the most important contribution to tZ j essentially vanishes for right-handed couplings only. To illustrate this point, assuming M Z = 400 GeV, we find that the tZ j cross-section for g L ut = 0.0, g R ut = 1.0 is approximately 40 times smaller than for g L ut = 1.0, g R ut = 0.0.     Figure 6: Comparison between ATLAS data [3] and a) SM background prediction b) BSM signal (tZ + ttZ ) plus SM background prediction for the event yields in 2LSS and 3L channels in categories based on the total lepton charge and the b-jet multiplicity. The BSM signal corresponds to a M Z = 400 GeV and couplings g ut = 0.04, g ct = 0.0 and g tt = 0.4, which minimizes the NLL for the ATLAS data [3] for M Z = 400 GeV. In both cases, we show the SM background prediction with fitted NP to the data although the fitted values differ due to the presence of BSM events in b).

Global search and kinematic features
In Sect. 3 we show how a specific Z vector boson can affect ttH and four-top-quarks analyses. Although all BSM diagrams to some extent affect both observables, the charge asymmetric contribution mimics mainly ttW ± , whereas the charge symmetric contribution mimics mostly four-top-quarks behaviour in different searches and channels. In this section we show how a more global study of the Z could help to break this degeneracy and disentangle signal from SM background processes. Then, based on the studied discriminating observables, we define signal-enriched regions that have either ttW ± or four-top-quarks as the main SM background and thus are sensitive to different parameters of the model. With the help of a few kinematical variables in these regions we aim to improve signal vs background discrimination, showing the potential of a future optimised Multivariate Analysis (MVA) discriminant.

Expected sensitivity of the global analysis
The event selection used for the global analysis is similar to that of the four-top-quarks search [5], but with some modifications. We define the regions as:

2LSS:
Two same-sign leptons, at least 3 jets and at least 1 b-jet, 3L: (Exactly) Three leptons, at least 2 jets and at least 1 b-jet.
No H T cut is applied in either of the regions.
The event categories in the global analysis are defined for the 2LSS and 3L selection separately, based on the discriminating variables: number of jets (N j ), number of b-jets (N b ), and total leptonic charge (Q lep ). In Fig. 7 we show the distributions for these discriminating variables in each category.
In order to perform a global analysis we define bins as (N j , N b , Q lep ) with Q lep = ± as a shorthand for ±2(1) for 2LSS (3L) and ≥ N j (N b ) meaning at least N j jets (N b b-jets). There are in total 34 (30) bins for the 2LSS (3L) channel. Low N j and N b bins will have a larger contamination of ttW ± and higher N j and N b bins will be more sensitive to four-top-quarks. The expected events for each one of these categories in the global analysis are shown in Fig. 8, where we use as a signal benchmark the combined ttH and four-top-quarks best-fit point for M Z = 400 GeV from Fig. 4. We also plot the main irreducible backgrounds simulated at NLO accuracy with the same generator settings as for the signal, with cross-sections scaled to those reported in Ref. [5]. From the global picture we see how the signal mimics ttW ± in low N j , low N b bins and four-top-quarks in high N j , high N b bins. There are some bins that provide a clear distinction between signal and background. More insight into how ttZ can resemble SM four-top-quarks is shown in Fig. 9, where both SM and BSM processes have a similar H T kinematic distribution in the events selected for the global histogram.
The global analysis shown in Fig. 8 can be used to estimate a potential discovery or the exclusion limit for the BSM Z model. This is done for M Z = 400 GeV and each point in the (g ut , g ct , g tt ) parameter space. Using Eq. 2 as detailed in Sec. 3, we generate global histograms for an arbitrary point in the parameter space using a fixed set of simulated points. These processes are not independent of each other and thus need to be considered simultaneously. An Asimov dataset [64] is used to calculate the expected significance and exclusion limits. The systematic uncertainties included in the global fit consist of 20% overall normalization uncertainty on the ttH, ttW , ttZ, four-top-quarks processes, and signal, as well as N j -and N b -dependent uncertainties on ttH, ttW , and ttZ to account for larger modelling uncertainties in the production of additional light/heavy-flavour jets. We show in Fig. 10 the expected significance and the expected exclusion limits for the case where M Z = 400 GeV and g ct = 0 for the corresponding integrated luminosity of the Run-2 dataset (139 fb −1 ) and of the expected Run-2 plus Run-3 datasets (300 fb −1 ).

Disentangling signal and background in signal-enriched regions
In order to define signal-enriched regions both for tZ and ttZ we require N b ≥ 3. Moreover, given the top-quark content in each of these processes we can use Q lep and N j to further separate them. We therefore define two signal-enriched regions: Region I is more sensitive to tZ while Region II is more sensitive to ttZ . We show the composition of the events that fulfil the selection of each Region in Fig. 11. We aim to examine observables that can disentangle signal from background in each one of the signal-enriched regions.
In the following, we plot only the fraction of events for the relevant processes in each region, i.e. ttW and tZ for Region I and four-top-quarks and ttZ for Region II. We show for Z the relevant masses 400 GeV and 600 GeV. We also find compelling to include in the study a BSM model in which we replace the spin-1 Z with a spin-0 scalar field Φ in the Lagrangian in Eq. 1, while coupling to the full chirality fermions u and t. For this scalar BSM model we use a mass M Φ = 400 GeV and the same couplings as the best-fit point of Z for this mass. The two relevant processes in this model are then tΦ and ttΦ.
To compare the relevant studied processes we calculate the separation between processes, defined as [65]: where f sig i (f bkg i ) is the signal (background) fraction of events in bin i.
We plot in Fig. 12 the H T distribution for both regions. We can see in both regions how the low mass (400 GeV) BSM scenarios mimic the SM distribution regardless of their spin (Z or Φ), whereas the 600 GeV case begins to show some deviation from the SM expected distribution.
The above results show that disentangling the 400 GeV BSM from the SM is challenging. Considering tZ and ttW events, we see that in the diagrams the lines connecting both samesign leptons have different Lorentz structure. We therefore explore the azimuthal angular separation between the leptons ∆φ( ± , ± ), which is a spin-correlation-sensitive observable. We observe in Fig. 13 that this observable can distinguish to some extent the different contributions in Region I, whereas no distinction occurs in Region II.
We can further observe that all leptons in tZ come from top quarks, whereas in ttW there is one lepton coming from a prompt W boson. In contrast to leptons coming from a top quark decaying to a W boson and then to a lepton, the lepton coming from a prompt W boson is not closely connected to a b-jet. Motivated by this, we propose a new kinematic variable defined as: MaxMin( , b) = The maximum of the minimum ∆R-distances between the same-sign leptons and a b-jet.
In the extreme case of boosted top quarks, we expect the tZ signal to have smaller MaxMin( , b) than the ttW background. For less boosted top quarks we expect this qualitative behaviour to still hold, although at a lesser extent. We show the MaxMin( , b) distribution for both regions in Fig. 14, where the expected behaviour is verified. Moreover, we find a slightly larger separation than for the ∆φ( ± , ± ) observable. We have verified that the MaxMin( , b) observable has better separation for larger H T , as expected from its construction.
The previous observables, ∆φ( ± ± ) and MaxMin( , b), have shown a separation power between tZ /Φ and ttW . We are interested in understanding how independent is the separation power of these observables. We show therefore in Fig. 15 a two-dimensional event histogram where we plot the distribution over both observables to distinguish between signal and background. To avoid low Monte Carlo statistic, we use the global analysis selection, without further splitting in Region I and II. There, we count events for all background processes and for the three different signal processes. We see from the figure that the two studied observables are not strongly correlated and can therefore be exploited simultaneously in a future MVA analysis.            Figure 12: H T fraction of events for both signal-enriched regions. As it can be seen, the 600 GeV BSM H T distribution has a separation an order of magnitude larger than the 400 GeV BSM scenarios for Region I, whereas for Region II all BSM scenarios are hardly distinguishable from SM. Region II Figure 13: Fraction of events of the azimuthal angular distance between the same-sign leptons for the main processes in Regions I and II. We find it useful to disentangle tZ /Φ from ttW in Region I.  Figure 14: Fraction of events of the MaxMin( , b) observable. We find this observable to be useful for region I which is dominated by tZ /Φ and ttW , since the latter has a lepton coming from a prompt W boson.  Figure 15: Simulated event yields for different processes in the two dimensional space formed by the azimuthal angular distance between the two same-sign leptons (x-axis) and the MaxMin( , b) observable (y-axis). When comparing the reported separations to those in Figs. 13 and 14, one should take into account that this is for the global selection, not Region I.

Summary and outlook
We have addressed the pattern of mild but persistent anomalies in multilepton plus b-jets events at the LHC from a phenomenological point of view. More precisely, we have studied a BSM model with a Z in a mass range 200 GeV to 600 GeV that couples hierarchically to the right-handed up-quarks u R , c R and t R . The hierarchy differentiates the third generation from the other two, allowing for tt diagonal couplings and FCNC tū and tc couplings. Although we briefly discuss possible UV completions for this model, we consider that a more in-depth analysis in this direction, probing the additional expected new processes, would be very interesting and compelling.
The motivation for this BSM model is based on the multilepton plus b-jets anomalies pointing towards excesses in four-top-quarks analyses, as well as excesses in the positivelycharged multilepton asymmetries resembling the ttW ± background.
Throughout this article we focus primarily on the ATLAS ttH analysis in Ref. [3] in which the post-fit yields of ttW ± show an excess over the expected cross-section by a factor of µ ttW = 1.39 +0.17 −0.16 . We study the case in which the ttW ± cross-section is not free-floating in the fit, allowing the BSM model to fill this excess. We explore this scenario by implementing the BSM model, simulating relevant processes up to detector level following the ATLAS analysis in two lepton same-sign (2LSS) and trilepton (3L) final states, and finding the BSM parameters that best fit the ATLAS results. We perform a new simultaneous fit to the observed ATLAS data in the eight bins corresponding to combinations of dilepton charge and b-jet multiplicity in 2LSS and 3L final states (Fig. 2 in Ref. [3]). The region in parameter space that best fits the ATLAS data is shown in Figs. 3-5. We show in Fig. 6 how the BSM model improves the interpretation of the ATLAS results in Ref. [3] without the need of re-scaling the ttW ± cross-section. These results indicate that for this observable the data would prefer over the SM a Z with mass 400 GeV to 600 GeV and couplings of the order g tt ∼ O(10 −1 ), g tu ∼ g tc ∼ O(10 −2 ). When the BSM model is also tested to reproduce the slight excess in the ATLAS four-top-quarks analysis in Ref. [5], the preferred region for g tt is reduced to approximately half its value, since otherwise the model would populate this analysis with too many events.
We have verified that the best-fit regions are safe to other observables such as D-meson mixing, Z -induced top quark rare decays, tt production, same-sign top-quark pair searches, tj resonance searches in ttj events, and Z -mediated single top-quark production.
In the second part of this paper we explore how the proposed BSM model could be distinguished from similar SM physics processes in the multilepton plus b-jets final state. In a first stage we define 34 (30) kinematic regions for the 2LSS (3L) channel that have different signal-to-background ratios and therefore could be exploited to disentangle them. We perform a simultaneous fit of the regions in this global analysis search to the Asimov dataset assuming 139 fb −1 and 300 fb −1 integrated luminosity, and determine the prospects for discovering or excluding the proposed BSM model (Fig. 10). In a second stage, we define two kinematic regions tailored to increase the fraction of tZ and ttZ events, respectively, and study observables that could help discriminating signal from the main SM background processes in each region. We find that H T could potentially discriminate tZ and ttZ for M Z above ∼ 600 GeV, whereas the H T distribution of ttZ is very similar to that of the SM fourtop-quarks process. We also find that the azimuthal separation between same-sign leptons ∆ϕ( ± , ± ) is also a good variable to separate signal from its main background. We also include the possibility of replacing the Z with a scalar Φ field and the same flavour structure, and find that the ∆ϕ( ± , ± ) observable has discriminating power to differentiate the three models: SM, Z and Φ. Finally, we also propose a new kinematic variable MaxMin( , b) (c.f. Eq. 8), which is likely to have a larger value if there is a lepton coming from a prompt W boson in comparison to a W boson coming from a top quark, which has a close b-jet as the top quark is boosted. This variable has a good separation power between signal and background. We find that it is easier to distinguish signal from background in Region I, in which tZ is the main signal. As a last test, we investigate for all events in the previous global analysis whether ∆ϕ( ± , ± ) and MaxMin( , b) are correlated by plotting their distribution in a 2D coloured histogram, and find them to have little correlation and thus considering them together improves the separation power in comparison to each observable on its own.
In summary, we have proposed a phenomenological FCNC Z model that couples hierarchically to the up-type right-handed quarks to explain LHC discrepancies in multilepton plus b-jet final states. We have found regions in parameter space that fit the data better than the SM and proposed different ways to explore the data to test the BSM Z model. We find that a sophisticated experimental search, along the lines of our proposed global analysis could in the near future shed light on the existence of such a BSM scenario.

A Pythia and Delphes parameters
As detailed in Sect. 3, we employ tuned parameters for Pythia 8 and Delphes. For the former, we use the Monash tune [44] and we change the parameters shown in Table 4 [44].
In the Delphes case, we tune a Delphes card for the ttH search [3] in order to match the expected ttW ± event yields with the reference cross-section σ = 727 fb. We use almost the same Delphes card for the four-top-quarks search [5], with the only difference being the b-tagging efficiency. The detailed features of the first Delphes cards are: • We set the electron identification efficiency as c-jet misidentification rate : 0.20 · tanh(0.02 · p T ) · 1 1 + 0.0034 · p T b-jet efficiency : 0.80 · tanh(0.003 · p T ) · 30 1 + 0.086 · p T • We modify the τ -tagging module to match the medium working point in Ref. [3] module TrackCountingTauTagging TauTagging  c-jet misidentification rate : 8.1 4.0 · 0.20 · tanh(0.02 · p T ) · 1 1 + 0.0034 · p T b-jet efficiency : 77.0 70.0 · 0.80 · tanh(0.003 · p T ) · 30 1 + 0.086 · p T This Delphes card allows us to recover the expected yields for ttW ± , four-top-quarks production, ttH and ttZ reported in Table 3 of Ref. [5] to a very good degree of approximation.
B Z induced D 0 − D 0 mixing We compute the one-loop contribution to D 0 − D 0 mixing due to a generic Z with arbitrary flavour-changing non-zero couplings g L,R ut,ct . We observe that the general tree level computation has already been performed in [67], whereas the one-loop contribution with equal couplings for both chiralities g L ut,ct = g R ut,ct can be found in [68]. Along this appendix we follow and expand the |∆C| = 2 calculations in previous references to the one-loop case with different couplings for each chirality.
The amplitude corresponding to the Feynman diagrams in Fig. 16 reads where P ut = g L ut P L + g R ut P R (10) P ct = g L ct P L + g R ct P R (11) and P L,R = (1 ∓ γ 5 )/2. We work in the low energy limit for u and c.
In order to perform the k-integral in Eq. 9 we observe that by symmetry arguments any integral with odd number of k's in the numerator vanishes. We are then left with the two relevant terms proportional to k k and m 2 t . We can then use the identity 1 A m B n = Γ(m + n) Γ(m) Γ(n) ∞ 0 dλ λ m−1 (λA + B) m+n (12) to perform k-integrals through the usual d-dimensions integrals such as setting d = 4 in all cases. After some algebra and some Fierz transformations, we obtain an explicit expression for M without any integral. By relating M terms to four-fermion effective operators Q 1...8 (see Ref. [67]) we obtain an expression for the effective Lagrangian L ef f . Using Q i = D 0 |Q i |D 0 and the modified vacuum saturation hypothesis as defined in [67], we obtain an explicit expression for D 0 |L ef f |D 0 . Utilizing the definition for we obtain the one-loop D-meson mixing parameter for chiral off-diagonal couplings with the top quark are the functions that appear when integrating in λ in the above scheme.
As it can be easily seen, one retrieves the results in Ref. [68] if Left and Right couplings are set equal. [63] CDF Collaboration, Search for a heavy particle decaying to a top quark and a light quark in pp collisions at √ s = 1.96 TeV, Phys. Rev. Lett. 108 (2012) 211805, [arXiv:1203.3894].