1 Introduction

The standard model (SM) is considered to be well established yet incomplete: It does not explain the puzzling structure of masses and mixings of elementary fermions. It also displays meta-stability in the Higgs sector, and ultimately loses control towards (very) high energy as the Higgs and hypercharge coupling run into Landau poles. It is therefore commonly accepted that the SM has to be extended into a more complete one, with guidance from both data and top-down theory frontiers.

The concept of asymptotic safety [1,2,3] opens up new directions [4,5,6,7] to build models that remain both fundamental and predictive at highest energies. Concrete models which extend the SM into asymptotically safe ones include vector-like fermions as well as additional scalars [7, 8], which allow for phenomenological signatures that can be probed at colliders. A crucial difference to common extensions of the SM is that the scalars, which are singlets under the SM gauge group \(SU(3)_C \times SU(2)_L \times U(1)_Y\), form a matrix in flavor space. This enhances the impact of Yukawa interactions and allows for new, flavorful signatures once the SM and BSM flavor sectors are connected.

In this work, we study concrete such models with new vector-like leptons (VLLs) and flavor portal couplings which link SM fermions to the new matter particles by renormalizable Yukawa interactions. Interestingly, while being asymptotically safe, or safe up to the Planck scale, such models allow to explain discrepancies between the SM and current data of anomalous magnetic moments (AMMs) [9,10,11]. Presently, the electron and the muon AMMs deviate from the SM by \(2.4~\sigma \) [12, 13] and \(3.5~\sigma \) [14], or 4.1 \(\sigma \) [15, 16], respectively. Note also recent debates [17,18,19,20]. In order to explain the current values of the AMMs, or values in a similar ballpark, the VLLs can be as light as a few hundred GeV, and thus can be probed at the LHC.

Early searches for VLLs at LEP excluded heavy leptons lighter than \(\sim 100~\)GeV [21]. At the LHC, ATLAS measurements excluded VLLs transforming as singlets under \(SU(2)_L\) in the range 114–176 GeV at 95% CL [22]. A recent CMS study based on \(77.4 \, \text {fb}^{-1}\) at 13 TeV searching for doublet VLLs coupling to third-generation leptons only [23] excluded VLLs in the mass range 120–790 GeV at 95% CL [24]. However, these VLL limits have been obtained within simplified models. In this work we study collider signatures of both the flavorful singlet and the doublet model in Refs. [9,10,11] in the multi-light lepton channel and confront them to the CMS search [24]. We further highlight how the specific features of the models such as lepton-flavor-violating-like decays suggest new observables with null test potential. These are worked out for the LHC full Run 2 150 \(\text {fb}^{-1}\) data set and high luminosity (HL)-LHC with 3000 \(\text {fb}^{-1}\) at \(\sqrt{s}=14\) TeV [25].

This paper is organized as follows: In Sect. 2 we present the BSM model frameworks and key parameters. In Sect. 3 we discuss the production and decay properties of the VLLs relevant for LHC phenomenology. We give the settings used for the event generation in Sect. 4. In Sect. 5 we compare distributions to the CMS measurements [24] to obtain constraints on the models’ parameter space. In Sect. 6 we construct new observables that allow for null tests of the SM and work out projections for the full Run 2 data set. Perspectives for the higher luminosity scenario of the HL-LHC are worked out in Sect. 7. We discuss signatures allowing for more general parameter regions beyond the AMMs in Sect. 8. In Sect. 9 we summarize. A comparison of resonance heights before and after including detector effects for Run 2 and the HL-LHC is provided in the Appendix.

2 BSM framework and setup

We start with the models of [9,10,11], which contain three generations of VLLs denoted by \(\psi _{L,R}\). These are either colorless \(SU(2)_L\) singlets with hypercharge \(Y=-1\) (singlet model) or colorless \(SU(2)_L\) doublets with \(Y=-1/2\) (doublet model). The \(SU(2)_L\)-components in the latter read \(\psi _{L,R} = (\psi ^0_{L,R}, \ \psi ^-_{L,R})^T.\) For the three generations of left-handed and right-handed SM leptons we use \(L=(\nu , \ \ell _L)^T\) and \(E=\ell _R\,\), respectively, and denote the Higgs doublet by H. All SM leptons and VLLs carry a lepton flavor index \(i=1,2,3\), which is often suppressed to avoid clutter. Both models also contain complex scalars \(S_{ij}\), with two flavor indices \(i,j=1,2,3\), which are singlets under the SM gauge interactions. In the interaction basis, the models’ BSM Yukawa sectors read

$$\begin{aligned} \begin{array}{l} \!\!{\mathcal {L}}^{\mathrm{{singlet}}}_{\mathrm{{Y}}} = -\kappa {\overline{L}}_i H \psi _{Ri} - \kappa '{\overline{E}}_i (S^{\dagger })_{ij}\psi _{Lj} - y\, {\overline{\psi }}_{Li} S_{ij} \psi _{Rj} + \mathrm {h.c.}, \\[1ex] \!\!{\mathcal {L}}^{\mathrm{{doublet}}}_{\mathrm{{Y}}} = -\kappa {\overline{E}}_i {H}^{\dagger } \psi _{Li} - \kappa '\, {\overline{L}}_i S_{ij} \psi _{Rj} - y\, {\overline{\psi }}_{Li } S_{ij} \psi _{Rj} + \mathrm {h.c.},\nonumber \end{array}\\ \end{aligned}$$
(1)

where the contraction of gauge indices is assumed. Here, we followed [9] and identified SU(3)-flavor symmetries of the leptons with ones of the VLLs. This identification has important consequences for phenomenology: Each lepton flavor is conserved and leptons couple universally within (1), and the BSM Yukawas \(y, \kappa , \kappa ^\prime \) become single couplings, instead of being tensors. While y is key in variants of the asymptotically safe framework [4, 7], in models like (1) with mixed SM-BSM Yukawas its presence is not required to achieve a controlled UV-behavior [9]. As in addition the phenomenological implications of y are less relevant we do not consider it in the numerical analysis.

After spontaneous symmetry breaking, the vector-like fermions and the leptons mix, see [9, 10] for details. To be specific, we denote the lightest three mass eigenstates by leptons and the others by VLLs, and continue to use the notation as introduced above. \(Z\rightarrow \ell \ell \) data [14] constrains the mixing angles \(\theta \) of left-handed (right-handed) leptons in the singlet (doublet) model as \(\theta \simeq \kappa v_h/\sqrt{2}M_F < {\mathcal {O}}(10^{-2})\), see [26] for recent fits. Here, \(v_h \simeq 246 \) GeV is the Higgs vacuum expectation value (vev), and we denote by \(M_F\) the common mass of all flavor and \(SU(2)_L\)-components of the VLLs. We learn that \(\theta , \kappa \ll 1\), which allows for a small angle approximation. At first order in \(\kappa \) and \(\kappa ^\prime \), the interactions in the mass basis in the singlet model read

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\text {int}}^{\text {singlet}}=&{\, -e\, {\overline{\psi }}\gamma ^\mu \psi A_\mu } {+\, \frac{g}{\cos \theta _w} \, {\overline{\psi }}\gamma ^\mu \psi Z_\mu } \\&+ \Big [-\frac{\kappa }{\sqrt{2}}\, {\bar{\ell }}_L \psi _R\, h -\kappa ^\prime {\bar{\ell }}_{R} S^\dagger \psi _L +g_S\, {\bar{\ell }}_R S^\dagger \ell _L \\&+ g_Z \, {\overline{\ell }}_L\gamma ^\mu \psi _L\, Z_\mu + g_W\, {\overline{\nu }} \gamma ^\mu \psi _L \, W_\mu ^+ + \text {h.c.}\Big ], \end{aligned} \end{aligned}$$
(2)

where \(A_\mu \) denotes the photon, h corresponds to the physical Higgs boson with \(M_h = 125\) GeV and \(e,g,\theta _w\) are the electromagnetic coupling, the \(SU(2)_L\) coupling and the weak mixing angle, respectively. The remaining couplings fulfill

$$\begin{aligned} \begin{aligned} g_S = \frac{\kappa ^\prime \kappa }{\sqrt{2}} \frac{v_h}{M_F},\ \ \ g_Z = {-\frac{\kappa g}{2\sqrt{2}\cos \theta _w} \frac{v_h}{ M_F}}, \ \ \ g_W = {\frac{\kappa g}{2} \frac{v_h}{M_F}} .\nonumber \end{aligned}\\ \end{aligned}$$
(3)

For the doublet model, we find

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_\text {int}^{\text {doublet}}=&{\, -e\, \overline{\psi ^{-}}\gamma ^\mu \psi ^- A_\mu } +\, \frac{g}{2\cos \theta _w}\Big [ \overline{\psi ^0} \gamma ^\mu \psi ^0 \\ {}&+ (2\sin ^2\theta _w - 1) \overline{\psi ^{-}}\gamma ^\mu \psi ^{-} \Big ]Z_\mu \\&{ + \Bigg [\frac{g}{\sqrt{2}}\, \overline{\psi ^{-}} \gamma ^\mu \psi ^0 W_\mu ^-} -\frac{\kappa }{\sqrt{2}}\, {\bar{\ell }}_R \psi _L^{-}\, h \\ {}&-\kappa ^\prime {\bar{\ell }}_{L} S \psi _R^- \ -\kappa ^\prime {\bar{\nu }} S \psi _R^0 +g_S\, {\bar{\ell }}_L S \,\ell _R \\&+ g_Z \, {\overline{\ell }}_R\gamma ^\mu \psi ^-_R\, Z_\mu + g_W\, {\overline{\ell }}_R \gamma ^\mu \psi _R^0 \, W_\mu ^- + \text {h.c.}\Bigg ], \end{aligned} \end{aligned}$$
(4)

with couplings

$$\begin{aligned} \begin{aligned} g_S = \frac{\kappa ^\prime \kappa }{\sqrt{2}} \frac{v_h}{M_F}, \ \ \ g_Z = {\frac{\kappa g}{2\sqrt{2}\cos \theta _w} \frac{v_h}{ M_F}}, \ \ \ g_W = {- \frac{\kappa g}{2} \frac{v_h}{M_F} } .\nonumber \\ \end{aligned}\\ \end{aligned}$$
(5)
Fig. 1
figure 1

Dominant pair (diagrams a and b) and single (diagrams c and d) production channels of vector-like leptons at pp colliders. Diagrams c and d involve the couplings \(g_Z\) and \(g_W\), which are induced by fermion mixing and \(\kappa \) (2), (4). In diagram d, the final states \(\psi ^-{\overline{\nu }}\) (\(\ell ^-{\overline{\psi }}^0\)) are only possible in the singlet (doublet) model

Full size image

The vertex \({\overline{\nu }} \gamma ^\mu \psi _L^- \,W_\mu ^+ \) arises only at higher order, see [9] for details, and can be safely neglected for the purpose of this analysis. The parameters we are concerned with in the two models are therefore \(M_F, \kappa , \kappa ^\prime \) and the common mass of the BSM scalars, \(M_S\). Addressing the muon AMM anomaly \(\varDelta a_{\mu } \equiv a_{\mu }^{\text {exp}} - a_{\mu }^{\text {SM}} = \ 268(63)(43)\times 10^{-11}\) [14] at one loop eliminates \(\kappa ^\prime \) in terms of the BSM masses by [9, 10]

$$\begin{aligned} \varDelta a_{\mu } =\frac{\kappa '^2}{32\pi ^2} \, \frac{m_{\mu }^2}{M_{F}^2} \, f\left( \frac{M_S^2}{M_F^2}\right) , \end{aligned}$$
(6)

with \(f(t) = (2 t^3+3 t^2-6 t^2 \ln t-6 t+1)/(t-1)^4\) positive for any t, and \(f(0)=1\). For example, for \(M_S=500\) GeV and \(M_F = \{100, 500, 1000 \}\) GeV one obtains \(\kappa ' \simeq \) \(\{3.6, 6.5,\) \( 10.4\}\), respectively. The coupling \(\kappa \) is needed to account for the electron AMM anomaly \(\varDelta a_e\) [12, 13], however, it can be chosen with some freedom since parameters of the scalar sector also play a role here, see [9, 10] for details. For simplicity, we fix \(\kappa /\kappa ^\prime = 10^{-2}\), consistent with Z-decay constraints and both AMMs.

Let us briefly comment on the scalar sector of the BSM framework [9, 10, 27]. The presence of H and the flavorful scalar matrix field \(S_{ij}\) allows for a substantial scalar potential with in total three quartic couplings plus a Higgs portal one \(\delta S^\dagger S H^\dagger H\). In addition to successful electroweak symmetry breaking, the diagonal entries of S can acquire a non-trivial vev, \(v_s\). Interestingly, two different configurations exist: A universal ground state in which diagonal entries have the same vev, and one in which the vev points in a single flavor direction, breaking universality spontaneously. Both the portal \(\delta \) and \(v_s\) are instrumental in achieving the chirally enhanced 1-loop contributions explaining AMMs. However, as discussed in the next section, the impact of \(\delta \), \(v_s\) on the present collider study is negligible, and we do not consider them in this work.

To summarize, in the following we perform a collider analysis in the two models, one with three flavors of singlet VLLs (2), and one with three flavors of doublet VLLs (4), featuring nine flavored scalar singlets and with parameters

$$\begin{aligned} M_S, M_F \, , \kappa /\kappa ^\prime =10^{-2}\, , \kappa ^\prime =\kappa ^\prime (M_S,M_F) \, . \end{aligned}$$
(7)

Using (6) together with the muon AMM data to express \(\kappa ^\prime \) in terms of \(M_S\) and \(M_F\) renders the numerical predictions for production and decay of the BSM sector in terms of the latter two masses. We discuss more general settings in Sect. 8.

The reduction of the models’ BSM parameter space (i.e. BSM masses, Yukawas, and quartic couplings) onto the set (7) is sufficient to catch the leading degrees of freedom for a pp-collider study, and to validate the models up to the TeV energy range. To further demonstrate that all couplings reach the Planck scale without poles or instabilities requires a complete renormalization group analysis. This has been done previously for a wide range of BSM parameters using \(M_F=2 M_S=1\) TeV [9, 10]. In general, not every point in the BSM parameter space is guaranteed to be Planck safe. Still, since the widening of the mass range towards \(M_S, M_F\) within 0.1 to 1 TeV has only a minor effect on the RG running, we can find Planck safe trajectories within suitable ranges for the remaining BSM Yukawas and quartics, using the methods of [9, 10]. This completes the discussion of our setup.

3 LHC production and decay

At the LHC, VLLs can be produced in pairs (upper plots) or singly (lower plots) in quark fusion through electroweak interactions shown in Fig. 1. Pair production occurs through s-channel photon or Z (Fig. 1a), and in the doublet model additionally through s-channel W exchange (Fig. 1b). Single production is induced by the Yukawa portal coupling \(\kappa \) through fermion mixing and Z, W-exchange (Fig. 1c, d). All three flavors are produced universally. Additional contributions to VLL production arise through s-channel Higgs and BSM scalars \(S_{ii}\) induced by Higgs-scalar mixing (not shown). Due to both quark-Yukawa and parton-luminosity suppression these contributions to matrix elements are suppressed by at least two orders of magnitude with respect to electroweak contributions and thus negligible. Further production channels through Yukawa interactions open up at lepton colliders, briefly discussed in [9].

In Fig. 2 we show pair- and single-production cross sections for a single species \(\psi _i\) – with lepton flavor index i fixed – at the LHC with \(\sqrt{s}=13\) TeV as a function of \(M_F\) for \(M_S = 500\) GeV and the procedure in (7). In both the singlet (left) and doublet model (right) the pair-production cross section is roughly three orders of magnitude larger than the single production cross section. This is due to the fact that single production is only induced by mixing between SM leptons and VLLs, while dominant pair production of VLLs at the LHC occurs through electroweak gauge interactions. \(\kappa ^\prime \), the larger of the BSM Yukawa couplings, is irrelevant also for single production, but turns out to be important for BSM sector decays.

In the singlet model, the decay rates of the possible decay channels of the VLLs are

$$\begin{aligned} \begin{aligned} \varGamma (\psi _i \rightarrow h \ell _i ^-)&=\kappa ^2 \frac{M_F}{64 \pi } (1 - r_h^2)^2 , \\ \varGamma (\psi _i \rightarrow S^{\,*}_{ij}\, \ell ^-_j)&= \kappa ^{\prime \,2} \frac{M_F}{32 \pi } (1 - r_S^2)^2 ,~~~~~( j ~\text{ fixed}) \\ \varGamma (\psi _i \rightarrow W^- \nu _i )&= g_W^2 \frac{M_F}{32 \pi } (1 - r_W^2)^2(2+1/r_W^2),\\ \varGamma (\psi _i \rightarrow Z \ell ^-_i )&= g_Z^2\, \frac{M_F}{32 \pi } (1 - r_Z^2)^2(2+1/r_Z^2), \end{aligned} \end{aligned}$$
(8)

where \(r_X = M_X/M_F\). For large values of \(\kappa '\) the decay \(\psi _i \rightarrow S^{\,*}_{ij}\, \ell ^-_j\) dominates if kinematically allowed, as seen in Fig. 3 (left). Quantitatively, in the large-\(M_F\) limit, the decays through the \(S_{ij}\) dominate over Higgs-mediated decays (decays through weak bosons) for \(\kappa '\gtrsim \kappa /\sqrt{6}\) (\(\kappa '\gtrsim \kappa / \sqrt{3}\)). In the doublet model, we obtain the decay rates

$$\begin{aligned} \begin{aligned} \varGamma (\psi _i^- \rightarrow h \ell _i ^-)&=\kappa ^2 \frac{M_F}{64 \pi } (1 - r_h^2)^2 , \\ \varGamma (\psi _i^- \rightarrow S_{ji}\, \ell ^-_j)&= \kappa ^{\prime \,2} \frac{M_F}{32 \pi } (1 - r_S^2)^2 ,~~~~~( j ~\text{ fixed}) \\ \varGamma (\psi _i^0 \rightarrow S_{ji}\, \nu _j)&= \kappa ^{\prime \,2} \frac{M_F}{32 \pi } (1 - r_S^2)^2 ,~~~~~( j ~\text{ fixed}) \\ \varGamma (\psi _i^- \rightarrow Z \ell ^-_i )&= g_Z^2\, \frac{M_F}{32 \pi } (1 - r_Z^2)^2(2+1/r_Z^2),\\ \varGamma (\psi _i^0 \rightarrow W^+ \ell ^-_i )&= g_W^2 \frac{M_F}{32 \pi } (1 - r_W^2)^2(2+1/r_W^2).\\ \end{aligned} \end{aligned}$$
(9)

The corresponding branching ratios of the \(\psi ^-\) and the \(\psi ^0\) decays are shown in Fig. 3 (right). As in the singlet model, the decays to BSM scalars dominate for large \(\kappa '\) if allowed by the mass hierarchy of the BSM sector. As already stated in the previous section we assume that the \(\psi ^-\) and \(\psi ^0\) are degenerate in mass; we therefore neglect small isospin splitting induced by electromagnetic interaction \(\varDelta m = M_{\psi ^{-1}} -M_{\psi ^{0}}=g^2/(4 \pi ) \sin ^2 \theta _W M_Z/2 \simeq 0.4\) GeV, that also allows for rare inter-multiplet decays \(\psi ^- \rightarrow \psi ^0 W^{- *}\). The smallness of the splitting prohibits that for instance searches in R-parity violating SUSY models into four light leptons [28] apply to the VLL models.

Fig. 2
figure 2

Cross sections for \(\psi _i\) pair production (top) and single production (bottom) at \(\sqrt{s} = 13~\) TeV for different vector-like lepton masses in the singlet model (left) and the doublet model (right) for \(M_S=500\) GeV and the procedure described in (7)

Full size image

The decays of VLLs to \(S_{ij}\) plus lepton, with flavorful scalars \(S_{ij}\), are a singular feature of the models, which distinguishes them from other theories with VLLs, such as [23, 29,30,31,32]. Moreover, the \(S_{ij}\) can decay through fermion mixing to lepton final states, in which they can be searched for. Specifically, the singlet model features the cascade decays

$$\begin{aligned} \psi _i \rightarrow S_{ij}^*\, \ell ^-_j \rightarrow \ell _i^- \, \ell _j^+ \, \ell ^-_j. \end{aligned}$$
(10)

Similarly, in the doublet model the decays of the charged and neutral VLLs proceed as

$$\begin{aligned} \begin{aligned}&\psi _i^- \rightarrow S_{ji}\, \ell ^-_j \rightarrow \ell _i^- \, \ell _j^+ \, \ell ^-_j,&\psi _i^0 \rightarrow S_{ji}\, \nu _j \rightarrow \ell _i^- \, \ell _j^+ \, \nu _j. \end{aligned} \end{aligned}$$
(11)

These processes preserve flavor; however, the scalar decay yields a dilepton pair with different-flavor charged leptons for \(i\ne j\), which looks as if lepton flavor has been violated and cleanly signals new physics. While the scalars may also decay to dibosons through triangle loops, or to two VLLs through the coupling y if \(M_S< 2M_F\), see (1) and [9] for details, here we assume that these rates are negligible.

Fig. 3
figure 3

Branching ratios of the on-shell decays of the VLLs as a function of their mass in the singlet model (left) and the doublet model (right) for \(M_S =500\) GeV and (7). Larger ratios \(\kappa /\kappa '\) would enhance the branching ratios of the electroweak decays

Full size image

In this work, we are interested in final states with at least four light leptons (4L), where a light lepton is an electron or a muon, as in [24]. When the \(\psi _i\) are pair-produced and decay through Eqs. (10) and (11), only certain flavor final states of each single decay can contribute to a 4L final state. These are given in Table 1. For instance, pair production of negatively charged third-generation VLLs and their subsequent decay through the \(S_{ij}\) can only give rise to 4L final states through the decays

$$\begin{aligned} \begin{aligned} \psi _3^{(-)}{\overline{\psi }}_3^{(-)}\ \rightarrow \ \&e^-e^+e^-e^+\tau ^-\tau ^+,\ \ e^-e^+\mu ^-\mu ^+\tau ^-\tau ^+,\\&\mu ^-\mu ^+\mu ^-\mu ^+\tau ^-\tau ^+. \end{aligned} \end{aligned}$$

A more detailed account of the flavor configurations of the 4L final states, including those arising from weak and Higgs-mediated decays, is provided in Sect. 5.1. Notice that the decay chains (10), (11) allow to observe resonance structures from \(S_{ij}\)-decays in clean different-flavor dilepton invariant mass distributions if the VLLs are sufficiently heavy, \(M_F > M_S\). We exploit this possibility in Sect. 6.

4 Event simulation

In this section, we describe the procedure used to generate a sample of events with 4L final states at the LHC.

We employ FeynRules [33] to compute the Feynman rules at leading order (LO) for the models in Eqs. (2) and (4). The particles and Feynman rules are then implemented into UFO models [34]. These UFO models are interfaced to the Monte Carlo generator MadGraph5_aMC@NLO [35] to compute production cross sections of single and pair production of VLLs at the LHC as well as distributions of final state particles at parton level. For each process we generate \(5\times 10^4\) events. The decay of particles is handled with MadSpin [36]. For the event generation the NNPDF3.0 [37] PDF set is used. PDF and scale variation uncertainties are computed within MadGraph5_aMC@NLO for each of the PDF sets. The scale variation uncertainties are computed by varying factorization and renormalization scales independently between \(0.5\mu _0 \le \mu _F/R \le 2\mu _0\), where the scale \(\mu _0\) is computed in MadGraph5_aMC@NLO with different schemes. For the theory uncertainty we add PDF uncertainties, scale variation uncertainties and scheme variation uncertainties in quadrature.

Table 1 Decay modes of the VLLs through the S scalars, (10) and (11), giving rise to a final state with four light leptons after \({\bar{\psi }}^{(-)}_i \psi ^{(-)}_i\) or \({\bar{\psi }}^0_i \psi ^-_i\) pair production. For the third generation \(\psi _3^0\) no corresponding 4L final states arise
Full size table

For the event generation we adapt settings similar to the CMS study [24]. We focus on the final states with at least four light leptons, that is, muons and electrons and require the missing transverse momentum, \(p_T^\text {miss}\), to be smaller than 50 GeV. This cut serves to resemble the signal region considered by CMS and to suppress contributions from neutrinos in the decay of the electroweak bosons. Electrons and muons are required to have a minimum transverse momentum of \(p_T^\ell \ge 20\) GeV. Similarly to the CMS analysis, we neglect all events with a light-lepton invariant mass, \(m_{\ell \ell }\), smaller than 12 GeV for all flavor and charge combinations. This cut serves to suppress resonances in the low-mass region. For the event generation we fix \(\kappa =10^{-2}\kappa ^\prime \). Masses of the new scalars and VLLs are varied between \(M_S=300 {-} 1200\) GeV and \(M_F = 100{-}1000\) GeV with \(\kappa '\) computed according to Eq. (6). We also consider ttZ, triboson and ZZ production as these processes contribute to the SM background for the distributions studied in [24]. The ZZ production includes contributions from virtual photons via \(p p \rightarrow \gamma ^* \gamma ^*,\gamma ^* Z\). ZZj final states are included via multijet merging in PYTHIA8 [38]. We also take into account in the cross section gluon-fusion contributions \(g g \rightarrow ZZ\), where the lowest order is induced at 1-loop. SM background processes are computed at LO within MadGraph5_aMC@NLO using the same set of PDF sets. Higher order corrections to SM production cross sections are taken from literature [39,40,41,42,43,44,45,46] and are taken into account by applying k factors to the LO distributions. To perform a simulation of the detector response we shower and hadronize the events with PYTHIA8 and use DELPHES3 [47] for the fast detector simulation, yielding events at particle level. Jets are clustered with the anti-\(k_t\) algorithm [48] with a radius parameter \(R=0.5\) applying the FastJet package [49]. All criteria for the analysis are taken from the CMS default card for simplicity. In Table 2 we summarize the values for the parameters and signal selection cuts used in the event generation and detector simulation.

Table 2 Parameters used in the event generation, detector simulation and the reconstruction algorithm
Full size table
Fig. 4
figure 4

Examples of signal channel Feynman diagrams with at least four light leptons in the final state in the singlet (2) and doublet (4) model. Only first- and second-generation vector-like leptons can contribute via the diagram with jets in the final state (a) and the single production diagram (c), while all generations can contribute to diagram (b)

Full size image

5 Constraints from CMS data

Here we confront the models (2), (4) to the CMS search [24] using the 4L final state, which is expected to be the channel most sensitive to contributions stemming from three generations of VLLs. In Sect. 5.1, we study decay chains into 4L final states and their multiplicities. In Sect. 5.2 we compare the distributions of the scalar sum of transverse momenta of the four light leptons \((e, \mu )\) with the largest transverse momenta, \(L_T\), with CMS data to obtain constraints on BSM masses.

5.1 4L multiplicities

4L final states stem from both single and pair production of VLLs. Due to the flavor structure of the BSM sector, the following decay chains include 4L final states in the singlet model:

$$\begin{aligned} \begin{aligned} p p&\rightarrow \psi _i {\bar{\psi }}_i \rightarrow \ell _i^- \ell _i^+ \ell _j^+\ell _j^- \ell _k^+\ell _k^- \quad \hbox {for}\ i,j,k=1,2,3 , {\ (20)}\\ p p&\rightarrow \psi _i {\bar{\psi }}_i \rightarrow \ell _i^- \ell _i^+ q_j {\bar{q}}_j \ell _k^+\ell _k^- \quad \hbox {for}\ i,k=1,2, {\ (15\times 4)}\\ p p&\rightarrow \psi _i {\bar{\psi }}_i \rightarrow \ell _i^- \ell _i^+ \ell _j^+\ell _j^- \nu _k {\bar{\nu }}_k \quad \hbox {for}\ i,j=1,2 , k=1,2,3 ,{\ (12)}\\ p p&\rightarrow \psi _i {\bar{\psi }}_i \rightarrow \nu _i \ell _i^+ \ell _j^+\ell _j^- \ell _k^- {\bar{\nu }}_k \quad \hbox {for}\ i,j,k=1,2, {\ (8)}\\ p p&\rightarrow \psi _i {\bar{\psi }}_i \rightarrow \ell _i^- {\bar{\nu }}_i \ell _j^+\ell _j^- \ell _k^+ \nu _k \quad \hbox {for}\ i,j,k=1,2,{\ (8) }\\ p p&\rightarrow \psi _i \ell _i^+ \rightarrow \ell _i^-\ell _j^+\ell _j^- \ell _i^+ \quad \hbox {for}\ i,j=1,2, {\ (4)}\\ p p&\rightarrow {\bar{\psi }}_i \ell _i^- \rightarrow \ell _i^+\ell _j^+\ell _j^- \ell _i^- \quad \hbox {for}\ i,j=1,2, {\ (4)} \\ \end{aligned} \end{aligned}$$
(12)

where i, j, k are flavor indices and \(q_i={u,d,c,s,b}\). We also indicate the values that the lepton flavor indices can take, and between parentheses the number of 4L final states of each chain after summing over all indices. Note that, for the first decay chain in (12), final states with four light leptons occur only when at most one of the three indices i, j, k is equal to 3. For explicit expressions of flavors in the decays see Table 1.

In the doublet model, the negatively charged state \(\psi ^-\) decays into 4L final states as in (12) with the exception of the decays with 8-fold multiplicity. These correspond to W-mediated decays \(\psi _i^- \rightarrow \nu _i \ell ^-_j {\bar{\nu }}_j\), which are subleading in the doublet model. Additionally, when the \(\psi ^0\) are produced, 4L final states arise through

$$\begin{aligned} \begin{aligned} p p&\rightarrow \psi _i^0 {\overline{\psi }}^0_i \rightarrow \nu _j {\overline{\nu }}_k \ell _j^+\ell _i^- \ell _i^+\ell _k^- \quad \hbox {for}\ i,j {,k }=1,2, {\ (8)}\\ p p&\rightarrow \psi _i^- {\overline{\psi }}^0_i \rightarrow \ell _i^- \ell _i^+ {\overline{\nu }}_j \ell _j^- \ell _k^+\ell _k^- \quad \hbox {for}\ i,j{,k }=1,2,{\ (8)}\\ p p&\rightarrow \psi _i^0 {\psi }^+_i \rightarrow \ell _i^- \ell _i^+ \ell _j^+\nu _j \ell _k^+\ell _k^- \quad \hbox {for}\ i,j{,k }=1,2,{\ (8)}\\ p p&\rightarrow \psi _i^- {\overline{\psi }}^0_i \rightarrow \ell _i^- \ell _i^+ {\overline{q}}_j q_j \ell _k^+\ell _k^- \quad \hbox {for}\ i{,k }=1,2,{\ (15\times 4)}\\ p p&\rightarrow \psi _i^0 {\psi }^+_i \rightarrow \ell _i^- \ell _i^+ {\overline{q}}_j q_j \ell _k^+\ell _k^- \quad \hbox {for}\ i{,k }=1,2.{\ (15\times 4)}\\ \end{aligned} \end{aligned}$$
(13)

The first decay chain in Eq. (12), involving a six charged-lepton final state, is the only one where production of the third generation \(\psi _3\) can give rise to a 4L final state. In all other cases, \(\psi _3\) production yields at most three light leptons in the final state, since a \(\tau ^+\tau ^-\) pair is always produced due to flavor conservation. In Fig. 4 we give examples of Feynman diagrams for the different decay chains, with jets (a) or without them (b), and from single production (c).

In Fig. 5 we show the cross section at the \(\sqrt{s} = 13\) TeV LHC for BSM production of at least four light leptons in terms of the VLL mass for the singlet model (left) and the doublet model (right) for \(M_S=500\) GeV, together with cross sections of the models in Ref. [23]. In general, our cross sections are larger by roughly two orders of magnitude. This enhancement stems from the first and second generation of VLLs, which present a richer multiplicity of decays into 4L final states than the \(\psi _{3}\). For \(M_F < M_S\), the enhancement originates predominantly from the additional final states with two jets and four light leptons, while for \(M_F > M_S\) cross sections increase further, up to a factor of approximately \(10^4\). This effect is caused by the VLLs decaying mainly through on-shell production of scalars \(S_{ij}\), see Fig. 3, and their subsequent decay into light leptons. In Fig. 5, the spikes in the 4L cross sections at \(M_F\sim M_S\) indicate resonant decays of the VLLs into the new scalars, signaling the onset of the on-shell decays of Eqs. (10) and (11).

Fig. 5
figure 5

Cross section for BSM production of at least four light leptons at a pp collider with \(\sqrt{s} = 13\) TeV in the singlet (left) and the doublet (right) model as a function of the VLL mass for \(M_S=500\) GeV. The red curves correspond to the VLL models (2) and (4), while the blue curves correspond to third-generation VLL models as in [23]. The band widths include uncertainties discussed in Sect. 4. The spikes in our VLL models at \(M_F\sim M_S\) arise due to the onset of on-shell decays of the VLLs through the new scalars

Full size image

5.2 \(L_T\) distributions and CMS constraints

Fig. 6
figure 6

Allowed (green and yellow points) and excluded (purple points) values of the VLL mass \(M_F\) and the BSM scalar mass \(M_S\) with \(\kappa '\) fixed (7). For the points marked as allowed, all bins in the sum of transverse momenta \(L_T\) of the 4L final sates fall within 1\(\sigma \) of central values measured by CMS [24]. For benchmark points marked as yellow circles we show the \(L_T\) distributions in Fig. 7. Above the green dashed curve \(\kappa ^\prime \) becomes non-perturbative

Full size image
Fig. 7
figure 7

\(L_T\) distributions in the singlet (left) and the doublet model (right) for SM background processes in our simulation (green shaded area) and for the different benchmark masses of vector-like fermions and new scalars (yellow circles in Fig. 6). The observables are shown for an integrated luminosity of 77.4 fb\({}^{-1}\) and subsequent detector simulation. Also shown are CMS data [24] (black points), including the range covered up to \(1 \sigma \) (hatched area), see text for details

Full size image

CMS has searched for VLLs employing the scalar sum of the leading four light leptons’ transverse momenta, \(L_T\), finding no significant discrepancies with the background [24]. To work out the implications of this analysis for the models (2) and (4) we compute the \(L_T\) distributions for different values of \(M_S\) and \(M_F\) and fixed BSM Yukawas (7). After performing the detector simulation we compare the distributions to CMS data for 4L final states. We also compute \(L_T\) distributions for the dominant SM background processes of ZZ, triboson and \(t{\bar{t}}Z\) production. We include the control region veto, two dilepton pairs with invariant masses \(76~\mathrm{GeV}<m_{2\ell }<106~\)GeV, and set the bin width to 150 GeV as in [24].

Since our simulation of the SM background is performed at LO and only a fast detector simulation is publicly available, differences with the one by CMS are expected. In contrast to CMS, we can not perform a fit of the background distribution to a control region. This prohibits a quantitative reinterpretation of the data at precision level, which could be obtained from an actual experimental analysis only. Still, we find that our background simulation is in reasonable agreement with the shape and the bin content of the \(L_T\) distribution. In view of the differences between our SM prediction and the one from CMS, and to make progress, in the following we refer to benchmarks as ’excluded’ if the BSM distribution overshoots the CMS data in at least one of the bins by more than one sigma.

Our findings are summarized in Fig. 6, showing which masses are compatible with data (green and yellow circles) and which are not (purple circles) for the singlet (left) and the doublet (right) model. We scanned 40 points in the singlet and 20 in the doublet model, and expect these to indicate the main features of the parameter space in the \(M_S,M_F\)-plane. The purple hatched region is excluded, while the remainder is still to be probed. Fig. 6 also shows for which masses the coupling \(\kappa ^\prime \) required to accommodate the present \((g-2)_\mu \) anomaly becomes non-perturbative (above the green dashed curve).

We observe the following pattern: For both models the region \(M_F\sim M_S\) is excluded. This is a result of the 4L cross sections’ enhancement around the on-shell S-production threshold, as can be seen also in Fig. 5. In the singlet model large areas of parameter space outside of the \(M_F\sim M_S\) region remain unconstrained. Conversely, for the doublet model a significantly larger part of the parameter space is already probed due to the larger 4L cross sections, see Fig. 5. We find that values of \(M_F\) below 800 GeV are excluded, consistent with the CMS 95% CL limit of 790 GeV. Still, departing from the \(M_F\sim M_S\) region we find areas in the doublet model parameter space that are in agreement with the 4L data.

We choose three allowed benchmark points per model to illustrate the analysis strategy described in the following sections. These points are marked as yellow circles in Fig. 6. For the singlet model, one of the benchmarks features VLL masses as low as 300 GeV and \(M_S=800\) GeV, while another one presents \(M_F = 600\) GeV for the same scalar mass. The third benchmark, with \(M_F,M_S = 800,500\) GeV, displays the opposite mass hierarchy, allowing on-shell decays of the VLLs through the \(S_{ij}\). In the doublet model values below \(M_F=800\) GeV are excluded regardless of \(M_S\). We have chosen a benchmark which saturates this bound, with \(M_F,M_S = 800,1200\) GeV. The two remaining benchmarks present the inverse mass hierarchy, allowing for on-shell \(\psi \rightarrow S\ell \) decays. The chosen parameters are \(M_F,M_S = 850,500\) GeV and \(M_F,M_S = 1000,800\) GeV, which in both cases lie at the frontier of the probed parameter space (see again Fig. 6).

In Fig. 7 we show the \(L_T\) distributions for these benchmarks (long-dashed curves) to explicitly show that they pass 4L constraints, that is, are within the CMS plus \(1 \sigma \) range (hatched area).

6 Optimized observables and null tests

In this section we design novel observables which target specific flavor features of our models and can serve as null tests of the SM. These optimized observables consist of invariant mass distributions which aim at reconstructing the masses of the new scalars and the VLLs. The latter are reconstructed through their decays to electroweak bosons, h and S plus lepton; final states with neutrinos are mostly removed through cuts on the missing transverse momentum, see Table 2. Thanks to the large values of \(\kappa ^\prime \), VLL decays to \(S_{ij}\) plus charged lepton are dominant when the \(S_{ij}\) can be produced on-shell, see Fig. 3, but remain significant also for \(M_F\lesssim M_S\). Key modes to probe VLLs and scalars with flavor are the decays (10) and (11) of the negatively charged \(\psi _i\) into six leptons

$$\begin{aligned} \begin{aligned} \psi _i {\overline{\psi }}_i&\rightarrow \ell ^-_ j S_{ij}^* \ell ^+_k S_{ik} \rightarrow \ell ^-_j \ell ^+_j\ell ^-_i \ell ^+_k \ell ^-_k\ell ^+_i \quad \mathrm{(singlet),} \\ \psi _i^- \psi _i^+&\rightarrow \ell ^-_ j S_{ji} \ell ^+_k S_{ki}^* \rightarrow \ell ^-_j \ell ^+_j\ell ^-_i \ell ^+_k \ell ^-_k\ell ^+_i \quad \mathrm{(doublet)\,} \, , \nonumber \\ \end{aligned}\\ \end{aligned}$$
(14)

which enable to construct the \(S_{ij}\) out of two leptons with opposite charge and same or different flavor. Combining these two leptons with a third one carrying the same flavor and opposite charge as one of the leptons in the initial pair enables us to reconstruct the \(\psi _i\). We reconstruct masses of the Z- and Higgs-boson as well as the masses of the new scalars \(S_{ij}\) considering jets (for Z and h only) and charged leptons as the final states. These invariant masses computed from two final-state particles are referred to as \(m_{2\ell }\). The subset of invariant masses reconstructed from leptons with different flavor are called \(m_{2\ell }\_{\mathrm{diff}}\). Combining the reconstructed bosons with the remaining charged leptons gives the reconstructed masses of the VLLs, called \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff}}\).

For the \(m_{i}\_{\mathrm{diff}}\) observables, in order to reconstruct both the \(S_{ij}\) and the \(\psi _i\) out of leptons with different flavors (\(i\ne j\) and/or \(i\ne k\) in Eq. (14)) we require to first find two pairs of different-flavor leptons with the same invariant mass within a small mass window \(\varDelta M_S\) (see Table 2), assuming a narrow width for the \(S_{ij}\). This allows to search for S-mediated decays without any assumption on the scalar masses. When two candidates for the new scalars are found, the VLLs are reconstructed applying flavor-conservation conditions. The flavor and mass requirements sufficiently suppress SM background making the processes in Eq. (14) with \(i\ne j\) and/or \(i\ne k\) the ‘golden channels’ for our analysis.

In Sect. 6.1 we discuss the algorithm to construct the observables \(m_{2\ell }\) and \(m_{2\ell }\_{\mathrm{diff}}\), which could signal scalar resonances. In Sect. 6.2 we discuss how to obtain \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff}}\) distributions, which could signal VLLs. Projections for the full Run 2 data set are worked out in Sect. 6.3.

6.1 \(m_{2\ell }\) and \(m_{2\ell }\_{\mathrm{diff}}\)

For each event with at least four light leptons, we compute all possible sets of two dilepton invariant masses from leptons of opposite charge, where each lepton contributes to only one of the invariant masses in the pair. This step includes \(\tau \) leptons if present. If the event contains jets, we include all possible pairs of invariant masses where one of them is a dilepton invariant mass and the other is the dijet invariant mass. For each event, only one pair of invariant masses is added to the observable \(m_{2\ell }\). In order to be added, it must fulfill one of the following requirements:

(a):

Each invariant mass is equal either to \(M_Z\pm \varDelta M_Z\) or to \(M_H \pm \varDelta M_H\), according to the parameters in Table 2, and each dilepton pair contains two leptons of the same flavor. The states used to compute the masses are in this case \((\ell _i^+\ell _i^-)(\ell _j^+\ell _j^-)\), \((\tau ^+\tau ^-)(\ell _i^+\ell _i^-)\) or \((\ell _i^+\ell _i^-)(jj)\) with \(i,j=1,2\). This condition reconstructs Z and Higgs bosons.

(b):

The difference between both invariant masses is less than \(\varDelta M_S\) (see Table 2), while none of the other invariant mass pairs present a smaller difference, and each invariant mass is computed from same-flavored leptons. The states used to compute the masses are in this case \((\ell _i^+\ell _i^-)(\ell _j^+\ell _j^-)\) with \(i,j=1,2\). This condition reconstructs two scalars, \(S_{ii}\) and \(S_{jj}\).

(c):

Both invariant masses differ by less than \(\varDelta M_S\) (see Table 2), while none of the other invariant mass pairs present a smaller difference, and at least one of the invariant masses contains two leptons of different flavor. The states used to compute the masses are in this case \((\ell _i^+\ell _j^-)(\ell _k^+\ell _i^-)\) with \(i,j,k=1,2,3\) leading to a maximum of one \(\tau ^+\) and one \(\tau ^-\). This condition reconstructs two scalars, \(S_{ij}\) and \(S_{ik}\).

We check for these conditions in the above order (\(a\rightarrow b \rightarrow c\)) and stop when one of the requirements is fulfilled. We define the observable \(m_{2\ell }\_{\mathrm{diff}}\) as invariant mass pairs that only fulfill condition c), where two particles of approximately equal invariant mass are found and at least one of them is computed from different-flavor leptons. All SM contributions to this observable are purely statistical, and therefore any significant excess away from SM resonances is an indication of new physics which can be explained by the VLL models of Eq. (2) or (4).

Fig. 8
figure 8

Di- and trilepton invariant mass distributions \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\), \(m_{3\ell }\), and \(m_{3\ell }\_{\mathrm{diff}}\) (see Sect. 6 for details) for the singlet model for different benchmark masses of the VLLs and the BSM scalars at a luminosity of 150 fb\({}^{-1}\) and \(\sqrt{s}=13\) TeV. The coupling \(\kappa '\) is fixed according to Eq. (7)

Full size image
Fig. 9
figure 9

Di- and trilepton invariant mass distributions \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\), \(m_{3\ell }\), and \(m_{3\ell }\_{\mathrm{diff}}\) (see Sect. 6 for details) for the doublet model for different benchmark masses of the VLLs and the BSM scalars at a luminosity of 150 fb\({}^{-1}\) and \(\sqrt{s}=13\) TeV. The coupling \(\kappa '\) is fixed according to Eq. (7)

Full size image

6.2 \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff}}\)

The \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff}}\) observables are designed to reconstruct the invariant masses of the VLLs via their three-body decays. For each pair of two-particle invariant masses added to \(m_{2\ell }\), we look for the additional lepton which stems from the decay of each \(\psi \). We add to \(m_{3\ell }\) the pairs of three-particle invariant masses which fulfill one the following conditions:

(i):

For two-particle invariant masses which reconstruct to a Z or Higgs (condition a in the previous section) each two-particle invariant mass is paired with an additional lepton present in the final state. The resulting three-particle invariant masses are added to \(m_{3\ell }\) if their difference is smaller than \(\varDelta M_F\), and no other combination presents a smaller difference.

(ii):

For two-lepton invariant masses which reconstruct to \(S_{ii}\) and \(S_{jj}\) (condition b in the previous section) each two-lepton invariant mass is paired with an additional lepton present in the final state which has the same flavor of the two leptons in the two-lepton invariant mass. The resulting three-lepton invariant masses are added to \(m_{3\ell }\) if their difference is smaller than \(\varDelta M_F\), and no other combination presents a smaller difference.

(iii):

For two-lepton invariant masses which reconstruct to \(S_{ik}\) and \(S_{kj}\) (condition c in the previous section) if a two-lepton invariant mass contains two same-flavor leptons, it is paired with an additional lepton present in the final state which has the same flavor. If it contains two different-flavor leptons, it is paired with an additional lepton which has the same flavor but opposite charge of one of the leptons in the two-lepton invariant mass. For each event, we find at most one possible combination that fulfills this condition. The corresponding three-lepton invariant masses are added to \(m_{3\ell }\).

In the last two conditions, flavor requirements are designed to reflect flavor conservation in the decays of the \(S_{ij}\). We define the observable \(m_{3\ell }\_{\mathrm{diff}}\) as invariant mass pairs that only fulfill condition (iii). As it turns out, the selection of the third leptons via flavor rules allows to populate \(m_{3\ell }\_{\mathrm{diff}}\) even when the \(\psi \)’s do not have a narrow width, which happens when the \(\psi \) undergoes frequent on-shell decays to S, i.e., for \(M_F>M_S\) and \(\kappa ^\prime \) large. \(m_{3\ell }\_{\mathrm{diff}}\) is a clean null test of the SM.

Fig. 10
figure 10

As in Fig. 8 after detector simulation, see Sect. 6 for details

Full size image

6.3 Benchmark distributions for Run 2

We study the \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\), \(m_{3\ell }\), and \(m_{3\ell }\_{\mathrm{diff}}\) distributions for allowed benchmark values (yellow circles in Fig. 6) of the VLL mass \(M_F\) and the BSM scalar mass \(M_S\) for the full Run 2 data set. Results based on the algorithms described in Sects. 6.1 and 6.2 are shown in Figs. 8 and 9 for the singlet and doublet model, respectively. The distributions are computed from \(5\times 10^4\) generated events and rescaled to an integrated luminosity of \(150~\mathrm{fb}^{-1}\), therefore small statistical fluctuations are present in the plots.

We learn the following generic features:

  1. (i)

    The results are qualitatively similar for the singlet and doublet models, with a larger cross section for the latter yielding more populated distributions. The “diff”-observables (plots to the right) are cleaner than their non-“diff” variants (plots to the left), i.e., have more efficient SM background suppression. For instance, note how \(m_{2\ell }\_{\mathrm{diff}}\) reduces the SM background around the Z mass in comparison to \(m_{2\ell }\). The \(m_{3\ell }, \ m_{3\ell }\_{\mathrm{diff}}\) (lower plots) are cleaner than the \(m_{2\ell }, \ m_{2\ell }\_{\mathrm{diff}}\) (upper plots) spectra. The \(m_{3\ell }\_{\mathrm{diff}}\) observable is SM background free.

  2. (ii)

    Resonance peaks from S-decays appear in the \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\) spectra for the benchmarks with \(M_F > M_S\), that is, when on-shell production of the scalars takes place. Narrow resonance peaks from \(\psi \)-decays appear in the \(m_{3\ell }, \ m_{3\ell }\_{\mathrm{diff}}\) distributions for the other benchmarks, which present \(M_F < M_S\). The latter condition eliminates the rapid on-shell decays to the scalars through the large \(\kappa ^\prime \) Yukawa.

  3. (iii)

    Distributions can signal BSM physics also in the tails away from a narrow resonance peak, or if none is present; see for instance the black curves in Figs. 8 and 9. All benchmarks display an excess above the SM in all distributions, with the exception of light VLLs \(M_F=300\) GeV and heavy-ish scalars \(M_S=800\) GeV in the singlet model (blue curve), which are underneath the SM contributions in \(m_{2\ell }, \ m_{2\ell }\_{\mathrm{diff}}\) spectra, but do show up in the \(m_{3\ell }, \ m_{3\ell }\_{\mathrm{diff}}\) distributions.

Including the effects of hadronization and finite detector resolution we show in Figs. 10 and 11 the observables for the singlet and doublet benchmark scenarios, respectively, after showering the events and applying a fast detector simulation. As expected, we find that peaks become broader and event rates drop. In the \(m_{2\ell }, \ m_{2\ell }\_{\mathrm{diff}}\) distributions the number of events in the peaks is reduced by roughly one order of magnitude, leading to \({\mathcal {O}}(1)\) events in the peaks for on-shell S production in the case \(M_F=800\) GeV in the singlet and \(M_F=850\), 1000 GeV in the doublet model. In the case of the \(m_{3\ell }, \ m_{3\ell }\_{\mathrm{diff}}\) observables, only in the singlet model and for small VLL masses \(M_F=300\) GeV (blue curves) we find \({\mathcal {O}}(1)\) events in the peaks in the \(m_{3\ell }\) distributions, while in all other scenarios and the \(m_{3\ell }\_{\mathrm{diff}}\) distributions the number of signal events is below one. Scaling factors comparing the number of events in the peak bins before and after detector simulation are given in Table 3 in Appendix A, where we also discuss in more detail the effects of the detector simulation.

Fig. 11
figure 11

As in Fig. 9 after detector simulation, see Sect. 6 for details

Full size image
Fig. 12
figure 12

As in Fig. 8 but for higher luminosity 3000 fb\({}^{-1}\) and \(\sqrt{s}=14\) TeV

Full size image
Fig. 13
figure 13

As in Fig. 9 but for higher luminosity 3000 fb\({}^{-1}\) and \(\sqrt{s}=14\) TeV

Full size image
Fig. 14
figure 14

Di- and trilepton invariant mass distributions \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\), \(m_{3\ell }\), and \(m_{3\ell }\_{\mathrm{diff}}\) for the singlet model with \(\kappa ^\prime =1\), for the full Run 2 luminosity 150 fb\({}^{-1}\) and \(\sqrt{s}=13\) TeV

Full size image

As we argued, the new observables have great sensitivity to flavorful BSM physics, and would benefit from higher luminosity. In the next section, we discuss perspectives for the HL-LHC.

7 Implications for the HL-LHC

As shown in Sect. 6.3 the discovery of a BSM sector consisting of VLLs and new scalars with a non-trivial flavor structure remains a challenging task at Run 2. Here we study the new observables \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\), \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff}}\) for the benchmarks scenarios of Sect. 6 at the HL-LHC, at higher luminosity \(3000~\mathrm{fb}^{-1}\), for upgraded detectors and \(\sqrt{s} = 14\) TeV [25]. For the detector simulation we employ DELPHES3 with the HL-LHC default card instead of the CMS default card. The distributions before detector simulation are shown in Figs. 12 and 13 and after detector simulation in Figs. 16 and 17 for the singlet and doublet models, respectively.

The HL-LHC setup enhances event rates relative to Run 2, in both models and benchmarks, both signal peaks and the SM background, according to \(\sim 3000/150=20\). Despite the different detector settings, and the increased center of mass energy, the corresponding distributions from Run 2 and the HL-LHC are very similar. For example, the singlet model 3000 fb\(^{-1}\) plots in Fig. 12 essentially look like scaled-up versions of the 150 fb\(^{-1}\) ones shown in Fig. 8. The scaling factors between no detector simulations and including them given in Table 3 remain also very similar between the two LHC settings with larger scaling in the HL-LHC scenario due to improved detector settings, see Appendix A for details. For example, the singlet model 3000 fb\(^{-1}\) plots with detector simulation in Fig. 16 essentially look like scaled-up versions of the 150 fb\(^{-1}\) ones shown in Fig. 10, and similarly for the doublet model. At the HL-LHC the new observables continue to feature great separation of BSM signals from the SM background, just with (more) events.

In the \(m_{2\ell }\_{\mathrm{diff}}\) spectra the bins with \(m_{2\ell }\_{\mathrm{diff}} \gtrsim 500\) GeV allow to search for both on-shell and off-shell S-production. For the former, we find \({\mathcal {O}}(10^3)\) events (\({\mathcal {O}}(20)\) after detector simulation) in the peaks of the \(m_{2\ell }\_{\mathrm{diff}}\) distribution for both the singlet and doublet model.

Fig. 15
figure 15

Trilepton invariant mass distributions \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff}}\) for the singlet model with \(\kappa ^\prime =1\) after detector simulation, for the full Run 2 luminosity 150 fb\({}^{-1}\) and \(\sqrt{s}=13\) TeV

Full size image

The \( m_{3\ell }\_{\mathrm{diff}}\) distributions have in both models \({\mathcal {O}}(10{-}10^2)\) events per bin (\({\mathcal {O}}(1{-}10)\) events per bin after detector simulation) with the exception of the doublet benchmark with light VLLs \(M_F \ll M_S\) (red curve in doublet model). Here, the \( m_{3\ell }\) distribution turns out to be powerful and produces up to \({\mathcal {O}}(10{-}10^2)\) events per bin after detector simulation. The \( m_{3\ell }\) distribution enhances also the peak in the singlet model with light VLLs and hierarchical spectrum \(M_F=300\) GeV, \(M_S=800\) GeV (blue curves in singlet model) up to this level.

We conclude that at the HL-LHC the (\(M_F\), \(M_S\))-parameter space consistent with the \(g-2\) anomalies can be probed and mass hierarchies extracted.

8 Beyond \(g-2\)

In the previous sections we have studied the VLL models in Eqs. (2) and (4) in the parameter space where the coupling \(\kappa ^\prime \) is fixed by the BSM masses \(M_F,M_S\) (7). In this section, we analyze the model space beyond the \(g-2\) constraint, entertaining the possibility of a shift in \(\varDelta a_\mu \) due to improved data and theory.

In general, the coupling \(\kappa \) remains limited in magnitude from above by Z decays, inducing small but relevant effects in fermion mixing. On the other hand, \(\kappa ^\prime \) is unconstrained by electroweak data. As it is already rather sizable in the benchmark (7), we investigate the implications of a reduced \(\kappa ^\prime \). The latter implies a suppression of \(\psi \) to S plus lepton decays. Since these modes are the dominant ones for \(M_F >M_S\), see Fig. 3, the width of the \(\psi \) in this region is proportional to \(\kappa ^{\prime 2}\). We expect therefore narrower resonances in \(M_F >M_S\) and a suppression of events in the region \(M_F <M_S\).

One may wonder what happens if \(\kappa ^\prime \) vanishes. In this case the models could still produce lepton flavor violation-like signals for \(y\ne 0\), with \(\psi \rightarrow S\ell \) happening at order \(y\theta \) and \(S \rightarrow \ell \ell ^{(\prime )}\) at order \(y \theta ^2 v_h/M_F\) times the lepton Yukawa with the Higgs. Due to the Z-constraints on the mixing angle \(\theta < {\mathcal {O}}(10^{-2})\) the “diff”-observables would be strongly suppressed up to some statistical noise. This outcome holds also for other UV-safe models with flavorful VLLs in representations of \(SU(2)_L \times U(1)_Y\) which do not allow for a \(\psi \)-S-lepton Yukawa coupling (“\(\kappa ^\prime =0\)”) [9]. Furthermore, if \(M_S\) were very heavy, in all models with mixed SM-BSM Yukawas the “diff”-observables would be SM-like. Note that the phenomenology of VLLs without any mixed SM-BSM Yukawas (“\(\kappa =\kappa ^\prime =0\)”) is markedly different and discussed for various exotic representations in [7]. In the following we focus on the singlet model, since the allowed parameter space of BSM masses is larger, see Fig. 6. For simplicity we use \(\kappa ^\prime =1\). We find that the benchmark \(M_F=800\) GeV, \(M_S=500\) GeV and \(\kappa ^\prime =1\) is excluded by CMS data, even though the \(g-2\) benchmark scenario with a larger coupling but the same BSM masses was found to fall within the allowed region (see Fig. 7). Nevertheless, we observe that larger VLL masses \(M_F\sim 900\) GeV are allowed for \(M_S\sim 500\) GeV and \(\kappa '=1\). Therefore, in this section we take \(M_F=900\) GeV, \(M_S=500\) GeV as one of our benchmarks. We study as well the \(\kappa '=1\) counterparts of the two remaining benchmarks considered in previous sections, which we find to be allowed.

The corresponding distributions of the new observables at the \(\sqrt{s}=13\) TeV LHC and the full Run 2 data set are shown in Fig. 14. We observe that for \(m_{2\ell }\) and \(m_{2\ell }\_{\mathrm{diff }}\) (see Fig. 14, upper row), the patterns are very similar to the \(g-2\) benchmarks shown in Fig. 8, with resonance peaks when the scalars can be produced on-shell (black curves). For \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff }}\) (see Fig. 14, lower row), however, reducing \(\kappa ^\prime \) leads to qualitatively different effects. As anticipated, these are mainly due to the VLLs’ narrower widths, which can be seen in all benchmarks. \(m_{3\ell }\) is the observable with more distinctive peaks regardless of the BSM mass hierarchy, with resonances above the SM background and reaching at least \({\mathcal {O}}(10)\) events in the peak bin for all benchmarks. For \(M_F>M_S\) (black curves) \(m_{3\ell }\_{\mathrm{diff }}\) is the optimal observable, since all SM background is suppressed and the number of events per bin barely decreases with respect to \(m_{3\ell }\). For \(M_F<M_S\), the \(m_{3\ell }\_{\mathrm{diff }}\) distributions are substantially depleted with respect to \(m_{3\ell }\), and higher luminosities would be beneficial.

In Fig. 15 we give the \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff }}\) distributions after hadronization and detector simulation The \(m_{2\ell }\) and \(m_{2\ell }\_{\mathrm{diff }}\) spectra are very similar to the ones in Fig. 10 and therefore not shown. The \(m_{3\ell }\) distributions show peaks for all three benchmarks, with \({\mathcal {O}}(1)\) events per bin for both \(M_F=300\) GeV (blue) and \(M_F=900\) GeV (black). For the latter, the \(m_{3\ell }\_{\mathrm{diff }}\) distribution allows for a null test, as the SM background is sufficiently suppressed while we find a peak with few events per bin in the \(M_F=900\) GeV distribution.

9 Summary

We investigated opportunities at the LHC and the HL-LHC to search for flavorful vector-like leptons \(\psi _i\) and new scalar singlets \(S_{ij}\). Such BSM sector (1) occurs in novel model building frameworks with favorable UV behavior [4,5,6], and particle physics phenomenology [7, 10].

We considered two explicit BSM models of this kind, featuring three generations of either \(SU(2)_L\) singlet or doublet VLLs, which can also accommodate present data of the muon and electron \(g-2\). Key ingredients for flavor phenomenology are the mixed SM-BSM Yukawa couplings, the flavor matrix structure of the BSM scalars, the identification of lepton and VLL flavor, and fermion mixing after electroweak symmetry breaking. Although all BSM interactions (2) and (4) are flavor-conserving, the decays of the VLLs through the \(S_{ij}\) and their subsequent decay \(S_{ij}\rightarrow \ell _i^+\ell _j^-\) lead to production of different-flavor lepton pairs, a signature we exploited to construct novel null tests of the SM: The dilepton invariant masses \(m_{2\ell }\) and \(m_{2\ell }\_{\mathrm{diff}}\), which permit to look for scalar resonances, and the three-lepton invariant masses \(m_{3\ell }\) and \(m_{3\ell }\_{\mathrm{diff}}\), which are designed to reconstruct VLL masses, as described in Sect. 6. The \(\_{\mathrm{diff}}\) distributions are populated exclusively by invariant masses which contain at least two leptons of different flavor and opposite charge, which results in a strong suppression of the SM background and targets models with a non-trivial flavor structure affecting the charged lepton sector. The background suppression is especially efficient for \(m_{3\ell }\_{\mathrm{diff}}\), which makes this an excellent null test.

For our study we implemented the models into UFO models using FeynRules. Predictions for observables including dominant SM background processes at pp-colliders are computed with MadGraph5_aMC@NLO together with PYTHIA8 and DELPHES3. We worked out constraints from a CMS search in final states with at least four light leptons (electrons, or muons) [24]. Results are summarized in Fig. 6, showing allowed regions (green and yellow circles) of VLL and scalar masses while accommodating \(g-2\) of the muon, (7). We find that, in general, regions around the \(M_S = M_F\) line are excluded by data due to the underlying enhancement of cross sections via on-shell S production. Lower limits for the VLL masses are around 300 GeV in the singlet model and 800 GeV in the doublet model. As such, our findings offer new constraints for the Planck safe models put forward in [10].

Fig. 16
figure 16

Di- and trilepton invariant mass distributions \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\), \(m_{3\ell }\), and \(m_{3\ell }\_{\mathrm{diff}}\) after detector simulation for the singlet model. The observables are shown for different benchmarks of the VLLs and BSM scalar masses at a luminosity of 3000 fb\({}^{-1}\) and \(\sqrt{s}=14\) TeV. The coupling \(\kappa '\) is fixed according to Eq. (6)

Full size image
Fig. 17
figure 17

Di- and trilepton invariant mass distributions \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\), \(m_{3\ell }\), and \(m_{3\ell }\_{\mathrm{diff}}\) after detector simulation for the doublet model. The observables are shown for different benchmarks of VLLs and BSM scalar masses at a luminosity of 3000 fb\({}^{-1}\) and \(\sqrt{s}=14\) TeV. The coupling \(\kappa '\) is fixed according to Eq. (6)

Full size image
Table 3 Scaling factors \(f = N_{\mathrm{peak, det}}/N_{\mathrm{peak}}\) for the observables of Sect. 6 for different benchmarks, with \(N_{\mathrm{peak, det}}\) \((N_{\mathrm{peak}})\) denoting the number of events at the peaks per bin after (before) detector simulation for \(\sqrt{s} = 13\) TeV and a luminosity of 150 fb\({}^{-1}\) and in parentheses for \(\sqrt{s} = 14\) TeV and 3000 fb\({}^{-1}\). We marked with * (**) the cases where the peaks fall under SM background (resonances are broad)
Full size table

Predictions for the new observables \(m_{2\ell }\), \(m_{2\ell }\_{\mathrm{diff}}\), \(m_{3\ell }\), and \(m_{3\ell }\_{\mathrm{diff}}\) after detector simulation are shown for several allowed benchmarks in Figs. 10 and 11 for the singlet and doublet model, respectively, for the full Run 2 data set with 150 fb\({}^{-1}\). The distributions exhibit a highly discriminating power on the BSM mass hierarchy, but suffer from marginal event rates and therefore would extremely benefit from higher luminosity. At the HL-LHC, for \(\sqrt{s}=14\) TeV and a luminosity of 3000 \(\mathrm{fb}^{-1}\), we obtain \({\mathcal {O}}(10^2)\) events after detector simulation in some bins, see Figs. 16 and 17 for the singlet and doublet model, respectively. Hence, these new, optimized observables are very promising for higher luminosity runs at the LHC, to discover and discern hierarchies in flavorful models with multi-lepton final states.

Studying more general versions of our models in Sect. 8 we reduce \(\kappa ^\prime \), the key Yukawa for filling the “diff”-distributions. Results are shown in Fig. 15 for the invariant mass distributions after detector simulation. We again observe striking BSM signatures with diagnosing power. Let us also mention that the other colorless models with effectively \(\kappa ^\prime =0\) put forward as asymptotically safe extensions of the SM [9] are not contributing significantly to the “diff”-observables, but could be probed using \(m_{2 \ell }\) and \(m_{3 \ell }\) or conventional VLL search strategies.