The minimal stealth boson: models and benchmarks

Stealth bosons are relatively light boosted particles with a cascade decay $S \to A_1 A_2 \to q \bar q q \bar q$, reconstructed as a single fat jet. In this work, we establish minimal extensions of the Standard Model that allow for such processes. Namely, we consider models containing a new (leptophobic) neutral gauge boson $Z'$ and two scalar singlets, plus extra matter required to cancel the $\text{U}(1)'$ anomalies. Our analysis shows that, depending on the model and benchmark scenario, the expected statistical significance of stealth boson signals (yet uncovered by current searches at the Large Hadron Collider) is up to nine times larger than for the most sensitive of the standard leptophobic $Z'$ signals such as dijets, $t \bar t$ pairs or dibosons. These results provide strong motivation for model-independent searches that cover these complex signals.


Introduction
New heavy resonances are easy to spot at the Large Hadron Collider (LHC) when they decay into charged leptons, e.g. Z → e + e − /µ + µ − , W → eν/µν, but they are quite more difficult to detect in hadronic final states, since the production of quarks and gluons by QCD interactions has a very large cross section. Still, heavy resonances decaying into boosted hadronically-decaying W , Z or Higgs bosons, or top quarks, may be separated from the background. In the last decade, great progress has been made in this direction with the development of jet substructure techniques [1][2][3][4][5][6][7][8][9][10][11][12][13] and grooming algorithms [1,[14][15][16]. These tools allow to distinguish jets originating from boosted hadronically-decaying bosons and top quarks from the Standard Model (SM) background, composed mainly by quark and gluon jets produced in QCD processes. In this way, searches for diboson [18][19][20][21][22][23][24], 1 jaas@ugr.es 2 filipe.joaquim@tecnico.ulisboa.pt tt [25] and tb [26,27] resonances in purely hadronic channels have been performed, with a sensitivity that turns out to be competitive with channels involving leptons in the final state. Nevertheless, these searches are insensitive to heavy resonances decaying into more complex hadronic final states, giving rise to multi-pronged jets.
One example of a multi-pronged jet signature is given by the 'stealth bosons' introduced in ref. [28], which are boosted particles S decaying hadronically into four collimated quarks, via two (equal or different) intermediate particles A 1 , A 2 , namely The particles A 1,2 in the above decay chain may be SM weak bosons W , Z, a Higgs boson, or new relatively light (pseudo-)scalars. When S is produced in the decay of a much heavier parent resonance R, (with X an additional particle) its experimental signature is a fat jet with four-pronged structure. Jet substructure observables designed to distinguish two-pronged Z, W and Higgs decays from the QCD background, for example the so-called D 2 [11] and τ 21 [7,10] variables respectively used by the ATLAS and CMS Collaborations, classify four-pronged jets as QCD-like. Therefore, should a new resonance involve one or more decay products of this type, it would be very hard to identify it in current searches. On the other hand, generic searches that use a multivariate tool like an anti-QCD tagger [29] to pin down multi-pronged jets from the QCD background are sensitive to this type of signals. Notice that if S weakly couples to SM particles, for example if it is a neutral scalar, its direct production cross section may be too small for this particle to be directly observed.
The aim of this paper is to investigate the minimal SM extension in which stealth boson signals may appear, and contextualise the relevance of these signatures as a discovery channel for new leptophobic resonances, when compared to the usual decay modes searched for at the LHC, like dijets, dibosons or tt pairs. In section 2 we find, following a bottom-up approach, that the minimal additional content that allows for the cascade decays in (1) is a Z boson and two scalars that are singlets under the SM group but charged under the extra U(1) (further details of the models are given in appendix A). An extension with one Z boson, a new scalar doublet and a scalar singlet, which is also attractively simple, does not serve our purposes, as briefly discussed in appendix B. Section 3 is devoted to the discussion of how benchmark scenarios for the scalar masses and mixings are tested, to ensure that the reconstructed model parameters lead to an absolute minimum of the scalar potential. In section 5 and appendix C those benchmark scenarios are studied in detail performing fast simulations of the various Z signals in the decays into stealth bosons, dijets, and tt, as well as their SM backgrounds. We discuss our results in section 6.

The minimal stealth boson models
Following a minimalistic approach, we assume that the heavy resonance R in (2) is a neutral colour-singlet Z boson, so that the gauge symmetry of the SM is extended by an extra U(1) . 3 We require the Z boson to be leptophobic, i.e. the left-handed lepton doublets L and right-handed singlets e R have zero hypercharge Y = Y e = 0 under the new U(1) . Otherwise, the leptonic signals Z → e + e − and Z → µ + µ − would be easy to observe at the LHC. Gauge invariance of the Yukawa couplings with the SM Higgs doublet Φ, then requires that Y Φ = 0, and that the left-handed quark doublets q L and right-handed quark singlets u R , d R have the same hypercharge Y q = Y u = Y d ≡ z (for simplicity we omit generation indices and generically denote the Yukawa couplings by y u , y d , y e ). The hypercharge assignments are collected in table 1, where z is unspecified.
Cancellation of the anomalies associated to U(1) requires introducing extra matter, which we assume to be vector-like under the SM gauge group, to preserve SM anomaly cancellation. Two simple choices for these extra degrees of freedom, which we will denote as model 1 and 2, are: • Model 1: One set of vector-like quarks, comprising a doublet (T 1 B 1 ) with SM hypercharge Y = 1/6 plus vector-like singlets T 2 , B 2 of charge 2/3 and −1/3, respectively.
The hypercharge assignments for these fields is summarised in table 2. We note that model 2, with z = 1/3, has been considered in previous literature [33], motivated by the search for an anomaly-free Z dark matter mediator (the dark matter particle corresponds to the singlet N 2 ) with weak constraints from direct detection experiments. A similar model, 3 Alternatives in the context of left-right models can easily be worked out from the results in Ref. [30]. In that work we focused on the 'resolved' signatures where three or four well-separated bosons are produced from the cascade decay of a W boson into new scalars. When the masses of the intermediate particles are lighter, their bosonic decay products are merged giving rise to signatures such as in (1). Cascade decays can also be produced in a variety of other non-minimal scenarios, see for example refs. [31,32].  with a three-fold replication of the new lepton set, (N i E i ) L , (N i E i ) R , N jL , N jR , E jL , E jL, , with i = 1, 2, 3 and j = 4, 5, 6, also preserves the cancellation of anomalies. In this case, the lepton hypercharges are 1/3 of the values quoted in table 2 for model 2. The phenomenology, except for the signals associated to the new fermions which we do not address here, is the same of model 1. Other more baroque possibilities exist for the choice of new fermions by, for example, introducing vector-like quark doublets with Y = 7/6 or Y = −5/6 and O ∼ 10 quark singlets. Notice that kinetic mixing would modify our hypercharge assignments but, since both the U(1) coupling g Z and the hypercharge parameter z are unspecified, it has no effect in our analysis and we do not consider it.
The scalar sector of the SM must be extended in order to break the U(1) symmetry and generate the Z boson mass. The simplest possibility is to consider a neutral complex SU(2) L singlet χ with non-zero hypercharge Y χ under U(1) . Having the Z mass generated by a higher SU(2) L multiplet is problematic, as its vacuum expectation value (VEV) would also contribute to the weak boson masses. Notice that the heavy fermion masses can also be generated with the same singlet provided Y χ = 3z (for model 1) or Y χ = 9z (for model 2), and that the new fermions do not have Yukawa interactions with the SM ones. Moreover, in model 2, if the lightest new fermion is a neutral singlet, it may possibly be a dark matter particle [33], while in model 1 the lightest new quark would have exotic signatures [34], not addressed here.
Additional scalars, besides this singlet, are required to yield the cascade decays (1) and (2) (ref. [33] only considers one scalar singlet). As discussed in appendix B, adding a second scalar doublet is not a viable option. Thus, we instead consider a scalar sector comprising the SM doublet Φ and two complex singlets χ 1 , χ 2 with the same hypercharge. Further extensions of the scalar sector that allow interactions of the new fermions with the SM ones are possible, but they are not required for our purposes. The most general gauge-invariant scalar potential is where V Z 2 (V Z 2 ) contains the terms which are invariant under (break) a Z 2 symmetry for which χ 2 → −χ 2 and all remaining fields transform trivially. Among all the above parameters, m 2 0 , m 2 11 , m 2 22 , λ 0−3 and λ 5,6 are real, while m 2 12 , λ 4 and λ 7−9 can be, in general, complex. We write the neutral scalar in Φ = (φ + φ 0 ) T and the singlets χ 1,2 as such that the VEVs are Rephasing χ 2 to χ 2 = e −iϕ χ 2 (which has real VEV χ 2 = u 2 / √ 2), the potential V can be written in terms of χ 2 as in (4) with the replacements while the remaining parameters stay invariant. Therefore, in the general complex case one can always assume ϕ = 0 without loss of generality in order to simplify the expressions, with a possible non-vanishing phase absorbed by the above redefinition. On the other hand, if the (unprimed) parameters in the potential are all real this cannot be done, and a non-zero phase ϕ could break CP spontaneously.
There are four minimisation conditions corresponding to the four parameters v, u 1,2 and ϕ, Since we will be interested in those vacuum configurations with u 1,2 = 0, we adopt the common definitions: The minimisation conditions are used to express m 2 0 , m 2 11 , m 2 22 and Im(m 2 12 ) as functions of the remaining potential parameters and VEVs. The two would-be Goldstone bosons are G 0 1 = η 0 and G 0 2 = cos β η 1 + sin β η 2 . The orthogonal state A 0 = − sin β η 1 + cos β η 2 is CP-odd, being a mass eigenstate if the parameters in the scalar potential are real and ϕ = 0. In the basis H i = (ρ 0 ρ 1 ρ 2 A 0 ), the squared mass matrix, denoted as M ij , has elements (with M ij = M ji ) We remark again that the minimum conditions and the expressions for M ij are valid both for a general potential with complex parameters, in which case one can just drop the primes and assume ϕ = 0, or for a potential with real parameters, in which case the primed parameters are defined by (7). We also note in passing that, should we have chosen χ 1 and χ 2 with different hypercharges, m 12 and λ 4−6 would vanish, in which case M i4 = 0 and A 0 would be massless.
The scalar interactions with the Z boson field originate from the term Since for the scalar doublet Y Φ = 0, there is no Z − Z mixing and B µ ≡ Z µ is a mass eigenstate, with mass We express the scalar weak eigenstates as H i = O ij H j , where H i are the mass eigenstates with mass M H i , and O is a 4 × 4 real orthogonal matrix. The Z H i H j (i < j) couplings are then with mixing factors Notice that R ij are anti-symmetric and therefore R ii = 0, reflecting the fact that Z → H i H i is forbidden. Also, it can be shown that i<j R 2 ij = 1 due to the orthogonality of the mixing matrix O. The Lagrangian for the interaction of the SM Z boson with two neutral scalars is with g W and c W being the weak coupling and the cosine of the weak angle, respectively. This term does not yield interactions among the Z and two physical neutral scalars since G 0 1 = η 0 . The three-scalar couplings in the mass eigenstate basis are complicated functions of the potential parameters, the VEVs and the mixing matrix O, but can generically be written as The constants λ ijk are totally symmetric under interchange of two indices, and their expressions are collected in appendix A. For convenience we introduce symmetry factors S ijk , obeying S ijk = 1 for three different indices, S ijk = 2 if two of them are equal, and S ijk = 6 if i = j = k. The coupling of the scalar mass eigenstates to the SM fermions f and weak bosons V = W, Z arise from their ρ 0 component, Notice that the interaction of H i to fermions is purely vectorial, even if they have a non-vanishing CP-odd component. This is so because the ρ 0 interaction with fermions is vectorial, and the matrix O is real.
By inspection of the mass matrix (10) one sees that it is rather easy to make our model compatible with experimental data. Taking λ 5,6,9 small, the mixing of SM-like Higgs boson H ≡ H 1 with the new singlets (given by O i1 , i = 1) can be made as small as desired, in particular fulfilling the current constraints [35]. The masses and mixing of the additional scalars depend on m 12 , u, v, λ 1−4 , λ 7−9 and β. If u is at the TeV scale, masses for the new scalars around the electroweak scale can naturally be obtained with small λ i couplings, without the need of fine-tuned cancellations. Mixing between the CP-even states ρ 1 , ρ 2 and the CP-odd one A 0 is possible with complex λ 4−6 without affecting the properties of the SM-like Higgs boson. If these parameters are real, A 0 is itself a mass eigenstate and H 2,3 are CP-even and an admixture of ρ 1 and ρ 2 .
The framework described above naturally accommodates the Z cascade decays we are seeking for. The decay widths of the Z boson into SM quarks and scalars are where N c = 3 is the quark colour factor and Y χ = 3z (9z) in model 1 (model 2). We have defined the kinematical function Neglecting the masses of the decay products, the total width of the Z boson into scalars is In the last equation, the symmetry factor (1+δ jk ) accounts for the presence of two identical particles in the final state when j = k. Since O 1j 1, the scalars will dominantly decay to lighter scalars, if kinematically allowed. Otherwise, they will decay into W + W − , ZZ or fermion pairs, like a Higgs boson with a mass M H i (decays such as H i → H j Z are absent). For lighter masses the dominant mode will be H i → bb.

Methodology for the parameter space scan
The number of parameters in the scalar potential (4) is large enough to reproduce any pattern of scalar masses and mixing. Thus, we will focus on setting benchmark scenarios representative of the signals we are interested in. Within this approach, and with the goal of reducing the number of parameters, we will consider a simpler version of the model with λ 7−9 = 0 in the scalar potential. This corresponds to having a softly broken Z 2 symmetry under which χ 2 → −χ 2 , i.e. the only term remaining in V Z 2 = 0 is the bilinear m 2 12 soft-breaking term. We are then left with twelve real parameters in the scalar mass matrix: Re(m 2 12 ), λ 0−3,5,6 , Re(λ 4 ), Im(λ 4 ), β, v and u (determined by the Z mass through eq. (12)), of which only ten are independent due to the relations These ten parameters match the four scalar masses and six independent parameters of the (real) 4 × 4 orthogonal scalar mixing matrix. The first equation in (22) determines one of the masses M H i through while the second one determines tan β. Taking O ij , v, M Z , λ 2 and three of the masses M H i as inputs, the remaining parameters in the potential can be determined as where M ij are expressed in terms of M H i and O ij using The 4 × 4 mixing matrix is paremeterised by the product of 2 × 2 rotations as being O kl a rotation in the (k, l) plane by an angle θ kl , whose (i, j) matrix element can be written as The constraints on the couplings of the SM 125 GeV Higgs boson (H ≡ H 1 in our models) obtained by the ATLAS Collaboration [35] imply O 2 21 +O 2 31 +O 2 41 ≤ 0.05 at 95% confidence level (CL), 4 in which case O 1j and O j1 are small for j = 1. This, in turn, implies that the mixing angles θ 12 , θ 13 and θ 14 are small. The strength parameter µ = 1.11 +0.09 −0.08 corresponds to O 2 11 in our notation, from which limits on the other matrix elements can be obtained using unitarity and approximating the Gaussian distribution by a symmetric one. Ignoring the fact that O 2 11 ≤ 1 due to unitarity, one arrives at the bound O 2 21 +O 2 31 +O 2 41 ≤ 0.05 at 95% CL. By restricting ourselves to the physical region O 2 11 ≤ 1, the 95% CL limit is relaxed to 0.11. In all our benchmark scenarios the matrix elements O 1j and O j1 (j = 1) are more than two orders of magnitude below these bounds.
we first determine M H 2 and tan β using eqs. (22). Afterwards, the parameters in (24) are computed using also eq. (25). At this point, for the chosen set of inputs, the potential is completely defined. It now remains to check whether stability and perturbative unitarity criteria are fulfilled, and ensure that our VEV corresponds to the global minimum of the potential. For the stability analysis we follow the method of refs. [36] and [37] based on requiring copositivity of the quartic coupling matrix Λ. Parameterising the field bilinears as defined in the basis (h 2 , h 2 1 , h 2 2 ). The stability of the potential is ensured if the above matrix is copositive, i.e., if the following conditions hold [36]: In the above relations ρ = 0, 1 depending on whether Re(λ 4 ) > 0 or Re(λ 4 ) < 0, respectively. Since in the cases we are interested in the quartic couplings λ i are typically very small (due to the fact that M H i /u 1), perturbative unitarity is automatically ensured. Thus, we do not perform the complete analysis of S-matrix unitarity for elastic scattering of two scalar boson states.
It now remains to check whether our minimum is the global minimum of the potential. For that, we must compare the value of the potential at our minimum, with the values at any other minima obeying the minimisation conditions (8). The most straightforward alternative solutions correspond to vacua with vanishing VEVs. Namely, we have where V i corresponds to the value of the potential at the corresponding set of VEVs. Notice that these two solutions must be discarded as being the global minimum of the potential since they imply no spontaneous symmetry breaking of the SM and/or the U(1) gauge symmetry. A class of nontrivial minima with ϕ = 0 and v, u 1,2 = 0 may exist. Since for these cases the analytical treatment is quite involved, we use a numerical routine to spot those solutions. For all alternative minima we check positivity of the scalar masses and if V i < V 0 . At the end, only those sets of input parameters which lead to a stable potential and to a global minimum are considered viable in the parameter space scan performed in the next section.

Stealth boson benchmarks
We are interested in scenarios with Z and H i masses close to those used for the anti-QCD tagger in ref. [29]. We remark that this assumption is done only for the sake of simplicity, and with the purpose of using the performance for signals and backgrounds obtained in previous work without the need of training neural networks for new taggers. Thus, we restrict our study to benchmarks where the Z boson decays into two stealth bosons, that is for example with H 2 subsequently decaying into quark pairs. Scenarios with Z → H 3,4 H 2 , H 3,4 → H 2 H 2 are also interesting and lead to signals that are quite elusive as well, since a light boosted H 2 → qq produces two-pronged jets that closely resemble one-pronged QCD jets. However, their analysis requires the development of new taggers, which is out of the scope of this work.
In the following we will identify three representative scenarios.

Scenario 1
GeV, as in one of the benchmark points used in the tagger labelled as std1000 in ref. [29]. The coupling is  (1) coupling is such that g Z z = 0.1. set to g Z z = 0.1. 5 We perform a scan of the allowed parameter space by varying θ 23 , θ 24 and θ 34 with a flat distribution (keeping the other mixing angles small as required by constraints on the couplings of the SM Higgs) and compute the Z branching ratio to scalars. The results for model 1 are presented in figure 1. The branching ratio for quark pairs (not summed over flavours) is included for comparison. For model 2 the branching ratios for scalars trivially scale by a factor of 4.8, and the branching ratios to quark pairs by 0.53.
The results show that the decay Z → H 3 H 4 can be dominant in wide regions of the 5 The results for decay branching ratios are quite independent of the actual value used for the coupling, which determines u for fixed Z mass and sets the scale for the λ 1−6 couplings. parameter space, as long as θ 23 and θ 24 are not close to π/2. In particular, if these two angles are small, the mixing matrix is approximately with ε ij 0.01. Neglecting these small parameters, the Z couplings to the mass eigenstates are with sin β cos β, from which it can be clearly seen that the dominant Z decay to scalars is Z → H 3 H 4 . This is seen in the bottom right panel of figure 1, where we restrict |θ 23,24 | ≤ 0.01.
For the unrestricted scan (with θ 23 , θ 24 and θ 34 free), Br(H 3,4 → H 2 H 2 ) 1 in most of the parameter space of model 1, as can be seen in the left panel of figure 2 (for model 2 the results are similar). This is expected in the sense that these decays are controlled by the couplings λ 223 and λ 224 , respectively, which are not suppressed. The small mixing with ρ 0 allows the decay of H 2 into quarks, but has little effect on the decay of the heavier particles. 6 The same holds when |θ 23,24 | ≤ 0.01, as shown in the right panel of the same

Scenario 2
The masses are set to M Z = 3.3 TeV, M H 3 M H 4 400 GeV and M H 2 80 GeV as in the tagger benchmark labelled as std1500 in ref. [29]. The U(1) coupling is set to g Z z = 0.2. For the Z branching ratios the results obtained from the parameter space scan are the same as in scenario 1, and are omitted for brevity. This is so because the scalars are much lighter than the Z boson, and kinematical effects are unimportant. The decay Z → H 3 H 4 can dominate the scalar decays of the Z boson, in particular when θ 23 and θ 24 are small.
The results for the decay of the heaviest scalar H 4 are shown in figure 3 (left panel). For H 3 , which has nearly the same mass, the outcome is similar. In most of the parameter space Br(H 3,4 → H 2 H 2 ) 1, and the same happens for model 2.

Scenario 3
The masses are set to M Z = 3.3 TeV, M H 3 M H 4 400 GeV as in the tagger becnhmark labelled as std1500 in ref. [29] but, in contrast with scenario 2, M H 2 300 GeV in order to forbid the decay H 3,4 → H 2 H 2 . The coupling is set to g Z z = 0.2. The results of the parameter space scan are the same as in scenarios 1 and 2 for the Z branching ratios and, thus, they are not presented. Regarding H 4 decays, the results are shown in the right panel of figure 3 (for H 3 , with nearly the same mass, the outcome is similar). Both scalars decay into pairs of SM particles with branching ratios that are nearly independent of the mixing angles. The partial widths are determined by the matrix element O 1i and the small triple couplings λ 11i (i = 3, 4), as seen from eqs. (21).

Stealth boson signals
The potential relevance of stealth boson signals as a discovery channel is assessed in this section by a comparative study of the sensitivity of three searches, • Z → jj.
• A generic search, using the efficiencies for signals and background previously obtained for the anti-QCD tagger.
In addition, for scenario 1 we investigate whether the decay Z → H 3 H 4 is visible in a diboson resonance search. The various processes in our analysis are generated using Mad-Graph5 [38], followed by hadronisation and parton showering with Pythia 8 [39] and detector simulation using Delphes 3.4 [40] using the configuration for the CMS detector. The reconstruction of jets and their substructure analysis is done using FastJet [41]. For the signal processes the relevant Lagrangian is implemented in Feynrules [42] and interfaced to MadGraph5 using the universal Feynrules output [43]. As background processes we consider QCD dijet production pp → jj, with j being a light jet, pp → bb, and tt production. In order to populate with sufficient Monte Carlo statistics the entire mass and transverse momentum range under consideration, we split the samples in 100 GeV slices in the transverse momentum of the leading jet (or top quark), from 300 GeV to 2.2 TeV and above, generating 2 × 10 5 events for jj, 10 5 events for bb and 10 5 events for tt in each slice. The different samples are then recombined with weights proportional to the cross sections. Signal samples for Z → jj (including bb), Z → tt and Z → H 3 H 4 have 10 5 events each, except for Z → tt in scenarios 2 and 3 and Z → H 3 H 4 in scenario 3, which have 2 × 10 5 events.
The dijet and tt analyses are common to the two signal benchmark scenarios studied. We do not recast any specific experimental search but we choose event selections similar to the ones commonly adopted by the ATLAS and CMS Collaborations: • Dijet resonance analysis: jets are reconstructed with the anti-k T algorithm [44] using a radius R = 0.8, and groomed using Soft Drop [16] with the parameters z cut = 0.05 and β = 0. The use of large-radius jets is motivated by the possible presence of hard radiation accompanying the energetic decay products of a heavy resonance, and the grooming is implemented in order to clean the jets from pile-up and initial state radiation (see for example ref. [45]). The leading and subleading jets are required to have pseudo-rapidity |η| ≤ 2.5 and transverse momentum p T ≥ 500 GeV, while their pseudo-rapidity difference must satisfy |∆η| ≤ 1.1.
• tt resonance analysis: we use large-radius jets with R = 0.8 reconstructed and groomed as in the dijet analysis. The leading and subleading jets are required to have pseudo-rapidity |η| ≤ 2.5, |∆η| ≤ 1.1 and transverse momentum p T ≥ 500 GeV. For b tagging, a collection of 'track jets' of radius R = 0.2, reconstructed with tracks only, is used. A large-R jet is considered as b-tagged if a b-tagged track jet (using the 70% efficiency working point) within ∆R = 0.2 of its centre is found. The large dijet background is reduced by requiring that both leading and sub-leading jets are b-tagged, have a (groomed) mass 100 ≤ m J ≤ 220 GeV, and a value of the (ungroomed) subjettiness variable [7] τ 32 ≤ 0.7.

Scenario 1
In this scenario we have M Z = 2.2 TeV and we focus on the decay Z → H 3 H 4 , with M H 3 M H 4 80 GeV. For the signal coupling we choose g Z z = 0.15, yielding a total Z production cross section of 142.6 fb. 7 The total Z width is Γ = 26.5 GeV (model 1) and Γ = 50.2 GeV (model 2). The small Z width, compared to the experimental resolution, justifies using the same signal samples for both models. We set the mixing factor R 34 in eq. (14) to unity for simplicity -as seen in the previous section for the numerical examples provided, R 34 is often very close to one. Therefore, we obtain for the branching ratios It is expected that the tagger performance for bbbb, bbcc and cccc multi-pronged jets is similar, so we include both channels. With four scalars H 2 from the Z cascade decay, the branching ratio factor is Br(H 2 → bb, cc) 4 = 0.824.
The event selection for the generic and diboson analyes is the same as for the dijet search, but requiring groomed jet masses 40 ≤ m J ≤ 100 GeV, and • For the generic search we apply the tagger performance efficiency factors obtained in ref. [29] of 0.01 for QCD jets and 0.31 for the signal jets (M H 3,4 = 80 GeV, M H 2 = 30 GeV).
• For the diboson search we require τ 21 ≤ 0.4 for both jets. The mass window is wider than the usual ones for diboson resonance searches (for example, 65 ≤ m J ≤ 105 GeV is commonly used for W and Z jets by the CMS Collaboration) but we prefer to keep the same window as in the generic search for better comparison.
The presence of the heavy Z resonance can be detected as a bump in the dijet or tt invariant mass distribution. We present in figure 4 these distributions for the generic (top panels), tt (middle panels) and dijet (bottom panels) analyses. For model 1 we use an integrated luminosity of L = 15 fb −1 , while for model 2, for which the signal is much larger, we take L = 2 fb −1 . The signal and background cross sections for the different event selections are collected in table 3. Although the resonance is relatively narrow, the detector resolution effects yield a wider distribution, especially for the decays into scalars which produce four-pronged jets. As it has previously been shown [17,28], standard grooming algorithms are not adequate for multi-pronged jets, shifting jet masses and momenta from their original values. The effect can clearly be seen in the signal profile for the dijet analysis, which is much wider in model 2, where more than half of the dijet events are actually Z → H 3 H 4 .
The expected significance of the Z signal in the different searches is computed by using the Monte Carlo predictions for signal plus background as pseudo-data, performing likelihood tests for the presence of narrow resonances over the expected background, using the CL s method [46] with the asymptotic approximation of ref. [47], and computing the p-value corresponding to each hypothesis for the resonance mass. The probability density functions of the potential narrow resonance signals are Gaussians with centre M (i.e. the resonance mass probed) and standard deviation of 0.065M . The likelihood function is where i runs over the different bins with numbers of events n i , b i is the predicted number of background events and s i the predicted number of signal events in each bin, and µ a scale factor. We do not include any systematic uncertainty in the form of nuisance parameters. For each mass hypothesis the value µ b that maximises the likelihood function (37) is calculated, and local p-value is computed as with The results, assuming luminosities of 15 fb −1 (model 1) and 2 fb −1 (model 2) are presented in figure 5. As it can be seen, stealth boson signals in the generic search are by far more significant than standard dijet signals. In terms of standard deviations, the significance in generic searches is 4 (8) times larger for model 1 (model 2). As we have noted at the beginning of this section, the actual values depend on the product g Z z, which is a free parameter. Several additional comments and clarifications are in order.
• In our generic search, sensitive to stealth boson signals, we have focused on a jet mass window 40 ≤ m J ≤ 100 GeV, adequate for the benchmark point of the anti-QCD tagger considered. In order to cover all masses for the new scalars, experimental searches should explore bidimensional phase space, also using varying jet mass windows (see for example ref. [48]). • The use of b tagging in the generic search would significantly improve the significance for stealth boson signals. Requiring one b tag in either jet enhances the ratio S/ √ B by a factor of 2, and requiring two b tags by 3.4, where S stands for signal and B for background cross sections. We have chosen not to make use of b tagging in our analysis because the signals are already quite conspicuous, especially for model 2, and the background is already small. In this regard, our results are quite conservative. The use of b tagging would be very useful for large luminosities, to further reduce the background keeping the same signal efficiency for the anti-QCD jet tagger.
• We have not considered systematic uncertainties in our estimation of the significance of the different signals. These uncertainties will be more important in the channels where the background is larger, that is, jj and tt. In the generic search, where the expected background lies between 0.01 − 0.1 event per bin near the resonance mass, the impact of systematic uncertainties is expected to be small.
• As aforementioned, existing grooming algorithms are not designed nor optimised for multi-pronged jets and may shift the mass peaks. (Several other grooming algorithms and parameters were explored in ref. [17] with similar results.) This is clearly observed in figure 5, where the maximum significance is for the stealth boson signal is near 2 TeV while the Z mass is 2.2 TeV. We have not attempted any mass recalibration because this small shift does not affect our results and conclusions.
• The relative (in)significance of the Z signal in dijet, tt and diboson searches depends on the model and the mass of the lightest scalar, as the signal for the dijet and diboson event selections receive contributions from various Z decay modes.
• For all the final states considered, the global significances of the deviations -which depend on the mass range studied in each experimental search -are smaller than the local significances in figure 5. In any case, for the point addressed in this section, namely to show that the sensitivity of a generic search is a factor of ten (in terms of significance) or more than for current searches, local significances suffice.
• For lighter H 2 , the four-pronged stealth boson jets have a more two-pronged structure, and the acceptance in diboson searches is slightly larger (see appendix C). Actually, for M H 2 = 30 GeV most of the signal that passes the diboson event selection is Z → jj and Z → tt, not Z → H 3 H 4 .
For completeness, let us comment about the direct production of the new scalars H 2−4 , which can be produced in the same processes (gluon-gluon fusion, associated production with a W/Z boson, etc.) as the SM Higgs H. The cross sections are the ones that would correspond to a SM Higgs of the same mass M H 2−4 , multiplied by the small factor O 2 1i . In the benchmark scenarios considered, with mixing angles |θ 1j | ≤ 0.01 (j = 2, 3, 4), all cross sections for processes mediated by SM particles are suppressed by a factor of O 2 1i 10 −4 , leading to unobservable signals.
The most stringent constraints on the lightest scalar H 2 result from the Higgs boson searches at LEP experiments [49]. For a mass of 30  In our benchmark, σ(H 3,4 ) × Br(bbµ + µ − ) 0.67 fb, three orders of magnitude smaller. At the Tevatron, the CDF Collaboration performed a search for pair production of new particles Y , each decaying into two jets, pp → Y Y → jjjj. The mass range explored M Y ≥ 50 GeV does not cover M H 2 = 30 GeV, but for illustration we can take the limit for M Y = 50 GeV, namely σ(Y Y → jjjj) ≤ 200 pb. In our benchmarks, the prediction is σ(H 3,4 ) × Br(jjjj) 0.1 fb, four orders of magnitude below that limit.  Table 4: Signal and background cross sections for the Z signal in scenario 2 and main SM backgrounds (in rows) under the three different event selections for generic, dijet, and tt resonance searches. The event selection for dijet and tt is the same as in scenario 1, and the quoted background numbers are the same as in table 3.

Scenario 2
This scenario is similar to scenario 1 but with heavier masses, M Z = 3.
The combined branching ratio factor for the four H 2 decays into quark pairs is Br(H 2 → bb, cc) 4 = 0.843. The event selection for the generic analysis is the same as for the dijet search, but requiring groomed jet masses m J ≥ 250 GeV. We apply the tagger performance efficiency factors obtained in ref. [29] of 0.01 for QCD jets and 0.33 for the signal jets (M H 3,4 = 400 GeV, M H 2 = 80 GeV).
The dijet / tt invariant mass distributions are presented in figure 6 for the generic (top), tt (middle) and dijet (bottom) analyses. Given that the cross sections for Z production are smaller to those of scenario 1, we present our results for integrated luminosities L   The difference between stealth boson modes and standard tt and dijet decays is quite pronounced. The significance of the signals in the generic search (expressed in terms of standard deviations) is 3 and 6 times larger than in dijets, for model 1 and model 2, respectively. Still, we remind the reader that we have not taken advantage of b tagging, which would improve the signal significance by a factor of 2 − 5 in the generic search.
Let us also comment on signals of the direct production of the new scalars within this scenario. Direct production of the lightest scalar H 2 with decay into SM particles (with its decay branching ratios corresponding to a SM Higgs with a mass of 80 GeV) can be constrained from Higgs boson searches. At LEP, the non-observation of a signal constrains O 2 1i ≤ 0.04 for M H 2 = 80 GeV [49], two orders of magnitude above the bound O 2 1i 10 −4 in our benchmark. At the Tevatron, searches by the D0 and CDF Collaborations [54,55] only cover masses above 90 GeV, and also are less sensitive.
For the heaviest scalars the production and decay chain is pp → H 3,4 → H 2 H 2 , with subsequent decay of H 2 . A search for pair produced resonances decaying into quark pairs by the CMS Collaboration [56] is sensitive to (but no optimised for) this process. The limits corresponding to stop masses mt = 80 GeV are σ(tt * ) ≤ 400 pb, assuming 100% branching ratio for the R-parity violating decayt → bq. In our benchmark, the cross section times branching ratio for final states with b quarks is σ(H 3,4 ) × Br(bbbb) 0.20 fb for M H 3,4 = 400 GeV. This is six orders of magnitude below the above experimental limit.  Table 5: Signal and background cross sections for the Z signal in scenario 3 and main SM backgrounds (in rows) under the three different event selections for generic, dijet, and tt resonance searches. The event selection is the same as in scenario 2, and the quoted backgrounds are the same as in table 4. The signal in the tt selection is also the same.

Scenario 3
Here we take M Z = 3.3 TeV and M H 3 M H 4 400 GeV, as in scenario 2, and we keep the same signal coupling g Z z = 0.2. Therefore, the Z production cross section and width are the same, σ(Z ) = 20.1 fb, and Γ = 70.2 GeV in model 1, Γ = 127.7 GeV in model 2.
We set R 34 to unity, hence the branching ratios are the same as in scenario 2, The differences with respect to scenario 2 stem from the fact that H 2 is now heavier, namely M H 2 300 GeV. This forbids the decays H 3,4 → H 2 H 2 . We therefore focus on H 2 decays into gauge boson pairs, taking (see The event selection for the three analyses is the same as in scenario 2 (for the dijet and tt searches it is also the same as in scenario 1). The dijet invariant mass distributions are presented in figure 8    is 2.5 and 9 times larger than in dijets, for model 1 and model 2, respectively. Regarding the results for the generic analysis, it is worth remarking a few points.
• The requirement of jet masses m J ≥ 250 GeV in the generic analysis filters the hadronic decays of both gauge bosons, which have branching ratio of 0.45 for W W and 0.49 for ZZ, reducing the signal with respect to scenario 2. This explains why the expected sensitivity in the generic analysis is slightly worse, despite the larger efficiency of the tagger for jets corresponding to H 3,4 → W W (0.49 for the signal for a background rejection of 10 2 ) than for jets with H 3,4 → H 2 H 2 in scenario 2.
• The decays H 3,4 → HH have a branching ratio of 0.13, and the generic search would also be sensitive to jets containing two SM Higgs bosons. Because there is no benchmark working point of the tagger available, we have not included these signal contributions in our analysis.
• The tagger has an efficiency of 0.21 for jets containing two top quarks, resulting from H 3,4 → tt. For simplicity we have not included this signal contribution to the generic search, given the small branching ratio Br(H 3,4 → tt) = 0.04.
Therefore, although the significance of the signals simulated for this scenario is smaller than in scenario 2, it is expected that the significances become similar when all possible heavy scalar decay channels are included. Let us also comment on the direct production of the scalars. Searches for new scalars produced in gluon gluon fusion and decaying into V V pairs set a limit of approximately σ × Br(V V ) ≤ 350 fb for a scalar mass of 400 GeV. In this benchmark we have σ(H 3,4 ) × Br(V V ) ≤ 0.26 fb, more than two orders of magnitude below the limit. The lightest new scalar H 2 with M H 2 = 300 GeV is not covered by this analysis, which only considers masses above 400 GeV. Still, σ(H 2 ) × Br(V V ) ≤ 0.66 fb is well below potential constraints at this mass.

Discussion
The search for elusive new physics signals yielding various types of multi-pronged jets requires a model-independent approach, with the use of novel tools like the anti-QCD tagger [29] or non-supervised learning methods [58][59][60]. In order to contextualise the relevance of these signals as discovery channel for new (leptophobic) resonances, it is crucial to provide examples of consistent models that may produce them. We have done so for stealth bosons, boosted particles with a cascade decay giving a four-pronged fat jet. We have worked out the minimal implementation, adding to the SM a leptophobic Z boson, two complex scalar singlets and extra matter, either new vector-like quarks (model 1) or new vector-like leptons (model 2), to cancel anomalies. In these models one can compare the potential significance of stealth boson signals, still unexplored at the LHC, with the standard signals (dijets, top pairs and dibosons) already searched for. Depending on the model and benchmark scenario considered, the significance of the former may be up to 9 times larger than the most sensitive of the latter. Therefore, it is clear that stealth boson signals might well be hidden in LHC data, yet invisible to current searches. Besides, direct production of the new light scalars is suppressed by the square of the small mixing, and signals are too small to be observed.
In the two models considered in this work the branching ratios of Z decays into scalars are sizeable (around 10% in model 1 and 50% in model 2). Moreover, cascade decays of the new scalars are likely to happen, provided one of the following conditions are fulfilled: • There is a hierarchy among the masses of the new scalars, so that the decays of one into others are possible. These decays are not suppressed by mixing with the SM scalar doublet, and will therefore be dominant, as in our scenario 1.
• The scalars are heavy enough to decay into W + W − (and possibly ZZ, HH and tt).
If the decay into other new scalars is kinematically allowed, it will be the dominant channel, as in our scenario 2. Otherwise, decays into pairs of SM bosons will be dominant, as in our scenario 3.
Therefore, as it has been shown with a scan on parameter space, it is natural to have stealth bosons as decay products of the Z . For simplicity, we have restricted our detailed simulations to Z decays into a pair of stealth bosons giving two four-pronged jets. Still, those processes giving one stealth boson (four-pronged jet) and one scalar that subsequently decays into quarks (two-pronged jet) are also possible and interesting. A generic search would be sensitive to all these possibilities at once, and this is one of the main virtues of the anti-QCD tagger.
In conclusion, we stress that despite the fact that stealth bosons are rather stealth for current LHC searches, they would be quite conspicuous in a generic search. Moreover, these signals may well appear in decays of heavy Z resonances. These facts already provide a strong motivation for model-independent searches.

A Triple scalar couplings
In the weak interaction basis H i = (ρ 0 ρ 1 ρ 2 A 0 ) the trilinear scalar interactions can be expressed in the condensed form: with the sum over p ≤ q ≤ r running from from 1 to 4. The C ijk coefficients are explicitly given by: [ λ 5 cos β + Re(λ 9 ) sin β ] , [ Re(λ 9 ) cos β + λ 6 sin β ] , Im(λ 9 ) cos β , Im(λ 9 ) sin β , In the H a mass eigenstate basis, with where the sums over a, b, c and p ≤ q ≤ r run from 1 to 4.
When two of the indices i, j, k are equal, the sum (48) contains each of the three independent permutations twice, thus introducing a double counting. When the three indices are equal, i = j = k, this sum counts six times the single term H i H i H i present in the sum (46). One can take this fact into account by introducing a symmetry factor S ijk , which is one if the three indices are different, two if two of the indices are equal, and six if i = j = k. With this convention, the interaction is (no sum over indices) keeping the definition (47) for λ ijk and all permutations (48), even repeated ones. When deriving the Feynman rule for the three-scalar interaction, one has to multiply by a symmetry factor for the presence of identical particles, which is precisely S ijk . Therefore, the Feynman rule for the vertex is simply −iuλ ijk .

B Model with two scalar doublets and one singlet
An attractive SM extension which apparently could lead to stealth boson decays would be that with and extra scalar doublet and a scalar singlet. However, this model does not produce any of the desired processes For illustration and completeness we summarise why. We label the two existent scalar doublets as Φ 1 = (φ + 1 φ 0 1 ) T and Φ 2 = (φ + 2 φ 0 2 ) T , and the singlet as χ. The scalar potential compatible with the SU(2) L × U(1) Y symmetry is Among the parameters above, m 11 , m 22 , m 0 , λ 1−4 and λ 8−10 are real, while m 12 , λ 5−7 and λ 11 can be complex. Writing the neutral scalar fields in the usual way: the would-be Goldstone bosons associated to the breaking of the U(1) and electroweak symmetries are η 3 , and a combination of η 1 and η 2 (as in the usual two-Higgs doublet model), respectively. Therefore, we have three scalars and a pseudo-scalar, which can in principle mix.
At least one of the two doublets must have vanishing hypercharge Y , as required by the existence of Yukawa terms (3). We choose it to be Φ 1 . If the other doublet Φ 2 has hypercharge Y Φ 2 = 0, then invariance under U(1) requires m 12 = 0, λ 5−7 = 0, λ 11 = 0. After applying the potential minimisation conditions it is found that the physical pseudoscalar is massless, which is unacceptable. Besides this obvious drawback, we note that the vacuum expectation value φ 0 2 = v 2 / √ 2 contributes to Z − Z mixing [34], which is constrained to be very small. Because the Z − Z coupling to the scalars in Φ 2 is proportional to v 2 , the width for Z → H i Z, which would be characteristic for this model, is also very small. If both doublets have Y = 0, neither of the decays in (50) is present, the former because of the vanishing doublet hypercharges, and the latter because the only Z coupling to scalars is Y χ Z µ η 3 ← → ∂ µ ρ 3 and η 3 is not physical.

C Alternative scenario 1
We present here results for an alternative scenario 1 for model 2, with M H 2 = 15 GeV, in which the substructure of the four-pronged fat jets resulting from H 3,4 → H 2 H 2 → qqqq resembles more a two-pronged structure because of the lighter H 2 . For this mass, we have Br(H 2 → bb) = 0.81 , Br(H 2 → cc) = 0.12 .
so that Br(H 2 → bb, cc) 4 = 0.732, being the signals slightly smaller. The efficiency of the tagger for this signal is practically the same as in scenario 1. The dijet mass distribution for the diboson analysis is shown in figure 10 for M H 2 = 30 GeV (left panel ) and M H 2 = 15 GeV (right panel). Besides the large differences in the cross section, due to the larger acceptance for Z → H 3 H 4 , in the latter case we observe a resonant signal structure that is not present in the former. The p-value for the different Z signals is given in figure 11. Notice that, despite the fact that a possible signal would be more visible in the diboson resonance searches than in tt and dijet final states, it is by far surpassed by the signal that would be visible in a generic search using the anti-QCD tagger.