Doubly Charged Higgs Boson Production at Hadron Colliders II: A Zee-Babu Case Study

Motivated by searches for so-called leptonic scalars at the LHC and the recent measurement of the $W$ boson's mass at the Tevatron, we revisit the phenomenology of the Zee-Babu model for neutrino masses and the ability to differentiate it from the Type II Seesaw model at the LHC. We conclude that this task is much more difficult than previously believed. All inputs equal in the two scenarios, we find that total and differential rates for producing pairs of doubly and singly charged scalars are identical in shape and only differ in normalization. The normalization is given by the ratio of hadronic cross sections and can be unity. Differences in cross sections are small and can be hidden by unknown branching rates. This holds for Drell-Yan, $\gamma\gamma$ fusion, and $gg$ fusion, as well as observables at LO and NLO in QCD. This likeness allows us to reinterpret Run II limits on the Type II Seesaw and estimate projections for the HL-LHC. Using updated neutrino oscillation data, we also find that some collider observables, e.g., lepton flavor-violating branching ratios, are now sufficiently precise to provide a path forward. Other means of discrimination are also discussed. As a byproduct of this work, we report the availability of new Universal \texttt{FeynRules} Object libraries, the \texttt{SM\_ZeeBabu} UFO, that enable fully differential simulations up to NLO+LL(PS) with tool chains employing \texttt{MadGraph5\_aMC@NLO}.


Introduction
It is a fact of life that neutrino oscillation data support large mixing angles and at least two neutrinos having tiny masses [1,2], whereas the Standard Model of particle physics (SM) postulates that all three neutrinos are massless Weyl fermions. As the SM is an otherwise successful description of data across many scales, one is inclined to extend the model in order to reconcile this discrepancy and reproduce the well-established [3,4] Pontecorvo-Maki-Nakagawa-Sakata (PMNS) paradigm [5][6][7] from more fundamental principles. What is less clear, is which, if any, of the many neutrino mass models throughout the literature are at least partially correct.
Among the most studied neutrino mass models are those that conjecture the existence of new fermions, such as right-handed (RH) neutrinos (ν R ) and vector-like leptons. While motivated, the most minimal [8] incarnations of these, e.g., the Types I [9][10][11][12][13][14][15] and III [16] Seesaws models, introduce coupling and/or mass hierarchies. Subsequently, additional states and interactions are typically needed to soften these hierarchies or render them experimentally testable. In light of this and of the limited guidance provided by data and theory, it is important to take a broad, complementary approach to neutrino mass models, including exploring scenarios without ν R .
Along these lines, the Type II Seesaw model [15,[17][18][19][20] and Zee-Babu model [21][22][23] are typical scenarios that can reproduce oscillation data without invoking ν R . Instead, these models postulate the existence of so-called leptonic scalars, i.e., scalars that carry nonzero lepton number (LN), that also couple directly to electroweak (EW) gauge bosons. In the Type II case, the existence of a scalar SU(2) L triplet∆ is hypothesized and left-handed (LH) Majorana neutrino masses are sourced at tree level from the vacuum expectation value (vev) of∆. In the Zee-Babu case, the existence of two scalar SU(2) L singlets k, h, which carry hypercharge, are hypothesized and LH Majorana neutrino masses are generated radiatively at two loops. In both cases, neutrino masses are proportional to a parameter µ L that signals the scale at which LN is broken.
In this work, we revisit the phenomenology of the Zee-Babu model and the extent to which it can be distinguished from Type II Seesaw at the hadron colliders. This study is motivated, in part, by recent the measurement of the W boson's mass at the CDF experiment [41]. With high significance, the collaboration reports a mass larger than predicted by precision EW data, but also one that fits naturally in the Type II Seesaw [42][43][44][45][46][47]. Subsequently, the measurement is reigniting interest in searches for doubly and singly charged scalars at the LHC [48][49][50][51]. If a discovery of exotically charged scalars follows soon at the LHC, it will be paramount to exclude one or both models [52][53][54]. The aim of the work is to provide some guidance in this direction.
Comparative studies of doubly charged scalars in the Zee-Babu model and other scalar extensions of the SM have been conducted in the past [35,37,38,40,52,[55][56][57][58][59]. These, however, are largely restricted to comparing inclusive cross sections via Drell-Yan (DY) at leading order (LO) in quantum chromodynamics (QCD), or brute-force recasting using simulated events at LO with parton shower-matching at leading logarithmic accuracy (LO+LL(PS)). In this work, we report an interesting observation. Namely, for the DY, gluon fusion, and photon fusion mechanisms, and for all inputs equal, hadronic cross sections (Sec. 5.1) and differential (Sec. 5.2) distributions of scalars in the Zee-Babu and Type II models differ at most by a uniform scaling factor. That is to say, for fixed masses, production channel, etc., the shapes of kinematical distributions in the two scenarios are the same. Therefore, fiducial cross sections predicted by one theory can be obtained for the other by a naïve re-scaling. This holds for observables at LO+LL(PS) and next-to-leading order in QCD with PS-matching (NLO+LL(PS)). For some processes, the scale factor, which is given by the ratio of hadronic cross sections, is precisely unity at LO in the EW theory. For others, it is nearly independent of scalar mass but inherits a weak sensitivity via the dynamics of parton density functions (PDFs). In all cases, differences in cross sections between the two models are sufficiently small that they can be hidden by unknown branching rates.
This observation means that experimentally distinguishing the Zee-Babu and Type II Seesaw models will be more difficult than previously believed. However, our findings also show that LHC searches for charged scalars decaying directly to leptons in the Type II Seesaw can automatically be reinterpreted in the context of the Zee-Babu model. Using Ref. [51], we estimate (Sec. 5.3) that k ∓∓ masses as large as m k = 890 GeV and decay rates as small as BR(k ∓∓ → ∓ ∓ ) = 16% for ∈ {e, µ} are excluded by the ATLAS experiment with L ≈ 139 fb −1 of data at √ s = 13 TeV. Assuming constant analysis and detector performance, we project this can reach m k = 1110 GeV and BR(k ∓∓ → ∓ ∓ ) = 8% with L = 3 ab −1 at √ s = 13 TeV. Furthermore, with updated neutrino oscillation data, we find (Sec. 5.4) that predictions for flavor-violating h ∓ → ∓ ν decays in the Zee-Babu are now sufficiently precise to have discriminating power. As a byproduct of this work, we report the availability of new Universal FeynRules Object libraries [60], the SM ZeeBabu UFO, that enable simulations up to NLO+LL(PS) with Monte Carlo tool chains involving MadGraph5 aMC@NLO [61,62]. We highlight that previously published UFOs capable of simulating Zee-Babu scalars [52] can only achieve simulations up to LO+LL(PS).
This report continues according to the following: In Sec. 2, we describe the theoretical framework in which we work. In Sec. 3, we summarize our computational setup and the tuning of our Monte Carlo (MC) tool chain. In Sec. 4, we broadly revisit and update the phenomenology of Zee-Babu model, including the first NLO in QCD predictions for Zee-Babu scalars at the LHC. We present our main results in Sec. 5, where we discuss similarities and differences of Zee-Babu and Type II scalars at hadron colliders. We summarize and conclude in Sec. 6.

Theoretical framework: the Zee-Babu model
The Zee-Babu model [21][22][23] extends the SM by two complex scalars, k and h, with the quantum number assignments (1, 1, −2) and (1, 1, −1) under the SM gauge group G SM =SU(3) c ⊗SU(2) L ⊗U(1) Y . Neither carries color or weak isospin but both are charged under weak hypercharge. k and h are assigned lepton number L = +2, which is normalized such that SM leptons carry L = +1. In terms of the SM Lagrangian (L SM ), the Lagrangian of the Zee-Babu model (L ZB ) is L ZB = L SM + L Kin. + L Yuk. + L ZB scalar + δL ν . (2.1) The kinetic part of the Lagrangian for k and h is given by the following covariant derivatives Here, the weak hypercharge operator is normalized such that the electromagnetic charge operator isQ =T 3 L +Ŷ , and Y k = −2 (Y h = −1). The weak hypercharge coupling is denoted by g Y ≈ 0.36. As neither k nor h mix with SM states, the mass eigenstates, denoted by k −− and h − , are aligned with their gauge states and carry the electric charges Q k = −2 and Q h = −1, respectively. (In the following, we often omit the charge symbols ∓ when referring to mass eigenstates.) After EW symmetry breaking (EWSB), the hypercharge field B µ mixes with the weak isospin field W 3 µ and can be decomposed in terms of the usual mass eigenstates B µ = cos θ W A µ − sin θ W Z µ , where θ W is the weak mixing angle. In the mass basis, Eq. (2.2) becomes As usual, e = g Y cos θ is the electromagnetic coupling. The list (S, S † , Q s ) is shorthand for the mass eigenstates and charge (k −− , k ++ , Q k ) and (h − , h + , Q h ). In anticipation of Sec. 5.1, we highlight that h and k couple to the Z boson only through B − W 3 mixing since they are SU(2) L singlets. And since the weak mixing angle (modulo running) is only about θ W ≈ 29 • , h and k inherently couple more to the photon than to the Z.
Here, (L i ) T = (ν i L , e i L ) is the SM LH lepton doublet with generation index i = 1, . . . , 3;L i ≡ iσ 2 (L c ) i is the usual rotation in SU(2) L space but of L's charge conjugate; and e i R is the SM RH charged lepton. The Yukawa couplings to h are given by f ij , a 3 × 3, complex matrix that is anti-symmetric, i.e., f ij = −f ji . The Yukawa couplings to k are given by g ij , a 3 × 3, complex matrix symmetric, i.e., g ij = g ji . After EWSB, the chiral states e i L/R can be rotated trivially into their flavor/mass eigenstates = e, µ, τ . At this point neutrinos are still massless, meaning that their gauge and flavor states are also aligned. Formally, this leads to the redefinition of g = R * i g ij R j and f = R * i f ij R j , where R and R are the identity matrix. The Yukawa couplings g and f induce LFV in decays and transition of τ and µ. While this is constrained by experimental searches for LFV, k and h with masses below 1 TeV are still allowed [28,35].
The scalar potential of k and h, including couplings to the SM Higgs doublet Φ, is given by The states H 0 , G ± , and G 0 are the usual SM Higgs and Goldstone bosons, with After EWSB, one has in the mass basis The physical masses of k and h are, respectively, Demands for a first-order EW phase transition favors lighter masses, with m k , m h < 300 GeV [63]. A few comments: (i) In this work, we adopt the conventional assignments of LN wherein leptons carry L = +1, h and k carry L = +2, and antiparticle states carry L < 0. This implies that the three-point vertex h−h−k, which is proportional to dimensionful parameter µ L , violates LN explicitly by ∆L = ±2 units. We identify µ L as the scale of LNV. In the µ L → 0 limit, LN is conserved in the Zee-Babu model; conversely, in the limit where µ L is fixed but either k or h is infinitely heavy, i.e., the decoupling limit [64], the h − h − k vertex vanishes and leads to LN conservation. (ii) The assignment also implies that the Yukawa interactions in Eq. (2.4) conserve LN. However, in the absence of µ L , g, f , or the Yukawa couplings y between the SM Higgs and charged leptons, LN can be redefined such that it is conserved [35]. (iii) In L ZB scalar , the coupling normalizations follow Refs. [35]. In this convention, the h − h − k Feynman rule, Γ h−h−k = −i(2!)µ L , carries a factor of (2!) that would otherwise cancel in other normalizations. (iv) In this model, neutrinos are massless at tree level after EWSB, i.e., δL ν = 0. Unlike the Types I-III Seesaws, they are generated radiatively. Discussion of this is postponed to Sec. 4.1.

Computational setup and Monte Carlo tuning
In order to assist reproducing our results, we now document our Monte Carlo tool chain (Sec. 3.1), our SM inputs (Sec. 3.2), and our benchmark Zee-Babu inputs (Sec. 3.3).

Monte Carlo tool chain
To study the Zee-Babu model numerically, we transcribe the tree-level Lagrangian with Goldstone boson couplings in Eq. (2.1) into FeynRules v2.3.36 [65,66]. For the SM Lagrangian, we use the implementation available in FeynRules, the file sm.fr v1.4.7. We phenomenologically parameterize the Lagrangian and set δL ν = 0. QCD ultraviolet and R 2 counter terms up to O(α s ) are computed using NLOCT v1.02 [67] and FeynArts 3.11 [68]. Feynman rules up to one loop in α s are then generated and packaged into a series of Universal FeynRules Output (UFO) libraries [60] that we collectively label the SM ZeeBabu UFO libraries 1 . In this work, we use the SM ZeeBabu NLO UFO, which enables the computation of tree-induced processes up to NLO+LL(PS) in QCD and QCD loop-induced processes up to LO+LL(PS). We have checked that our implementation of the model agrees with partonic expressions for k and h pair production at LO [55], as well as with hadronic-level rates for k pair production at LO [52][53][54].

Standard Model inputs
For numerical results, we assume n f = 5 massless quarks and the following SM inputs [80]: We otherwise assume charged lepton are massless 3 . We approximate the Cabbibo-Kobayashi-Maskawa matrix by the identity matrix. For hadronic cross sections we use the MMHT 2015 QED NLO (lhaid=26000) and next-to-next-to-leading order (NNLO) (lhaid=26300) PDF sets [81]. Both PDF sets employ to the LUXqed formalism to determine the photon PDF [82,83] and use α s (M Z ) ≈ 0.1180. PDFs and α s (µ r ) are evolved using LHAPDF v6.3.0 [84]. PDF uncertainties are extracted using eigenvector sets [85] as implemented in LHAPDF. For all DY processes we use the NLO PDF set; for all non-DY processes we use the NNLO PDF set.
For DY calculations at LO and NLO, we set the central collinear factorization (µ f ) and renormalization (µ r ) scales to be half the sum of transverse energies of final state particles: For all other calculations, we set the two scales equal to the scale of hard scattering process: The 9-point scale uncertainty is obtained by varying ζ over the discrete range ζ ∈ {0.5, 1.0, 2.0}. Finally, the shower scale µ s is set according to its default prescription [62]. To steer the shower, we use the MSTW 2008 LO PDF set (lhaid=21000) [85] and the ATLAS A14 central tune (Tune:pp = 20) [86]. We do not estimate the uncertainty associated with PS modeling.

Zee-Babu inputs
As neutrinos are effectively massless on momentum scales observed at the LHC and to minimize potential theoretical biases, we take a phenomenological approach and neglect neutrino masses for collider computations. In practice, this means that the relationships in Eqs. (4.6), (4.9), and (4.10) are not imposed. This allows us to vary nonzero f and g freely and independently. Unless specified, we assume the following model benchmark inputs

Neutrino masses
In the Zee-Babu model, there are no ν R and therefore no Dirac neutrino masses. Likewise, k and h cannot contract with L and Φ so as to generate LH Majorana masses at tree level. Instead, the µ L term in L ZB scalar induces LH Majorana masses at two loops. In the flavor basis and with flavor indices i, j, a, b ∈ {e, µ, τ }, neutrino masses are described by Lagrangian [22,23] The integral factor I ab (r) is approximately given by the expression [35,87] Here, r = (m 2 k /m 2 h ) and M max = max(m k , m h ). Rotating neutrinos from their flavor states (ν ) into their mass states (ν m ) via the PMNS matrix [5][6][7], i.e., allows one to diagonalize the mass matrix. The result is For Yukawa couplings f ∼ g ∼ O(0.1), and charged lepton masses (m a m b ) ∼ O(0.1) GeV 2 , neutrino masses that are naturally O(1) eV can be obtained from Zee-Babu mass scales µ L , M max that are O(100) GeV. Given presently available oscillation data [4], the above expression impose meaningful constraints on f ij and g ij . However, exploring the rich complementarity of oscillation data, low-energy flavor data, and high-energy collider data for the Zee-Babu model is outside the scope of this work. Such studies have been conducted in Refs. [28,35,37,39]. Nevertheless, we comment on a nontrivial correlation between oscillation data and the Yukawa couplings.
The antisymmetric nature of f implies a zero determinant: and subsequently that det(M mass ν ) = 0. This forces at least one neutrino to be massless. Thus, the mass spectrum is fixed by the measured atmospheric and solar mass splittings, up to the mass-ordering ambiguity. Moreover, one can build eigenvector equations that relate elements of f and U PMNS . For the normal ordering (NO) of neutrino masses, one has [35] f eτ f µτ = tan θ 12 cos θ 23 cos θ 13 + tan θ 13 sin θ 23 e −iδ , (4.9a) f eµ f µτ = tan θ 12 cos θ 23 cos θ 13 − tan θ 13 sin θ 23 e −iδ , (4.9b) while for the inverse order (IO), one has Relationships for g also exist but are more complicated. For both mass orderings, consistency between oscillation data and the hierarchy of charged lepton masses leads to scaling [35]:

Decay channels of k ∓∓ and h ∓
We now comment on the leading and sub-leading decays of k ∓∓ and h ∓ . While formula for two-body partial widths (Γ) have been documented before, not all of the following properties have been reported. For an n f -body state f , the i → f partial width is given by the formula Here, M is the i → f matrix element; the summation in the first line is over discrete degrees of freedom (dof); e.g., helicity, S i is the spin-averaging multiplicity; and dP S n f is the phase space integration measure. The total width (Γ Tot. i ) of i and the i → f branching rate (BR) are then: Starting with the state k ∓∓ , the leading two-and three-body decay channels include The first proceeds by the symmetric Yukawa coupling g and conserves LN since k ∓∓ carries L = ±2. The second proceeds through µ L and violates LN since h ∓ also carries L = ±2. The last three are radiative corrections to the first two and are coupling or phase-spaced suppressed. The k ∓∓ → ± ± partial width with full lepton-mass dependence is given by The Kronecker δ accounts for 1/(2!) symmetry factor for identical particles in the final state. Assuming m k > 2m h , the k ∓∓ → h ∓ h ∓ partial width is given by Unusually, this decay is inversely proportional to the mass of k; normally, partial widths grow as a positive power of a parent particle's mass. Subsequently, the k ∓∓ → h ∓ h ∓ partial width can  h (column 4-5), assuming g , f = 1, as well as the normalized partial widths for k ∓∓ → ± ± (column 6), k ∓∓ → h ∓ h ∓ (column 7), and h ∓ → ± ν (column 8). In parentheses are the branching rates. be suppressed if the scale of LNV is much smaller than m k . At the same time, the branching rate can be competitive, or even dominant, if the couplings g and f are sufficiently small.
For the benchmark inputs in Eqs. (3.2) and (3.4), one finds the following branching rates: For our values of m k and m h , the two-body k −− → h − h − decay is kinematically forbidden. However, the largeness of µ L enhances the three-body, LN-violating k −− → − ν h − decay to the per mil level. The hierarchy displayed by three-body decays k −− → − ν W − reflects the fact the rates are proportional to charged lepton masses. More specifically, the k −− → − ν W − decay proceeds through the intermediate step k −− → − − * → − ν W − , which is mediated by the coupling of k −− to a RH lepton − R and the coupling of W − to LH leptons. This implies that − * must propagate in its RH helicity state, and hence that the amplitude scales with its mass. We caution that these rates are only illustrative. They assume that all nonzero g and f are unity; realistic values must be more varied to satisfy oscillation and flavor data [28].
Assuming that m k < m h , the leading two-and three-body decay channels are In analogy to k ∓∓ , the first channel proceeds through the antisymmetric Yukawa coupling f . The last three can be classified as being radiative corrections to the first channel and therefore are suppressed. The final proceeds through the LN-violating k − h − h vertex. The h ± → ± ν partial width with charged lepton mass dependence is given by For the benchmark inputs listed in Eqs. (3.2) and (3.4), one finds the following branching rates These channels proceed through LH chiral states and so little dependence on m is observed. For representative masses m k , m h and coupling µ L (columns 1-3), we summarized in Table 2 the total widths Γ Tot.
h (column 4-5), assuming g , f = 1. We also summarize the (normalized) partial widths for the processes k ∓∓ → ± ± (column 6), k ∓∓ → h ∓ h ∓ (column 7), and h ∓ → ± ν (column 8). In parentheses are the corresponding branching rates assuming the benchmarks listed in Eq. (3.4). For these inputs, we find that the total widths of k and h span approximately Γ Tot.  We postpone further exploitation of correlations among k and h decays to Sec. 5.4.

Constraints from partial wave unitarity
Before presenting predictions for the LHC, we consider constraints from partial-wave unitarity. For some Seesaw scenarios, the J = 0 partial wave in vector boson scattering (VBS) is known to be interesting at high energies [89][90][91][92]. For the Zee-Babu model, we focus on the channels As depicted in Fig. 1, the Born-level matrix element for ZZ → k ++ k −− is facilitated by a 4-point term (M 4 ), as well as t-(M t ), u-(M u ), and s-channel (M H ) terms. The first thee diagrams are determined entirely by EW couplings while the last is set by the Similarly, the W + W − channels are mediated each by the three diagrams (not shown) with s-channel γ/Z/H 0 appearing as intermediaries. The γ/Z diagrams are determined entirely by EW couplings while the third is set by λ kZBh and λ hZBh for k −− k ++ and h − h + , respectively.
Notably, the helicity amplitude for Z 0 Z 0 → k −− k ++ , where Z 0 is a longitudinally polarized Z boson, undergoes strong cancellations. Without any approximations, one finds the scaling (powers of s, m 2 k , and M 2 Z ) (higher powers of s, m 2 k , and M 2 Z ) . This means that in the high-energy limit, where s M 2 Z , m 2 k , the pure gauge contribution vanishes. Taking the (M 2 Z /s), (m 2 k /s) → 0 limit before summing diagrams leads to separately divergent terms. Some of this behavior can be attributed to the structure of longitudinal polarization vectors, which scale like ε µ (q, λ = 0) ∼ q µ /M Z + O(M Z /q 0 ). Intuitively, scalars in the Zee-Babu model carry only hypercharge, not weak isospin. Therefore, in the unbroken phase, they should decouple from the weak sector. The same behavior is found in the other channels.
Summing over all diagrams, followed by taking the high-energy limit leads to the simple expressions: For completeness, the partonic cross sections in this kinematic limit simplify to the expression where M V V is the invariant mass of the (V 0 V 0 ) system, and {λ} is either λ kZBh or λ ZBh . (No spin averaging is needed as V 0 V 0 are polarized.) Following the procedure of Ref. [93], we checked that our implementation of the Zee-Babu model (see Sec. 3.1) reproduces this cross section. The J = 0 partial-wave amplitudes are obtained from Eq. (4.24) using where V 0 ∈ {W 0 , Z 0 } and S ∈ {k, h}. This results in the following partial-wave amplitudes: The perturbative condition of While these bounds are relatively weak, more aggressive restrictions on |a J | translate into more aggressive limits on λ kZBh , λ hZBh . Further considerations from VBS is left to future work.

k ++ k −− and h + h − pairs at the LHC
We now turn to the production Zee-Babu scalars at the √ s = 13 TeV LHC. We focus on k −− k ++ and h − h + pair production through a variety of processes that are depicted at the Born level in Fig. 2. As an outlook, we also consider a hypothetical Very Large Hadron Collider (VLHC) at √ s = 100 TeV. Our results are summarized in Fig. 3, where we plot as a function of scalar mass (m k = m h ) the total inclusive cross section (σ) for these processes. For some channels these are the first LHC predictions that have been made in the context of the Zee-Babu model. In high-p T hadron collisions, the inclusive production rate of final-state F is given by [94][95][96] σ(pp → F + anything) = i,j∈{q,q,g,γ} whereσ ij→F is the partonic ij → F scattering cross section as obtained from the formulâ Here, M ij→F is the partonic matrix element calculable using perturbative methods and the Feynman rules of Sec. 2; N c and S are, respectively, the color and spin multiplicities of i and j; the summation is over all discrete dof / multiplicities; and Q 2 = (p i + p j ) 2 > M 2 (F) is the (squared) hard scattering scale. Q must exceed the invariant mass of F in order for the process to proceed. The phase space volume element for an n f -body final state is defined in Eq. (4.12b). In Eq. (4.29), the f are the collinear PDFs that represent the likelihood of finding partons i, j ∈ {q, q, g, γ}, for light quark species q, in proton p carrying particular longitudinal momentum fractions. ∆ ij describes the likelihood of soft radiation emitted in the ij → F scattering process. The symbol ⊗ denotes the convolution of these probabilities. Cross sections throughout this section assume the SM inputs of Sec. 3.2 and the Zee-Babu inputs of Eq. (3.4). For select (gluon fusion) computations, the Zee-Babu inputs of Eq. (4.42) are used and will be discussed below.
Drell-Yan: We begin with pair production via the Drell-Yan (DY) mechanism, i.e., quarkantiquark annihilation, which at the Born level is depicted in Fig. 2(a) and given by (4.31) To simulate this at NLO with the SM ZeeBabu libraries and mg5amc, we use the commands 4 : 4 See also Ref. [62] for instructions on operating mg5amc. For both k −− k ++ (black) and h − h + (teal) production, we show in Fig. 3(a) the inclusive production cross sections for √ s = 13 TeV at NLO in QCD and with residual scale uncertainties (band thickness). While predictions at NLO in QCD for Type II scalars have been available for some time [97,98], this is the first for the Zee-Babu model. Over the mass range m k , m h = 50 GeV − 1400 GeV, we find that the cross sections for the two processes span approximately with residual scale uncertainties spanning about δσ DY (NLO) 13 TeV ≈ ±2% − ±5% for both channels. An interesting observation is that ratio of the two DY rates is constant and equal to four. This can be attributed to the couplings of k and h to γ and Z. As shown in the Lagrangian of Eq. (2.3), the three-point S − S † − V vertices are proportional to the electric charge of the scalar but are otherwise the same for k and h. Hence, the rates differ by the square of the charges: Further investigation of the matrix element shows that in the large m k , m h limit, i.e., where O(M 2 Z /Q 2 ) terms can be neglected, a strong destructive cancellation occurs between the γ and Z diagrams. For both the uu and dd parton channels, the γ −Z interference term is negative and slightly larger in magnitude than the Z channel. (In the dd channel, the γ, Z, and interference terms are actually all comparable in size.) In essence, the DY channel is driven by the γ diagram.
To quantify the size of QCD corrections, we show the NLO in QCD K-factor (K NLO ) in the lower panel of Fig. 2(a) as a function of mass. For both DY channels, K NLO spans roughly These numbers are comparable to exotically charged scalar production in the Type II Seesaw [97,98] when adjusted for PDF and scale choices. In the lower panel, we also show the residual scale uncertainty at NLO (darker inner bands) and PDF uncertainty (lighter outer bands). PDF uncertainties roughly span δσ ≈ ±2% − ±8%. The curves and bands for k −− k ++ and h − h + pair production overlap almost perfectly. This follows from the DY rates for k and h differing by a constant at both LO and NLO, which cancel when taking the respective ratios.
Zee-Babu scalar mass [TeV] Photon Fusion: Next we consider inclusive pair production from γγ fusion (AF), given by and shown diagrammatically in Fig. 2(c). This process can be simulated at LO using the syntax generate a a > kk kk QED=2 QCD=0 output DirName3 generate a a > hh hh QED=2 QCD=0 output DirName4 Over the mass range investigated, the cross sections for the two processes approximately span with scale uncertainties reaching δσ AF (LO) 13 TeV ≈ ±20% (±5%) at low (high) masses for both channels. These moderate scale uncertainties are QED uncertainties. They are common to γinduced processes at LO (when the photon is inelastic) and can generically [99,100] be attributed to logarithmic and power-law terms in real-radiation matrix elements at NLO in QED, i.e., logarithmic and power-law terms in the tree-level qγ → qSS † matrix element. More specifically, the corrections are associated with tree-level q → qγ splittings and, after phase space integration, have the forms O log(M γγ /p γ T ) and O(p γ T /M γγ ). Like in QCD, real radiation diagrams can alternatively be described by employing so-called multi-leg matching (MLM) techniques. These would systematically combine, for example, the matrix elements for PDF uncertainties are about δσ AF (LO) 13 TeV ≈ ±2% over the mass range. As in the DY case, the ratio of the two AF rates are proportional to ratio of electric charges: This follows from the fact that the three point S − S † − γ vertex and the four point S − S † − γ − γ vertex are each proportional to the charge of k and h. Due to (a) the destructive interference between the γ and Z contributions in the DY channel, (b) the charge enhancement in the AF channel, and (c) the fact that the γγ luminosity is sourced from valance quark scattering whereas the DY process is sourced by valance-sea scattering, the γγ → k −− k ++ cross section consistently sits between the k −− k ++ and h − h + DY rates. The γγ → h − h + channel sits below all three curves. This suggests that a second k −− k ++ could be seen at the LHC before the h − h + channel. We caution that the similarity of the DY and AF rates for k −− k ++ production in the Zee-Babu model is mostly a consequence of the suppressed DY rate, not the "largeness" of the photon PDF. This is an uncommon occurrence but also not an artifact of the photon PDF. As discussed in Sec. 5, doubly charged scalars in other models carry different gauge quantum numbers; this often leads to larger DY rates in those models [55]. Claims that the cross sections for photoninduced processes readily exceed DY rates can usually be attributed to mis-modeling of photon PDFs or discounting large uncertainties. For dedicated discussions, see Refs. [98,99,101,102].
Gluon Fusion: We now consider k and h pair production from gluon fusion (GF): This loop-induced process is mediated by s-channel H 0 and Z bosons, as depicted in Fig. 2 With SM ZeeBabu NLO, we simulate the channel at LO, i.e., at one loop in α s , in mg5amc using We first note that the Z contribution vanishes due to two mechanisms: (a) Angular momentum conservation in the gg → Z * sub-graph causes the transverse component of the Z's propagator Π Z ρσ (q), i.e., the g ρσ term, to vanish when contracted with the quark loops. Table 3.
Ultimately, this can attributed to Z * behaving as if it is a pseudoscalar in the gg → Z * → X process [103], which is at odds with the parity conserving nature of the S − S † − Z coupling. A second comment is that the S − S − H 0 vertices differ only by a normalization. More specifically, from the Lagrangian in Eq. (2.7), the Γ H 0 −k−k and Γ H 0 −h−h vertices are The inputs of Eq.

4).
A third comment is that the QCD corrections to GF are typically large. This follows from both positive, virtual corrections and the opening of partonic channels in real corrections. To account for these, we apply a K-factor derived for the next-to-next-to-next-to-leading logarithmic threshold corrections (N 3 LL(thresh.)) to exotic scalar production in the Type II Seesaw [98] This scale factor captures the leading contributions to the GF channel at NNLO [104]. Using this K-factor is justified by the following: The processes in Eq. (4.39) and the analogous processes in the Type II Seesaw contain the same sub-graphs that are susceptible to QCD corrections; that is, the gg → H * /Z * sub-graphs are the same for all four processes. Therefore, the QCD corrections to the production cross sections are the same. Moreover, in Ref. [98], the K-factors at N 3 LL are reported to span K N 3 LL = 3.04 − 3.15 for scalar masses between 100 GeV and 800 GeV, with scale uncertainties reaching O(5% − 10%), when using the scale choices stipulated in Sec. 3.2. Approximating K N 3 LL ≈ 3.1 for this entire mass range leads to at most a ±2% over/underestimation of QCD corrections, which is within scale uncertainties.
Over the mass range m k , m h = 100 GeV − 1 TeV, the N 3 LL-corrected cross sections for k −− k ++ and h − h + pair production from GF span about At low masses, the k −− k ++ channel exhibits a cross section that is comparable to those of DY and AF. The rates for both GF channels quickly fall as masses increase. Beyond m k (m h ) ∼ 800 (400) [98,103]. PDF uncertainties in the LO rates are as low as δσ LO 13 TeV ∼ ±2% for low masses and as large as δσ LO 13 TeV ∼ ±20% for high masses.
Vector boson scattering: Finally, we consider VBS, which is given by the permutations of The diagrams for Zγ scattering are similar to those given for γγ and ZZ in Figs. 1 and 2. We include full interference by considering at the full, tree-level 2 → 4 process To simulate this at LO, the corresponding syntax in our Monte Carlo setup is generate qq qq > kk kk qq qq QED=4 QCD=0 output DirName7 generate qq qq > hh hh qq qq QED=4 QCD=0 output DirName8 While the photons in VBS are virtual and never on-shell, there is some phase space overlap with the AF channel due to PDF fitting and evolution. Moreover, as we include interference from all gauge-invariant diagrams, Eq. (4.46) formally includes diboson production, e.g., qq → SS † γ * → SS † qq, and associated Higgs processes, e.g., qq → ZH * → qqSS † . To minimize these additions and to regulate infrared divergences, we apply the following kinematic restrictions: Over the mass range investigated, the cross sections for the two processes roughly span The cross sections for both VBS channels sit well below those for AF but also differ qualitatively. At lower masses, the VBS rates are comparable to each other but bifurcate at higher masses. Summary: For representative scalar masses, we summarize the above cross sections, scale uncertainties, PDF uncertainties, and QCD K-factors at the √ s = 13 TeV LHC for all processes in Table 3. As an outlook for future experiments, we present the same results at a hypothetical √ s = 100 TeV VLHC in Fig. 3(b) and also in Table 3. (For the GF channel, we use K N 3 LL = 2.60.) For brevity, we do not comment much on cross sections at higher energies. Aside from the obvious jump in parton luminosities, which manifests as higher production rates, the cross section hierarchy does not qualitatively change from √ s = 13 TeV. Likewise, uncertainties at 100 TeV do not qualitatively differ from 13 TeV outside extreme values of m k and m h .
In these figures and table, we quantify for various production mechanisms the dependence on scalar masses, collider energy, PDF choice, factorization/renormalization scale choice, and some QCD corrections. For the DY and AF channels at their orders of perturbation theory, there are no additional dependencies in their cross sections on Zee-Babu model parameters since they are controlled entirely by SM gauge couplings. The GF channels, gg → k ++ k −− and gg → h + h − , are controlled additionally by the top quark mass, which is well-measured, as well as the scalar couplings λ kH and λ hH , respectively. These couplings appear quadratically in cross sections, i.e., σ GF ∼ |λ| 2 . And aside from these, the two GF cross sections are actually identical at LO in EW theory. (This degeneracy is broken at NLO in the EW theory since k and h have different weak hypercharges.) Therefore, in some sense, the inputs of Eq. (4.42) and the subsequent differences in cross sections illustrate the dependence on the scalar couplings λ kH and λ hH . Similarly, the differences in the VBF channels in the high-mass limit can be understood as varying λ kH and λ hH . As discussed in Sec. 4.3, the two VBF channels are driven by λ kH and λ hH in the high-energy limit; this limit is partially triggered by taking the masses of k and h to be much larger than those of the W and Z. Moreover, as only the production of k ++ k −− and h + h − pairs are being considered (and not, say, same-sign k ±± k ±± pairs), λ kH and λ hH are the only scalar couplings to appear. Other processes must be considered to explore other scalar couplings.
5 Distinguishing the Zee-Babu and Type II Seesaw models at the LHC Doubly and singly charged scalars are not unique predictions of any one model. They are in fact integral to several scenarios [13,15,[17][18][19][20][21][22][23]. Therefore, if doubly charged scalars are discovered at the LHC, work must be done to discern their nature, e.g., what are their gauge quantum   numbers, decay rates, and coupling strengths. In this section, we explore several ways one can potentially distinguish exotic scalars in the Zee-Babu model from those in the Type II Seesaw.
This section contains the main results of our study. In summary, we find that one must rely on decay correlations of exotic scalars to distinguish the models; most (but not all) productionlevel observables are too similar to be of use. We start with Sec. 5.1, where we compare production cross sections of exotic scalars. (Similar results at LO have been reported before [55]; however, our improved numerical analysis permits us to make new statements.) In Sec. 5.2, we compare predictions for differential distributions of doubly charged scalar pairs up to NLO+LL(PS) accuracy. We then reinterpret constraints on the Type II Seesaw from the LHC [51] in terms of the Zee-Babu model in Sec. 5.3. We discuss in Sec. 5.4 correlations in exotic scalar decays. And in Sec. 5.5, we discuss some criteria for establishing LNV in the two scenarios.

Total Cross Section
We start by briefly summarizing the relevant principles of the Type II Seesaw model. In the notation of Ref. [98], the model extends the SM by a single complex scalar multiplet∆ that carries the quantum number assignment (1, 3, +1) under the SM gauge group G SM (see Sec. 2). In terms of U(1) EM states,∆ and its vev v ∆ are given bŷ Conventionally,∆ is assigned lepton number L∆ = −2, making its Yukawa couplings to SM leptons LN conserving. The kinematic term and covariant derivative of∆ are given by After EWSB, one obtains the mass eigenstates ∆ ±± and ∆ ± . The state ∆ ±± is completely aligned with the gauge state∆ ±± , while ∆ ± is partially misaligned with the gauge state∆ ± . ∆ ± is an admixture of∆ ± and the Goldstone boson G ± , which decouples when v ∆ vanishes. Regarding quantum number assignments,∆ is defined with hypercharge Y∆ = +1. This means that the particles ∆ ++ and ∆ + carry weak isospin charges (T L ) 3 = +1 and 0, respectively, and antiparticles carrying the opposite charges. The gauge states k and h in the Zee-Babu model are defined with Y k = −2 and Y h = −1, so the particles k −− and h − carry negative electric charges, antiparticles carry positive electric charges, and all states have (T L ) 3 = 0. This distinction is why∆ is assigned L∆ = −2 while k and h are assigned L k , L h = +2.
Given this, we first show in Fig. 4(a) the production rate for k −− k ++ and h − h + pairs in the Zee-Babu model as a function of mass through various mechanisms at the √ s = 14 TeV LHC, i.e., the 14 TeV analogue of Fig. 3. In comparison, we show in Fig. 4(b) the √ s = 14 TeV LHC cross sections for ∆ ±± associated and ∆ ++ ∆ −− pair production in the Type II Seesaw at various accuracies, as adapted from Ref. [98]. Qualitatively, there are several similarities between the two models: (i) Both models feature pair production of doubly and singly (not shown) charged scalars through the neutral current DY mechanism. (ii) Both models feature pair production of charged scalars through AF. (iii) Both models feature pair production of charged scalars through GF. (iv) Both models feature pair production of charged scalars through VBF.
There are also three qualitative distinctions: (i) Unlike the Zee-Babu model, the Type II Seesaw admits associated production via the charged current DY process, i.e., qq → W ± * → ∆ ±± ∆ ∓ , since the ∆ ±(±) couple directly to the W boson. (ii) As the ∆ ±(±) belong to a multiplet in the Type II Seesaw, the existence of CP-even and -odd states ∆ 0 and ξ 0 implies additional production channels (not shown). (iii) Due to the different gauge couplings and multiplet states, many more sub-processes are present in VBS in the Type II Seesaw than in the Zee-Babu case.
The qualitative similarities listed above are also quantitative similarities. Explicit computation shows that the AF and GF cross sections for doubly and singly charged in the two scenarios are the same. (This assumes that masses S − S † − H 0 couplings, etc., are set equal.) Likewise, QCD corrections for these processes are the same; EW corrections are anticipated to be different due to the different gauge quantum numbers of the scalars but this must be investigated.
To further study the quantitative similarities of the neutral current DY channels, we define the scale factors ξ and χ as the ratio of the Zee-Babu and Type II hadronic cross sections: For all cross sections, we use the total rate at NLO in QCD. Conservative scale uncertainties are obtained by taking the ratio of extrema as allowed by scale variation, e.g., max(ξ) = max(σ DY (NLO) ). In Fig. 5, we plot the scale factors ξ (black band) and χ (teal band) with scale uncertainties (band thickness) as a function of scalar mass, fixing all equal, at Here, Q f and (T f L ) 3 are, respectively, the electric and weak isospin charges of f . The terms P γγ , P γZ , and P ZZ denote the pure photon, interference, and pure Z contributions to SS † production. In the following we work in the limit (M 2 Z /Q 2 ) → 0, and implicitly take M Z = 0.
First, note that the entire m S dependence in the partonic cross section is contained in the momentum/phase-space factor β. For fixed flavor f , scattering scale Q, and mass, the ratio of partonic cross sections, i.e.,σ DY (LO) (f f → k −− k ++ )/σ DY (LO) (f f → ∆ ++ ∆ −− ), is independent of m S . At the hadronic level, this holds under reasonable approximations: Assuming the production of TeV-scale SS † is driven by valence-sea annihilation, that at large momentum fractions the up-flavor PDF is twice as large as the down-flavor PDF, i.e., f u/p ≈ 2f d/p , and sea densities are equal, i.e., f u/p ≈ f d/p , then the ratio of hadronic cross section is proportional to Nominally, the dependence on m S in the ratio cancels. (There is a small dependence on m S for the doubly charged case that we address below.) The ratio is also stable at NLO in QCD. For instance: at NLO in QCD, virtual and soft corrections factorize (see, for example, Ref. [105]) and cancel. Under our assumptions, PDF subtraction terms also follow this behavior.
Moving onto the value of the ratios themselves, we note that the quantities (Q 4 P V V ), for V ∈ {γ, Z}, depend only on gauge quantum number when M Z can be neglected. For the uu → k −− k ++ channel, the γ − Z interference (Q 4 P γZ ) and pure Z (Q 4 P ZZ ) contributions strongly cancel, resulting in an O(−10%) correction to the pure γ contribution (Q 4 P γγ ). (This is a gauge-dependent statement and can be interpreted differently.) For the dd → k −− k ++ channel, the cancellation is stronger since Q 4 P γγ < Q 4 |P γZ | < Q 4 P ZZ . The cancellation between the γ−Z interference and pure Z contributions is an O(+6%) addition to the pure γ term. Essentially, f f → k −− k ++ can be treated as a QED process, which is consistent with the k − k − Z coupling being Weinberg angle-suppressed (see Sec. 2). For f f → ∆ ++ ∆ −− , the nonzero weak isospin charge induces constructive interference among the three terms for both up and down flavors.
All inputs equal, the ratio of partonic cross sections for doubly charged scalars simplify to ratios of "P V V -terms." For flavor combinations (uu) and (dd), these are given by We find good numerical agreement between this and the full matrix element calculation. Parameterizing the relative (uu) and (dd) contribution naïvely as (2/3) U and (1/3) D gives This agrees remarkably well with the numerical results reported in Fig. 5.
A closer inspection shows that the central value of the ratio ξ sits just below (above) the ξ ≈ 0.411 guideline at the lowest (highest) masses. This follows from PDF dynamics, namely deviations from our crude assumptions that f u/p = 2f d/p and f u/p = f d/p . For instance: masses that are O(100−300) GeV require momentum fractions that are O(0.01−0.05) at √ s = 13 TeV. This is where an asymmetry occurs between the u and d PDFs, with f d/p > f u/p . Also in this range, the u PDF is only O(20% − 40%) larger (µ f = 2m s ) than the d. This means that (dd) annihilation occurs more frequently than naïvely argued and pulls down the scale factor ξ. At larger momentum fractions, i.e., beyond O(0.1), the u − d asymmetry closes and the u/d ratio increases. Subsequently, (uu) annihilation occurs more frequently and pulls up ξ. The estimate ξ remains within or at the edge of the scale uncertainty band. PDF uncertainties for individual cross sections (lower panel) are at the ±5% − ±15% level for the masses investigated. Therefore, while deviations from ξ can evolve with improved PDF fits, the change will not be significant.
At 100 TeV (Fig. 5(b)), this behavior is unchanged for the masses under investigation since the same regions of individual PDFs are being probed. For example: the momentum fraction at √ s = 100 TeV. This is a manifestation of Bjorken scaling. Therefore, the (small) differences between the curves at 13 and 100 TeV are due to DGLAP evolution, which is logarithmic. (As stipulated in Eq. (3.3), PDFs are evolved up to factorization scales that scale as µ f ∼ 2m S .) Turning to singly charged scalars, note that the gauge charges for h − and ∆ − are both Q S = −1, Y S = −1, and (T S L ) = 0. Therefore, for fixed f and m S , the partonic cross sections for f f → SS † are the same. It follows that the ratio of hadronic rates is unity for all masses: This behavior is reflected in Fig. 5 at both collider energies. A caveat of this result is that cross sections are obtained at LO in the EW theory. Since h ∓ is an SU(2) L singlets but ∆ ± belongs to a triplet representation, it is likely that EW corrections can break this degeneracy.
Discussion: If k −− k ++ and/or h − h + pairs are discovered at the LHC, one could arguably use the observed cross section to extract gauge quantum numbers. However, k and h are shortly lived (see Sec. 4.2) and readily decay to SM particles. Therefore, what is actually measured is the combination of production and decay rates. As shown in Table 2, Eq. (4.17), and Eq. (4.20), the branching rates of k and h are sensitive to the relative sizes of masses and µ L , not just Yukawa couplings. Factors of ξ ∼ 0.4 can easily be absorbed by a branching rate. This last statement is also true for the Type II Seesaw. Decay rates of charged scalars in the Type II Seesaw are also correlated with neutrino oscillation parameters [106]. However, the uncertainty oscillation parameters remain sufficiently large to effectively absorb factors of ξ [98].

Kinematic distributions of doubly charged scalars
Beyond total cross sections, it is also possible to compare Zee-Babu and Type II scalars at a differential level. One argument goes that since the k ∓∓ and ∆ ∓∓ (or h ∓ and ∆ ∓ ) carry different gauge quantum numbers, and hence couple to the intermediate γ/Z with different strengths, then one may anticipate differences in differential distributions. We report that this argument does not work in the present case. Even at the differential level, we find that kinematic distributions of scalars produced by the DY process in the Zee-Babu and Type II models have the same shape and differ by only a normalization; the normalization is given by the ratio of hadronic cross sections. This finding also holds for the GF and AF channels.
To demonstrate this, we simulate at NLO+LL(PS) for √ s = 13 TeV the two DY processes 5 following the methodology in Sec. 3. This is the first time that kinematic distributions of the Zee-Babu model beyond LO+LL(PS) have been reported. We normalize the total Zee-Babu cross section to Type II cross section using the mass-dependent scale factor ξ(m s ): We also show (lower panel) the bin-by-bin ratio of the Type II and scaled Zee-Babu differential rates to the unscaled Zee-Babu prediction. The unscaled rate is our baseline. This is given by We report results for m S = 500 GeV in Fig. 6 and for m S = 1250 GeV in Fig. 7. The results are based on samples with N = 400k events each; statistical uncertainties are denoted by crosses +. Closed bars ][ denote scale uncertainties. We focus on k ++ and ∆ ++ but by momentum conservation the kinematics of k −− , ∆ −− mirror those for the positively charged states.  We start with Fig. 6(a). There we show (top) the transverse momentum distribution p S T of positively charged scalars S ∈ {k ++ , ∆ ++ }. As a function of increasing p T , the distributions rise to a maximum at about p T ∼ 400 GeV, or (p T /m S ) ∼ 0.8, and then fall with a power-law-like behavior. The spectra become small at small (p T /m S ) since SS † pairs are not produced precisely at threshold but rather with modest momenta. The qualitative behavior of all three distributions is the same but the unscaled Zee-Babu curve sits below the scale Zee-Babu curve and the Type II curve. Scale uncertainties (middle) for all three curves reach about ±4% at small and large p T ; at intermediate p T , scale uncertainties reduce to the sub-percent level. In comparison to the baseline (bottom), the bin-by-bin normalization of the Zee-Babu and Type II models are statistically indistinguishable. This is the first indication that kinematic distributions of doubly charged scalar production in the two models are identical, up to an overall normalization.
In Fig. 7(b) we show the rapidity y of S, defined by One sees that a bulk of the phase space sits in the range |y| < 1.5, indicating a longitudinal momentum small compared to the total energy carried by S. In other words, SS † pairs produced at the LHC are largely high-p T objects with only moderate longitudinal momentum. The scale uncertainty ranges from about ±2% at large rapidities (|y| > 1.75) but reduce to the sub-percent level at small rapidities (|y| < 0.5). In comparison to the baseline, the bin-by-bin normalization between the scaled Zee-Babu distribution and the Type II distribution are statistically indistinguishable. At low (high) rapidity, the Type II curve slightly juts above (under) the scaled Zee-Babu curve. However, the statistical uncertainty bars overlap for all y.
To explore angular correlations, we consider in Fig. 7(c) the azimuth separation (∆Φ) between S and the underlying hadronic environment. This is given by where p Had. is the vector sum over all hadrons within a maximum rapidity of y max = 10, i.e., We impose |y i | < y max to exclude beam remnant and simplify our generator-level analysis.
(Generally, such objects have little-to-no impact on transverse kinematics.) We observe that ∆Φ is a mostly flat distribution, with a slight monotonic increase as one goes from a parallel orientation (∆Φ = 0) to a back-to-back orientation (∆Φ = π). The means that it is more (less) likely for S and the hadronic activity to propagate in opposite (same) transverse direction. To understand this, note that the transverse part of p Had. is also the recoil of the (SS † ) system: Hence, the angle ∆Φ S Had can be interpreted as an azimuthal rotation of S relative to the (SS † ) system in the transverse plane. In the limit that | q T |/Q → 0, where Q = q 2 is the invariant  mass of the (SS † ) system, Born-like kinematics are recovered and the ∆Φ S Had distribution is flat, i.e., there is no dependence on ∆Φ S Had in 2 → 2 scattering. As | q T | grows, the (SS † ) system, and hence S, recoils more against the hadronic activity. Nonzero | q T | induces back-to-back separation in the transverse plane and simultaneously suppresses same-direction propagation.
Notably, ∆Φ S Had. is only accurate to LO+LL(PS) in our simulations. When the qq → SS † matrix element is known at LO, ∆Φ S Had is ill-defined as the (SS † ) system carries no p T . The p T of the (SS † ) system is generated first at LL accuracy by the parton shower and eventually at LO accuracy by the real radiation correction at NLO in QCD. Despite this formally lower accuracy, the scale uncertainties are at the sub-percent level. Again, the bin-by-bin normalizations of the scaled Zee-Babu distribution and Type II distribution are statistically indistinguishable.
In Fig. 6(d), we plot the polar distribution (cos β SQ ) of S in the frame of the (SS † ) system relative to the propagation direction of the (SS † ) system in the lab frame. Defining p to be the momentum of S in the (SS † ) frame, the observable is given symbolically by The distributions exhibit NLO+LL(PS) accuracy since q has longitudinal momentum in the lab frame, event at LO. We observe that all three curves obey a dσ ∼ (1 − cos 2 (β SQ )) distribution. This follows from angular momentum conservation: Imagining the decay of massive, virtual photon γ * → SS † , the Feynman rules of Eq. (2. 3) indicate that the corresponding helicity amplitude describes a p-wave process with −iM(γ * λ → SS † ) ∼ sin(β SQ ) for either transverse polarization of γ * λ , and where theẑ-axis is aligned withq. At the square level, one obtains λ=± |M λ | 2 ∼ sin 2 (β SQ ) = 1 − cos 2 (β SQ ) . (5.23) Scale uncertainties reach as large as ±1.5% in backward (cos β SQ = −1) and forward (cos β SQ = +1) regions, and are below ±1% in the central region (cos β SQ = 0). The bin-by-bin normalizations of the scaled Zee-Babu and Type II distributions are statistically indistinguishable. In Fig. 7, we plot the same observables as in Fig. 6 for the benchmark m S = 1250 GeV. Qualitatively, we find strong similarities between the two mass choices. Quantitatively, the absolute normalizations of the differential distributions are smaller than the previous case due to naturally the smaller production cross section. Beyond this, shape broadening or narrowing can be attributed to the larger mass scale. For the p T and y distributions, we find a slightly larger residual scale uncertainty, but also find that uncertainties stay below ±4%. Importantly, as in the low-mass case, the Zee-Babu and Type II distributions are statistically indistinguishable.
For completeness, we consider two measures of the hadronic activity in SS † pair production to demonstrate that the underlying event also remains unchanged between the Zee-Babu and Type II scenarios. For (a,b) m S = 500 GeV and (c,d) m S = 1250 GeV, we show in Fig. 8(a,c) the transverse momentum of the (SS † ) system, q T , and (c,d) the inclusive H T , defined by which is built directly from hadrons. Despite being LO+LL(PS) accurate, scale uncertainties reach only ±2%. As before, the scaled Zee-Babu and Type II distributions are indistinguishable.   Figure 9. (a) Estimated and projected cross section limits at 95% CL on k −− k ++ → 4 production, ∈ {e, µ}, in the Zee-Babu model at the √ s = 13 TeV LHC with L ≈ 139 fb −1 and L ≈ 3 ab −1 of data, respectively, as derived from constraints on ∆ ++ ∆ −− → 4 production in the Type II Seesaw by the ATLAS experiment with L ≈ 139 fb −1 [51]. (b) Same as (a) but for the k ∓∓ → ∓ ∓ branching rate.

Limits and projections for the LHC
As demonstrated above, total and differential cross sections for pair production of charged scalars in the Types II and Zee-Babu models differ at most by an overall normalization. Consequentially, their decay products will inherit this sameness and also exhibit nearly identical kinematics.
Despite this hardship, there a silver lining of this sameness: the selection (ε) and acceptance (A) efficiencies obtained by LHC experiments in searches for charged scalars in the Type II Seesaw are automatically applicable to the Zee-Babu model. This is a nontrivial conclusion. It implies that for common final states the two models can be tested simultaneously at the LHC without additional event generation or additional signal/control/validation-region modeling.
Normally, new signal events must be simulated to reinterpret or recast a collider analysis for one scenario in terms of a second scenario. That is not needed for the charged scalars in the Zee-Babu and Type II models. For example: suppose that after all acceptance and selection requirements the total number of Type II events at a given integrated luminosity L is (5.25) Since the selection and acceptance efficiencies are the same, which follows from final states having identical kinematics, the corresponding number of Zee-Babu events is  Table 4. For representative m k (top), estimated (middle) and projected (bottom) limits on the decay rate of k ∓∓ → ∓ ∓ set by ATLAS at 95% CL using L ≈ 139 fb −1 [51] and L = 3 ab −1 of data.
One can identify the cross-section ratio in Eq. (5.26) as the scale factor ξ(m S ) in Eq. (5.14).
Since the final states are assumed to be the same, the backgrounds are the same. And after acceptance and selection cuts, the number of background events are the same. Subsequently, the upper limit on a cross section derived for one scenario is the same as for the other. The specific parameter space excluded is, of course, different for the two scenarios.
To demonstrate this we consider the search for ∆ −− ∆ ++ pairs in the 4 ± channel, where ∈ {e, µ}, by the ATLAS experiment using the full Run II data set [51]. In Ref. [51], a mass-dependent 95% confidence level (CL) limit (σ ATLAS 95% CL ) is reported on the (unfolded) quantity (5.28) Applying the same limit to the (unfolded) quantity in the Zee-Babu model we obtain the limit shown in Fig. 9. Assuming branching rates of unity, k ∓∓ with masses m k < 890 GeV (5.30) are excluded by ATLAS at 95% CL using approximately L ≈ 139 fb −1 of data at √ s = 13 TeV. Assuming that sensitivity σ ATLAS 95% CL scales with the square root of integrated luminosity, i.e., and unchanged selection and acceptance efficiencies, we also show the projected sensitivity with L ≈ 3 ab −1 of data at √ s = 13 TeV. For branching rates of unity, k ∓∓ with masses m k < 1110 GeV (5.32) can be excluded at 95% CL. This sensitivity can be improved with higher collider energies, larger data sets, and improved analysis techniques. (This would equally benefit searches for ∆ ∓∓ .) Invert Eq. (5.29), we obtain limits on the k ∓∓ → ∓ ∓ branching rate. This is given by In Fig. 9(b) we show the estimated and projected limits for L ≈ 139 fb −1 and L = 3 ab −1 . As summarized in Table 4, decay rates as small as 16% (74%) can be probed for m k = 400 (800) GeV with the full Run II data set. Tentatively, this can be improved two-fold at the HL-LHC.

Correlations in flavor-violating decays of singly charged scalars
We now discuss a consequence of the relationships between Yukawa couplings f and oscillation parameters as summarized in Eqs. (4.9) and (4.10). As discussed in Secs. 4.2 and 5.1, k and h are shortly lived and decay readily to SM particles. In principle, production rates themselves are not observed at the LHC but rather the combination of production and decay rates. For the production of h − h + pairs, which preferably decay via the h ∓ → ∓ ν channel, neutrinos cannot be flavor tagged at the LHC. One can only identify the charged lepton ∓ . Subsequently, measurements of the pp → h − h + → − + ν ν process implicitly sum over all neutrino states. Guided by this, we define the following pair of inclusive partial widths for h: Recalling that f ee , f µµ = 0, and neglecting O(m 2 /m 2 h ) terms in expressions for total widths gives which is the e-over-µ branching ratio for h. Inserting Eqs. (4.9) and (4.10), we obtain for the NO and IO of neutrino masses. These simple, analytic expressions are a direct consequence of Eqs. (4.9) and (4.10) but have not previously been reported in the literature. The utility of the branching ratio R, as oppose to the branching rate BR, is its independence of h's total width. Unlike individual branching rates, the ratio R h eµ is measurable at the LHC. It can be determined, for example, by comparing the event numbers (N ) in the processes: where E T is the missing transverse momentum of the event. Accounting for selection and acceptance efficiencies, integrated luminosity (L), the hadronic cross section σ , and assuming the narrow width approximation, then the branching ratio (squared) is obtained from the ratio R h eµ can also be extracted from the ratio of (e + µ − + p T ) and (µ + µ − + p T ) events. The e-over-τ and µ-over-τ branching ratios can also be constructed in the same manner. The relative smallness of the IO's uncertainty is due to its dependence on only two oscillation angles, one of which (θ 13 ) is well measured. Fixing either θ 12 or θ 23 in the NO formula to its central value roughly halves the uncertainty. Given the certainty in oscillation parameters, the predictions are sufficiently robust to discriminate between the NO and IO. For completeness, we report in Table 5 the oscillation angles and phase at which the R h eµ are maximized/minimized.
New correlations in low-energy transitions from oscillation data: As stipulated in Sec. 4, a comprehensive discussion on the phenomenology of the Zee-Babu model in low-energy processes is outside the scope of this study. (These can be found in Refs. [28,29,35,37,39].) Nevertheless, in light of Eqs. (5.35) and (5.36), which have not previously been reported, it is worth commenting briefly that new connections can also be established between flavor-violating, low-energy transitions and oscillation data.
For instance: in the Zee-Babu model, → νν decays are additionally mediated at treelevel by the charged scalar h ± and the f couplings. Experimentally, this manifests as a Fermi constant (G ZB F ) that is shifted from its SM value (G F ). Explicitly, this shift is given by [35,107] where G F is assumed to be extracted from, e.g., decays of hadrons. This result holds so long as m h remains large compared to the vev of the SM Higgs. The extracted Fermi constant from various → νν decay channels then isolates the antisymmetric couplings: Taking this a step further, the ratio of these differences gives the ratio of f couplings: Once again, using Eqs. (4.9) and (4.10) for NO and IO, respectively, one obtains correlations for flavor-violating transitions that are fixed by oscillation data. Again, these expressions have not previously been reported in the literature. Such correlations can also be established in → γ transitions. However, further investigations, including numerical studies, are left to future work.

Establishing non-conservation of lepton number at the LHC
If neutrino masses are generated through the Zee-Babu mechanism [21][22][23], then neutrinos are Majorana fermions and LN is not conserved in scattering and decay processes. Whether this can be established through low-energy experiments, such as neutrinoless ββ decay, it is paramount to establish if LN is violated at other energies. Among other reasons, LNV is predicted in many neutrino mass models and correlating LNV across the different scales can provide critical guidance on the underlying theory. We now discuss a strategy for establishing LNV in the Zee-Babu model at the LHC. For reference, we start an analogous strategy in the Type II Seesaw.
Establishing LNV with Type II scalars: Under conventional quantum number assignments, LN is broken explicitly in the scalar potential of both the Type II and Zee-Babu models by a dimensionful parameter µ L . In both scenarios, the Yukawa couplings to exotically charged scalars are LN conserving. After EWSB, the triplet scalars of the Type II Seesaw acquire a vev v ∆ proportional to µ L , which manifests in the three-point ∆ − ∆ − W coupling. To establish LNV in the Type II Seesaw, one can search for the following LHC processes and work by contradiction. The argument goes as follows: Assuming that LN is conserved, the four-lepton channel establishes that the ∆ ∓∓ states each carry L = ±2 since each ∓ carries L = ±1. The second channel is mostly reconstructable, particularly if employing techniques pioneered for extracting neutrino momenta in leptonic decays of top quark pairs [108][109][110]. The channel establishes that the ∆ ∓∓ states carry of L = 0 since W ± carry L = 0. This leads to a contradiction that LN is conserved. (Technically, the v ∆ in the W − W − ∆ vertex carries away |L| = 2.) Along these lines, the Zee-Babu model does not contain the W − W − S vertex. Therefore, observing the channel would signal that the Zee-Babu model is not realized.
Establishing LNV with Zee-Babu scalars: The above strategy does not apply to the Zee-Babu model since neither k nor h couples to the W . And based on the size of µ L , LN-violating processes may be inaccessibly at the LHC even if possible elsewhere. Assuming relevant decay processes are accessible, one can establish LNV by searching for the following LHC processes The argument, which also works by contradiction, goes as follows: Suppose LN is conserved.
Observing the four-lepton channel (Eq. (5.46a)) establishes that each k ∓∓ carries L = ±2 since each ∓ carries L = ±1. Observing the opposite-sign dilepton channel (Eq. (5.46b)) establishes that h ∓ carries L = 0 or L = ±2. This assertion requires a few clarifying remarks. First, the signature pp → + − E T features a large multi-boson and top quark background, and will be difficult to observe. However, considering the extent to which SM simulations at NNLO+LL(PS) can describe data [111][112][113][114], the presumable mass difference between SM and Zee-Babu particles [28,35], and prospective search strategies [52,56,115], we premise this is attainable. Second, it is possible to show that the signature's E T is driven by two neutrinos since (a) the cross-section ratio of pp → + − E T with respect to different lepton flavors is fixed by oscillation data (see Sec.5.4), and (b) the kinematic distributions of , , and E T are constrained since h ∓ → ∓ ν is a two-body decay involving (approximately) two massless states. Therefore, extending NNLO+LL(PS) technology for pp → W + W − → + − E T in the SM [111][112][113] to h − h + pair production provides a (limited) means of checking whether the E T is driven heavier states, some light dark-sector fermion, or by more than two states. Failure to satisfy (a) and (b) would suggest that the Zee-Babu model is not realized. Even if satisfied, one can only assert that each h ∓ carries L = 0 or L = ±2 since it is impossible to check the LN of outgoing neutrinos.
Finally, since each h ∓ carries L = 0 or L = ±2, observing the four-lepton and E T channel (Eq. (5.46c)), which again can be checked using kinematic distributions and ratios of cross sections, establishes that each k ∓∓ carries L = 0 or L = ±4. This is in contradiction with the four-lepton channel (Eq. (5.46a)), which establishes L k = ±2, and implies that LN is not conserved. In the event that m h > m k , then the kinematically suppressed channel shows that LN cannot be conserved since the h ∓ → k ∓∓ h ± * splitting suggests L k = 0 or ±4.

Summary and Conclusions
With widely felt impact in nuclear physics, astrophysics, and cosmology, the origin of neutrinos' tiny masses and large mixing is among the most pressing mysteries in particle physics today.
Establishing whether neutrinos are their own antiparticles, implying that LN is not conserved, is also fundamental to model building. Naturally, there are numerous models of increasing complexity that answer these questions and are also testable at ongoing and near-future experiments. Among these scenarios are the Type II Seesaw and Zee-Babu models for neutrino masses, which, less commonly, can reproduce oscillation data without invoking sterile neutrinos. Both scenarios hypothesize the existence of exotically charged scalars that couple directly to the SM Higgs and SM gauge bosons, and therefore can be produced copiously at the LHC if kinematically accessible. In this study, we have revisited the phenomenology of the Zee-Babu model (Sec. 4) and focused (Sec. 5) on the ability to distinguish singly and doubly charged scalars from the two models at the LHC. We conclude that this task is much more difficult than previously believed.
After reviewing the tenets of the Zee-Babu model (Sec. 2), and after presenting updated cross section predictions for k −− k ++ and h − h + production at the LHC through various mechanisms up to NLO in QCD (Sec. 4.4), we compared total (Sec. 5.1) and differential (Sec. 5.2) predictions for the Zee-Babu and Type II Seesaw models. All inputs equal, we find that total and differential rates for producing pairs of doubly and singly charged scalars are identical in shape and differ by a normalization equal to the ratio of hadronic cross sections, which can be unity. This holds for the Drell-Yan, γγ fusion, and gg fusion, as well as observables at LO+LL(PS) and NLO+LL(PS) in QCD. Importantly, the differences in normalizations are sufficiently small that they can be hidden by unknown branching rates or unknown couplings to the SM Higgs. This similarity allows us to reinterpret LHC constraints and projected sensitivity on doubly charged scalars decaying to leptons from the Type II Seesaw in terms of the Zee-Babu model (Sec. 5.3).
Outlook: Despite potential hardships, there is some guidance on distinguishing the two models: Unlike the Type II Seesaw, the Zee-Babu model predicts one massless neutrino, and therefore features a clearer prediction for the rate of neutrinoless ββ decay. Aside from this, charged scalars in the Zee-Babu model do not couple to the W boson at tree level. This means that the Type II Seesaw predicts several associated production channels at the LHC not found in the Zee-Babu model. If the Zee-Babu model is realized by nature, then these channels are absent. Furthermore, neutrino oscillation parameters are now sufficiently precise to make clear predictions (Sec. 5.4) for branching ratios, i.e., ratios of branching rates, of charged scalars in the Zee-Babu model. Such observables are less sensitive to unknown decay rates of charged scalars and are presented for the first time in Sec. 5.4. Similarly, new correlations between oscillation data and searches for lepton flavor violation at low-energy experiments, such as → νν decays, can also be established. Finally, the inherent differences in the two models require different strategies for establishing LNV at the LHC (Sec. 5.5). Finally, it is also possible that NLO in EW corrections to the production rates of charged scalars at hadron colliders can help break degenerate predictions. In all these directions we encourage and anticipate future exploration.