Multi-track Displaced Vertices at B-Factories

We propose a program at B-factories of inclusive, multi-track displaced vertex searches, which are expected to be low background and give excellent sensitivity to non-minimal hidden sectors. Multi-particle hidden sectors often include long-lived particles (LLPs) which result from approximate symmetries, and we classify the possible decays of GeV-scale LLPs in an effective field theory framework. Considering several LLP production modes, including dark photons and dark Higgs bosons, we study the sensitivity of LLP searches with different number of displaced vertices per event and track requirements per displaced vertex, showing that inclusive searches can have sensitivity to a large range of hidden sector models that are otherwise unconstrained by current or planned searches.


Introduction
The nature of dark matter (DM) remains one of the most important outstanding questions in particle physics. There has been a recent shift in attention towards DM masses below the electroweak scale due to the relative lack of constraints compared to weakly interacting massive particles (WIMPs) and the promising prospects for new experimental studies [1][2][3][4]. For masses below about 5 GeV, it has long been established that thermal DM models require the existence of new mediators between DM and the Standard Model (SM) [5]. The reason is that the SM mediators in WIMP models, such as the Z and Higgs bosons, are much heavier than the DM mass and consequently DM annihilation is sufficiently feeble that the DM abundance cannot be depleted enough to match observations. Viable GeV-scale DM models necessitate the addition of new mediators [6], such as a dark photon [7,8] or dark Higgs [9][10][11][12], which appear in models where DM is charged under a force or its mass originates from spontaneous symmetry breaking. All of this suggests an entire new hidden sector of particles coupled to DM: not only does this make the physics of DM more analogous to that of the SM, but hidden-sector dynamics can also have other consequences, including the generation of SM neutrino masses, the creation of a baryon asymmetry, or self-interactions that affect DM structure formation.
Most searches for hidden-sector particles at intensity-frontier experiments have focused on a relatively small subset of minimal models. In most of these scenarios, there is only one new particle that interacts with the SM, and both its production from and decay to SM states is determined by a single coupling. Alternatively, this single particle can be treated as completely invisible, either because it is meta-stable or because it decays to invisible, stable DM particles. This approach allows for powerful searches for new hidden-sector particles by exploiting the very specific production and decay modes predicted by the model, and also permits the ready comparison of different experiments using a common set of benchmarks, such as those proposed by the Physics Beyond Colliders working group [4].
However, as soon as additional states are added to the hidden sector, the phenomenology can change dramatically (see, e.g., [2,13]). For example, if we consider the case of a dark photon that is kinetically mixed with the SM photon, there are strong constraints on the model if the dark photon decays either exclusively through its coupling to the SM [14][15][16][17][18], or invisibly to unseen hidden-sector states [19][20][21]. If, however, the dark photon decays to hidden-sector states that partially decay back to the SM (for example, in the case of inelastic DM [22][23][24][25][26][27][28][29][30][31]), then neither search strategy applies, and the parameter space can be relatively unconstrained. Indeed, truly model-independent constraints from electroweak precision tests [32] and deep inelastic scattering [33] are sufficiently weak that they do not completely exclude dark photons as an explanation for the 4σ discrepancy between theory and experiment in (g − 2) µ [34,35]. If the mediator particle decays in a partially visible manner, then we must change our search strategies: the mediator simply can't be discovered in the usual channels! One of the major challenges of probing multi-particle hidden sectors is the proliferation of signatures: if we try to constrain the hidden sector with exclusive, highly targeted searches, then a small variation in the nature of the hidden sector invalidates the searches. It is an open problem to determine the best way to probe the model space in a manner that covers many different scenarios while limiting the number of experimental searches. At the Large Hadron Collider (LHC), the standard approach is to perform relatively inclusive searches, with results presented in a simplified model framework that facilitates re-interpretation [36][37][38][39]. These are supplemented by targeted, exclusive searches for well-motivated signals that are otherwise limited by large backgrounds. By contrast, at B-factories, which offer relatively clean environments and high luminosities that are ideal for discovering hidden-sector particles with masses below 10 GeV, most hidden-sector searches exploit exclusive signatures that target only a single, specific hidden-sector model [40][41][42][43][44][45][46][47][48]. The B-factories' hermetic nature, and excellent particle identification and reconstruction capabilities, offer interesting parallels with high-energy collider experiments such as ATLAS and CMS, making B-factories promising sites for generic searches targeting multi-particle hidden sectors.
In this paper, we advocate for a program of inclusive searches at B-factories that can dramatically improve the sensitivity of existing experiments to multi-particle hidden sectors. We focus on signatures with long-lived particles (LLPs) for a number of reasons: LLPs are ubiquitous in hidden-sector models, and backgrounds are expected to be small. Furthermore, only a limited number of these searches currently exist at B-factories [42,47,49], and because a significant amount of work has been done at the LHC to generalize and expand the scope of LLP searches [39,50], much of this knowledge gained can be applied to B-factories as well. We propose several extensions of existing searches for LLPs at BABAR [49] and LHCb [51], both of which are inclusive with respect to production mechanism but require exclusive reconstruction of 2-body decays: our proposals include searching for displaced vertices (DVs) with more than two tracks or two DVs in an event. Because of these more stringent selection criteria, such searches can relax other aspects of the analysis, such as prompt objects associated with the specific LLP production mechanism, or the exclusive reconstruction of the LLP mass; this allows a significant broadening of the applicability of the searches. We find that the projected sensitivities of our proposed searches can surpass the sensitivities for minimal dark photons or dark Higgs bosons. Our work also complements recent proposals for new, model-specific LLP searches at B-factories [29,30,[52][53][54].
To motivate our approach to inclusive searches, we use SM hadronic physics as an analogy for the types of physics expected in the hidden sector [13]. Hadrons are produced copiously in pairs via the strong or electromagnetic interactions, but some can only decay through the weak force due to an approximate flavour symmetry, leading to a macroscopic lifetime due to a small ratio between the LLP and W masses. In the case of SM LLPs, the decay of the LLP is well-described by an effective field theory (EFT), namely the Fermi theory [55], while LLP production occurs through a low-mass mediator such as a photon, gluon, or pion.
We propose an analogous EFT approach for the study of hidden-sector LLPs. In particular, we explore several pair-production mechanisms for new LLPs, including dark photons and dark Higgs bosons. The decay of the LLP is in turn dictated by an EFT operator that couples a single LLP to SM states. This approach has several advantages: it systematically captures all hidden-sector scenarios where the LLP decays through a heavy, off-shell mediator back to SM states; the decay width of the lightest hidden-sector state is suppressed by powers of the ratio of its mass to the EFT cutoff scale, naturally leading to long lifetimes; and it generates a rich array of phenomenological signatures. We study a number of production and EFT decay modes to show that our minimal set of inclusive LLP searches can cover a wide array of possible signals. Although our EFT framework technically only applies to off-shell degrees-of-freedom, it also characterizes the gauge-invariant SM final states that can be produced from on-shell intermediate states.
Our paper is organized as follows: in Sec. 2, we introduce our EFT framework for hiddensector particle decays and argue that LLPs generically arise in such a scenario. We summarize our key results in Sec. 3, providing details of the particular production and decay modes in Sec. 4 and simulation details and quantitative results in Sec. 5. Although our results are based on the assumption of approximately background-free signal regions, we consider the leading backgrounds and mitigation strategies in Sec. 6, and our concluding discussion is found in Sec. 7.

Long-Lived Particles from a Hidden Sector
We consider the set of models where the hidden sector is composed of a mediator, φ, and an LLP, χ; φ and χ could have any spin consistent with the decay φ →χχ. We assume that both have masses 10 GeV so that they are accessible for on-shell production at B-factory energies, e + e − → φ + X, φ → χχ, where X is some unspecified SM states that could be produced in association with φ. Apart from the mediator, φ, the only other interactions between the SM and the hidden sector occur at a much higher energy scale, Λ, which we take to be the scale of new physics that mediates χ decays.
A good way to parametrize LLP decays is through effective operators [38,56]. 1 The reason is that the combination of mass ratios and couplings in LLP decays are typically tiny, so as long as the particles mediating the decay are at all heavier than the LLP, the EFT approach is reasonable. The operators mediating LLP decay are of then of the form χO SM ; O SM must be a singlet under U (1) EM and SU (3) C , and it is constrained by the spin of the LLP χ, but is otherwise undetermined. 2 We will argue below that this leads to the observation that χ typically decays to many electrically charged particles in the final state. dim χO SM leptonic semi-leptonic hadronic photonic The LLP can be of any spin, but we will exclusively focus on spin-0 and spin-1/2 LLPs because these cover most of the phenomenologically relevant signatures. 3 This set-up enables a quasi-systematic way of outlining the possible LLP decay modes. We simply list the possible operators O SM subject to the gauge invariance constraints and organize them by operator dimension. The latter is key in determining the resulting LLP lifetime τ χ ≡ Γ −1 χ , where n = dim χO SM , Λ is the scale suppressing the higher-dimensional operator in the effective action, and we dropped the phase space factor which depends on the precise final-state multiplicity. If we consider a width of Γ ∼ 10 −16 GeV (corresponding to a proper decay length of about 1 m, which is the characteristic length scale of the B-factory detectors) and m χ ∼ 1 GeV, then we have Λ ∼ 10 16/(2n−8) GeV.
For specific operator dimensions, we find: for n = 5, Λ ∼ 10 8 GeV; for n = 6, Λ ∼ 10 4 GeV; for n = 7, Λ ∼ 500 GeV; for n = 8, Λ ∼ 100 GeV. Thus, for n ≤ 6, it is very easy to get lifetimes spanning essentially the whole range from 1 mm to 1 m while ensuring the EFT cutoffs are above existing or prospective LHC limits. For n = 7, it is straightforward to get lifetimes ∼ 1 m, but shorter lifetimes may come into conflict with LHC measurements (particularly if the off-shell particles at scale ∼ Λ have quantum chromodynamics, or QCD, charge). Finally, n ≥ 8 always gives LLPs but the lifetime may be too long to detect at B-factories while simultaneously satisfying constraints from high-energy colliders. There is still a small acceptance for LLP decays inside the detector even for cτ m; in this case, missing-momentum searches also constrain the model, and these may or may not have better sensitivity than a displaced vertex search.
The task is now straightforward: list all SM gauge singlet operators O SM with where dim χ = 1, 3/2, and then classify the resulting χ products. We neglect operators involving electroweak bosons since these cannot be produced on-shell in B-factories, and the additional (m χ /m W ) 4 contribution to the χ decay rate from such operators pushes the lifetime to be well beyond that accessible at B-factories. Since O SM is a gauge singlet, one can find all independent 3 Spin-1 particles can be long-lived as well, including a weakly coupled fundamental gauge boson (such as a long-lived dark photon) and composite states like the ρ mesons of a confining hidden sector [59], but these typically couple to conserved currents leading to two-body decays. One exception to this is a dark photon lighter than 2me, which decays to 3 SM photons [60]. choices of this operator by harnessing the results of SMEFT (the collection of higher-dimensional operators constructed from SM fields alone) [61][62][63][64] if the LLP is a scalar, or the SM + sterile neutrino EFT (νSMEFT) [65][66][67] if the LLP is a fermion. The corresponding operators are summarized in Appendix A. In Tables 1 and 2, we collect the parton-level final states from χ decay enabled by these operators for spin-0 and spin-1/2 LLPs, respectively.
One pronounced feature of final states from the EFT framework is that the majority feature multiple charged particles. Note that when the operator involves quarks, the actual decay final states are generated by hadronization which may lead to a different number of charged particles than naïvely expected depending on the LLP mass. We use the above observation to develop a generic set of analysis proposals for multi-track displaced vertex searches that can probe many of these final states.

Summary of Main Results
Before turning to our detailed results, we provide a brief summary of our primary conclusions. We propose inclusive displaced vertex searches, studying the signal acceptance and efficiency as a function of the number of displaced vertices per event, N dv , and the number of tracks required per vertex, N tr . We expect that searches with N dv ≥ 2 and N tr ≥ 3 should be background free, although searches with lower track and/or vertex multiplicities are possible especially if additional selections (like vertex track mass cuts) are implemented. Our findings include: • A minimal set of searches has broad coverage of different production modes, decay modes, and LLP lifetimes. While the detailed acceptances depend on model specifics, all of the simplified production and decay modes we consider have appreciable acceptance for LLP searches with different N dv and N tr , including N dv ≥ 2 and N tr ≥ 3. In most cases, the sensitivity is sufficient to probe interesting parameter space motivated by (g − 2) µ or dark matter models. 4 Expected coupling sensitivities are competitive with or surpass existing searches for fully visible or invisible dark photons and dark Higgs bosons.
• Signal acceptance falls with stricter requirements on N dv and/or N tr , but not as severely as we expect backgrounds to be suppressed. For combinatoric and other rare backgrounds, increasing the track multiplicity, vertex multiplicity, track mass, etc., can effectively suppress the backgrounds. Meanwhile, the signal acceptance is always somewhat reduced by the tighter selections, but usually the residual acceptance is still sufficient to probe interesting parameter space. 4 While the (g − 2)µ anomaly can be explained by new particles, recent lattice calculations of hadronic vacuum polarizations in the SM alone can also bring the theoretical prediction into agreement with the BNL and FNAL results [68].
• LLPs that decay to partially or fully hadronic final states can give large multiplicities of charged tracks in the final state. Even when the quark-level operators suggest a low multiplicity of charged final states, hadronization effects (whether in a parton shower or chiral perturbation theory) can give higher multiplicities of charged pions or kaons in the final state. This allows for a significant signal acceptance even with tighter cuts N tr ≥ 3. Heavier LLPs more readily pass selections requiring more than two tracks per vertex due to more pronounced effects of parton showering or heavy intermediate QCD resonances, while lighter LLPs can suffer a severe acceptance penalty when requiring more than two tracks per vertex for hadronic decays.
• Hadronization uncertainties can be significant depending on the model, and are greatest for low-mass LLPs. When LLP decay operators resemble those in semileptonic decays of SM particles, the effects of hadronization are relatively well-understood. For other operators, the uncertainties are larger, especially for LLP masses such that the typical hadronic 4-momentum invariant is below 1 GeV. Even when different models of hadronization disagree, we still find a reasonable signal acceptance for all hadronization models that allows DV searches to probe interesting parameter space. Crucially, the search itself can be done without worrying about theory uncertainties, as it is only the mapping of the results of a search to a particular model that is sensitive to such uncertainties. For more details, see Appendix B.
• Exclusive signal reconstruction would reduce or eliminate the sensitivity to many of the models we consider. However, some model-dependent selections could help suppress backgrounds if needed. For several of our models, it is challenging to fully reconstruct the LLP decay and/or production signal, either due to invisible particles or the fact that there are many possible LLP decay modes. For examples, many modes include π 0 → γγ in the LLP decay chain that could be difficult to associate with the displaced vertex. These models are not covered by the existing model-independent search by BABAR which requires a reconstruction of the LLP mass in displaced track pairs [49]. However, if a totally inclusive search is limited by backgrounds, it is possible to improve sensitivity by applying a few additional selections that depend on the model. These could include particle identification requirements on one or more of the tracks at the displaced vertex, tagging a photon from initial-state radiation (ISR) associated with dark photon production in the radiative return process, and so on. We encourage experimentalists to apply as few of these requirements as is feasible to ensure the broadest coverage of any given search.

Production and Decay of Long-Lived Particles
The central observation of the previous sections is that many operators connecting the hidden sector to the SM induce decays to several charged particles, giving rise to multi-track events. Thus, it is useful to develop a widely applicable search strategy to target these generic final states. In this section, we identify several test cases that lead to events with multiple tracks emanating from displaced vertices and that can be probed with this approach. As a reminder we denote the mediator by φ and the LLP by χ, although these can each be scalars, fermions, or vectors depending on the interaction. We consider leptonic, semi-leptonic, and hadronic final states. While these final states arise from a multitude of operators, for simplicity of simulation we focus on the following operators for concreteness: • leptonic final states (scalar LLP) χμµēe; • fully and non-fully reconstructible semi-leptonic final states (fermion or scalar LLP) • hadronic final states (scalar LLP) We note that the specific flavour combinations chosen are somewhat arbitrary; we selected these either to allow for the exploration of the widest kinematic range in LLP decays by focusing on light flavours, or because some hadronic decays can be modelled by analogy with certain SM processes. We implement these models in the Universal FeynRules Output (UFO) format [69] using FeynRules [70]. These models are then used in MadGraph [71,72] to simulate the relevant processes. The operators containing quarks require special attention since they must be hadronized to generate physical final states. We describe our treatment in detail in Appendix B.
We must also specify a production mechanism for the LLPs, which can be parametrized by a few benchmark simplified models. For simplicity, we study two production modes: 1. Photon associated production, namely e + e − → γφ, φ → χχ. This occurs for vectorportal models through the radiative return process (either through kinetic mixing with the SM photon or direct coupling of the vector to electrons) [15,19], or scalars with direct couplings to electrons. It is also valid for axion-like particles that possess a φγγ coupling (and the production proceeds through e + e − → γ * → γφ), although the kinematics is distinct from radiative return. A common feature to both scenarios is the monochromatic photon that can, in principle, be used to discriminate against backgrounds with a bump hunt in the photon energy, provided it is tagged. The multiplicity of tracks and photons from the production process is low, which could allow us for discrimination against QCD processes. For concreteness, we use a dark photon model to implement this production mechanism (see Sec. 4.1 for more details).
2. Heavy-flavour production. Mediators lighter than the B meson can be produced in flavour-changing neutral current decays such as B → Kφ, φ → χχ. If necessary, it is possible to suppress continuum QCD backgrounds using event-shape variables and by requiring that the event (with B → Kφ candidate) matches the constrained kinematics of Υ(4S) → BB. This production mode is also interesting because there should be a double displaced vertex (one "vertex" from the K production and another from the LLP decay). These are generally high track-multiplicity events, which can result in larger backgrounds. For concreteness, we use a dark Higgs model to implement this production mechanism (see Sec. 4.2 for more details).
In both cases, we assume 100% decay of the mediator into LLPs. These production modes are well-motivated and also encompass a range of mediator and LLP kinematics.
We emphasize that our proposed searches are designed to be inclusive of the production mechanism. As a result, sufficiently inclusive searches can actually be sensitive to broader classes of scenarios, such as: • Heavy-flavour leptonic forces: If the mediator is an L µ − L τ gauge boson [73,74] or a leptophilic scalar [75,76], the dominant production of the LLPs will be in association with µ + µ − or τ + τ − pairs. This channel was recently studied for dielectron LLP decays in the BABAR leptophilic dark scalar search [47].
• Dark-flavour non-diagonal production: As an example, consider inelastic dark matter coupled to a dark photon, e + e − → γA , A → χ 2 χ 1 , where χ 1 is a stable invisible particle and χ 2 is the LLP [24][25][26][27][28][29][30][31]77]. This example is an illustrative model that predicts only a single LLP, showing the importance of single-LLP searches (N dv = 1) where possible. In this minimal model, the χ 2 exclusively decays through the same dark photon coupling that allows its production, but more general models could additionally include different decays of the χ 2 .
• Neutrino production: If a neutral fermion N mixes with the SM neutrino, then N is produced in meson decays such as B → N + X [78] (where X usually includes a D meson, but could also be "nothing", in which case we have exclusive production in leptonic decays). If N decays through the same neutrino mixing, then we recover the Majorana neutrino scenario that is the subject of a search at Belle [42], although that search only examined a single exclusive decay mode for the LLP, N → ± π ∓ . If the model additionally contains a singlet scalar S, the dominant decay could instead be N → Sχ, where either S or χ could be invisible stable particles or LLPs.
• Strongly interacting hidden sectors: An example of this is a confining hidden valley [13] where a dark photon decays into hidden-sector quarks that then shower and hadronize [79]. This could encompass Strongly Interacting Massive Particles (SIMPs) [80,81] In several of the above examples, the dark sector goes beyond the "next-to-minimal" set-up we consider (a mediator, LLP and a higher dimensional decay operator), and although each merits its own careful study, it is evident that inclusive LLP searches at B-factories should encompass a variety of signals in terms of number of charged particles per LLP decay and vertex multiplicity.

LLP Production: Radiative Return (Dark Photon)
The main mechanism by a which an electron-coupled vector or scalar mediator, such as a dark photon, is produced at B-factories is via radiative return. In the radiative-return process, ISR photons drive the e + e − center-of-mass (CM) energy to resonance. We will use the dark photon, A , with the kinetic mixing interaction (ε/2)F µν F µν as the specific mediator example for our implementation of radiative return. The corresponding Lagrangian that enables A production and decay into LLPs is where J µ EM (J µ χ ) is the electromagnetic (dark) current, and g D is the coupling between A and the LLP.
The total radiative return cross section is, in the narrow-width approximation (see, e.g., Ref. [82]), where In the above expression, the first factor is the resonant annihilation cross section for e + e − (divided by the Mandelstam variable s and with the delta function stripped off), while the second factor is the effective e + e − parton luminosity with f e being the electron distribution function which encodes the probability of finding an electron or positron with momentum fraction x. Note that Eq. (8) includes the part of phase space in which the photon is collinear with the beams and therefore unobservable; this cross section can be a factor of ln(s/m 2 e ) ∼ 20 larger than if a photon is required in the detector acceptance (a requirement used in existing dark photon searches at BABAR [15,19]). We show the radiative return cross section as a function of A mass in Fig. 1.
To simulate events with radiative return production of A , we sample x, the lepton momentum fraction, from the parton luminosity function using the Kuraev-Fadin distribution [83] from Ref. [84]. The kinematics of A and its decay products in the lab frame can then be reconstructed by noting that A is produced nearly at rest in the CM frame in the narrow-width approximation. We can then boost the combined decay chain of A to LLPs and LLPs to SM (generated using MadGraph) into the lab frame. While a publicly-available lepton ISR plugin exists that would have allowed us to include ISR directly in MadGraph [85,86], we found that it had difficulty populating the phase space around the narrow A resonance. We have verified that our procedure gives the correct parton luminosity, proportional to the integral piece in Eq. (8), by matching to the results of Ref. [82,87]. We also compared the results of Eq. (8) with MadGraph including ISR, finding that the cross sections match at those parameter points for which MadGraph events were successfully generated.
For m A < √ s the expected number of LLP pairs produced for the full Belle II integrated luminosity is where we normalized the rate to the kinetic mixing that saturates model-independent bounds from precision electroweak observables [32]; for m A ≈ 2 GeV, ε ≈ 0.03 also resolves (g − 2) µ . While for minimal assumptions (purely visible or invisible decays of A ) these "large" kinetic mixings are excluded by exclusive searches [15,19], this parameter space can be open if the A decays do not satisfy the selection criteria in those analyses.
The large number of allowed LLP events in Eq. (9) then suggests that novel analyses can offer significant discovery potential (including for (g − 2) µ -preferred parameter space) even with small signal acceptance. To illustrate this, we show in Fig. 2 the schematic background-free discovery potential in the dark photon parameter space as a function of the signal acceptance. We require that 10 signal events pass all selections with the indicated acceptance and the full Acceptance 10 −1 10 −3 10 −5 10 −7 Figure 2: Sensitivity to the dark photon mass and kinetic mixing, ε, under the assumption that selections with the indicated signal acceptance are background-free and can be discovered with 10 signal events after cuts with L = 50 ab −1 at Belle II. We also show model-independent limits on the dark photon parameters [32,88], as well as the band favoured by the (g − 2) µ anomaly [88].
integrated luminosity of Belle II, L = 50 ab −1 . It is evident that even if searches require very stringent selections to eliminate backgrounds, leading to a very low signal acceptance (∼ 10 −7 ), they can still be sensitive to parameters motivated by (g − 2) µ .
The coupling of A to the LLP current, g D A µ J µ χ , suggests that χ is charged under a dark U (1) gauge symmetry. This gauge symmetry must be broken to allow χ to decay as an LLP back to SM particles. This can be accomplished in several ways: for example, the χO SM /Λ n−4 operators allowing χ to decay could have an additional insertion of some dark Higgs vacuum expectation value, h /Λ. Alternatively, χ could mix with a heavy Majorana singlet state ξ after dark symmetry breaking, and it is secretly a ξO SM /Λ n−4 effective operator that allows χ to decay to the SM. We are agnostic about which mechanism allows χ to decay to the SM in the dark photon model, in part because we treat the χ lifetime as a free parameter and these model-dependent considerations simply lead to adjustments in the relationship between the χ lifetime and the EFT scale, Λ.
We finally note that the interactions of Eq. (7) imply an irreducible branching fraction of A into SM final states f , which scales as BR(A →ff ) ∼ ε 2 e 2 /g 2 D . For sufficiently large values of ε, this suggests that A can be probed with the usual visible dark photon searches such as Ref. [15]. Since the sensitivity of visible searches in this limit scales as ε 4 , the visible decays are negligible when ε 0.01 and g D ∼ 1 even with the expected full Belle II luminosity. We therefore focus on the parameter space where the dominant sensitivity is instead provided by the multi-track searches described in Sec. 5.

LLP Production: B Meson Decays (Dark Higgs)
We consider a Higgs-portal coupling between a scalar mediator, φ, and the LLP, χ. After electroweak symmetry breaking, the SM and dark Higgs states mix with angle θ, giving a tree-level coupling to quarks and LLPs: where the coupling y of φ to LLPs is assumed to be sufficiently large such that BR(φ → χχ) ≈ 100%. In loop-induced B decays, the dark Higgs benefits from the large coupling of φ to top quarks, and it therefore has a large production rate in flavour-changing decays of b quarks, b → sφ. For m φ m b , the inclusive φ production rate in B meson decays is [89] BR The cross sections for B + B − and B 0B0 production are each 550 pb at the Υ(4S) resonance. Before selections, and with 50 ab −1 of integrated luminosity, this gives a number of χχ pairs equal to For larger φ masses, it is more appropriate to sum over exclusive decay modes B → Kφ, where we include various kinematically accessible kaon resonances. For 4.3 GeV < m φ < 4.78 GeV, only a single decay into the lightest kaon occurs with a rate of about 10% of the inclusive rate from Eq. (11), while for somewhat smaller φ masses the sum of exclusive production modes is comparable to the rate predicted by the quark-level calculation. In our analysis, we use the rates from Ref. [89].
We have used FeynRules to implement a UFO model containing Υ(4S), B ± , and various K ± mesons, along with the dark Higgs φ, and we simulate Υ(4S) → B + B − , B + → K + φ, while the other B − decay mode is unspecified. 5 The effective flavour-changing φ − B − K coupling is implemented at tree level, and the relative couplings of the different kaon states are chosen to match the relative exclusive decay rates from Ref. [89]. Since we are interested in LLPs that are sufficiently heavy that they decay into multiple charged particles, we are mostly interested in the case where m φ is not much smaller than m B − m K and hence the production of φ is well approximated by summing over exclusive decay modes.

Simulation and Analysis
In this section, we demonstrate that searches for displaced vertices containing multiple charged tracks are broadly sensitive to the generalized LLP decays described by EFT operators as described in the previous sections. We focus on two production channels, radiative return and rare B meson decays, which we simulate as described in Sections 4.1 and 4.2. The resulting LLPs are decayed in MadGraph and hadronized in Pythia [90], or manually by matching to chiral perturbation theory as discussed in Appendix B. We analyze these events in the context of the Belle II experiment at the SuperKEKB facility.
The Belle II experiment studies e + e − collisions at √ s = 10.58 GeV (with E e + = 4 GeV and E e − = 7 GeV) using a detector described in detail in, e.g., Refs. [91,92]. Over its lifetime, Belle II will collect 50 ab −1 of integrated luminosity, yielding unprecedented sensitivity to rare processes. For our LLP searches, the most important detector subsystems are the vertex detector (VXD), with layers starting only 1.4 cm from the interaction point and going out to 17 cm, and the central drift chamber (CDC), which extends radially out to 1.1 m. Both the VXD and CDC provide polar angular coverage between 17 • and 150 • in the lab frame. The displaced vertex (DV) detection efficiency is transverse radius-dependent and we follow Ref. [29] in taking 6 for −55 cm < z < 140 cm and zero outside. We also account for the 42 mrad beam "halfcrossing" angle [91] as described in Ref. [29]. In identifying tracks, we demand charged particles to be within the geometric acceptance and to have transverse momentum p T > 100 MeV; we also assign a detection efficiency per track of 0.9. While track-finding efficiencies are strongly dependent on p T , impact parameter, and decay position, the 0.9 value we use is a simple way of incorporating tracking inefficiencies that gives, on average, results consistent with expectations from the Belle II detector [93]. Same-sign tracks with angular separation θ < 0.05 are merged and counted as one, although we find this has little effect on our results.
We use these definitions to find tracks and displaced vertices in our event samples and study the acceptance (and the resulting reach) for different choices of the minimum number of DVs per event, N dv , and the minimum number of tracks per DV, N tr .
We estimate the sensitivity of LLP searches assuming there are no backgrounds, and thus a 10-event signal detection is sufficient for discovery. We expect that this is a good approximation for searches requiring N dv ≥ 2 and N tr ≥ 3. While even this stringent selection is sensitive to many LLP production and decay channels, it is also interesting to consider looser selections, which may require additional background mitigation strategies. We will describe the relevant backgrounds and possible additional cuts in Sec. 6, but it is important to note that these cuts can be used to reduce backgrounds without significantly impacting signal acceptance across a broad range of LLP parameter space.

Radiative Return Production of a Dark Photon
We first discuss the sensitivity of multi-track displaced vertex searches to e + e − → γA with A → χχ. The signal acceptance and cross-section sensitivity are shown for the following LLP decays, expressed in terms of parton-level operators: • Fully leptonic decay of scalar LLP χ → µ + µ − e + e − in Figs. 3 and 4; • Semileptonic, fully visible decay of fermion LLP χ → edu in Figs. 5 and 6; • Semileptonic, partially invisible decay of scalar LLP χ →νedu in Figs. 7 and 8; • Fully hadronic decay of scalar LLP χ →ūuūu in Figs. 9 and 10.
Fully leptonic LLP decays. The most striking signals arise in fully leptonic decays. We consider for concreteness the decay of a scalar LLP, χ → µ + µ − e + e − ; other flavour combinations are equally motivated. In this case, each event contains two vertices with four tracks and the acceptance is mainly determined by the LLP decay length. In the left panel of Fig. 3 we show the acceptance as a function of LLP mass for cτ χ = 1.2 cm. The probability of observing two vertices is roughly the square of a single vertex, which is reflected in the relationship between the N dv = 1 and N dv = 2 acceptances. Finally, note that requiring additional tracks has a negligible impact on the acceptance, since each decay always produces four tracks. A similar pattern is evident in the right panel of Fig. 3, which shows the acceptance as a function of cτ χ for m χ = 2 GeV.
In Fig. 4 we estimate the sensitivity of the multi-track search for m A = 4 GeV and 8 GeV in the left and right columns, respectively; the upper (lower) row shows analyses with N dv ≥ 2, N tr ≥ 2 (N dv ≥ 2, N tr ≥ 3). Each panel shows the 10 signal event cross-section reach in the m χ − cτ χ plane, with darker colours corresponding to smaller cross sections. The dashed white lines correspond to σ = 10 −2 fb, which translates to a kinetic mixing reach of ε ∼ 10 −5 ; this can be compared with ε 10 −2 needed to explain the (g − 2) µ anomaly with m A GeV. We see that our proposed searches offer excellent sensitivity independently of N tr and N dv selections. Moreover, the reach extends even to very low LLP masses 1 GeV.
Semileptonic LLP decays. The first semileptonic decay mode we consider is of a fermion LLP χ → e −d u via the (ū R γ µ d R ) ¯ R γ µ χ operator. Like the fully leptonic decay above, this decay mode is fully reconstructible, although unlike the fully leptonic case the dominant decay modes for small LLP mass are χ → e − ρ + and χ → e − π + , which feature only two charged tracks. However, at higher LLP mass the hadronization ofdu can give rise to larger track multiplicities, especially due to intermediate vector resonances such as the a + 1 . We employ a τ -inspired exclusive meson decay model outlined in Appendix B.1 for m χ ≤ 2 GeV, and Pythia for showering and hadronization for larger masses. 7 We show the acceptances for this model in Fig. 5. The acceptances for the N dv ≥ 2, N tr ≥ 3 search are reduced by over an order of magnitude relative to the leptonic case due to the reliance on hadronization to produce extra tracks. The acceptance for this search is greatly reduced at low m χ because the dominant χ decays have two tracks in this mass range. As before, the acceptance of N dv ≥ 2 events is roughly the square of the N dv ≥ 1 case.
Although the acceptances are smaller than for fully leptonic decays, they are still not neg-    ligible. The cross-section sensitivities for several dark photon mass benchmarks are shown in Fig. 6. While the N dv ≥ 2, N tr ≥ 2 sensitivity shown in the upper panels is qualitatively similar to the fully leptonic case, the N tr ≥ 3 search is limited at low LLP masses as explained above. The signal acceptances for lower m A masses and N tr ≥ 3 (lower left panel) are more negatively affected since m χ is restricted to smaller values such that the A →χχ is still kinematically allowed. Despite this, we see that the multi-track search still offers sensitivity to tiny cross sections across a broad range of parameter space.
Next, we consider an even more challenging LLP semileptonic final state: χ →νedu. For this LLP decay mode, an inclusive search strategy is required because the decay is not fully reconstructible. The signature includes missing energy in addition to the lower track multiplicity fromdu and, once again, these challenges are clearly reflected in the acceptances in Fig. 7. The resulting search reach is shown in Fig. 8. However, the results are qualitatively similar between the two semileptonic χ decay modes. Note that there are larger theoretical uncertainties on the hadronization of this scalar operator compared to the vectorū R γ µ d R operator above. In Figs. 7 and 8 we use Pythia for showering and hadronization, and we provide an alternative model based on chiral perturbation theory (as well as comparisons with Pythia) in Appendix B.2.
Fully hadronic LLP decays. For completeness, we also consider fully hadronic LLP decays. As a benchmark, we study the final states with χ →ūuūu. Hadronization of this final state produces a pair of neutral pions most of the time, so the requirement of multiple tracks per vertex can penalize the acceptance. Branching fractions into charged states is also very sensitive to the hadronization model, so our predictions for these decays are subject to a large systematic uncertainty. We emphasize that this does not affect the experiment's ability to carry out the proposed searches. In Fig. 9 we show the acceptance as function of LLP mass in the left panel and as a function of its lifetime in the right panel for event samples generated using Pythia to hadronize the final states. As we show in Appendix B, this most likely provides an optimistic    estimate of the acceptance, particularly at low LLP masses 2 GeV. Fig. 10 shows the crosssection sensitivity for the same hadronization procedure. As in the previous examples, multitrack searches can be sensitive to cross sections far smaller than those motivated by solutions to (g − 2) µ . Note that we do not show the m A = 4 GeV case because of large hadronization uncertainties for m χ ≤ m A /2 = 2 GeV.

Dark Higgs Production in B Meson Decays
We now turn to LLP production in the decays of a dark Higgs. The full process is B → Kφ, φ → χχ, where we sum over various kaon resonances. The primary effect of the change in production mode is that it alters the boost of the LLPs, which shifts the values of cτ χ corresponding to the optimal decay length sensitivity of 1 mm-50 cm for a tracker-based displaced vertex analysis. However, the acceptances are otherwise uncharged, largely because the track p T selections are easily passed independent of the LLP boost for the χ masses we consider.
Because the sensitivities of our proposed displaced vertex searches are largely unchanged from the dark photon production mode, we only show results for one of the LLP decay modes, namely the fermion LLP χ → edu decay originating from the (ū R γ µ d R ) ¯ R γ µ χ operator. We have selected this decay mode for illustrative purposes and because of the robustness of our modelling of the hadronic part of the decay over a wide mass range; however, we emphasize that all possible χ decay modes allowed by the EFT are comparably motivated.
In Fig. 11, we show the acceptances for various requirements on N dv and N tr for the B decay production mode. Most striking is the fact that the right pane of Fig. 11 is nearly identical to the right pane of Fig. 5, which shows the acceptance dark photon production and the same decay mode. The value of cτ χ for which the acceptance peaks is shifted between the two production modes, reflecting the different characteristic boosts for the LLPs, but the curves are otherwise very similar. We show in Fig. 12 the sensitivity to the branching fraction of B ± → K ± φ assuming a 100% branching fraction of φ →χχ, finding once again similar sensitivity to the dark photon production mode in terms of cross section. We again show the 10-event sensitivity     We fix m A = 8 GeV. The white dashed contours correspond to a production cross section of 10 −2 fb, which is roughly equivalent to ε ∼ 10 −5 for m A ≤ 8 GeV.   to a cross section of 10 −2 fb, which in the dark Higgs model corresponds to a mixing with the SM Higgs of θ ∼ 10 −4 for m φ = 4.5 GeV. This analysis demonstrates that sufficiently inclusive searches can have comparable acceptances regardless of LLP production mode.

Backgrounds
Backgrounds for LLP searches are notoriously difficult to estimate [39], particularly in the absence of access to a full detector simulation. Data-driven methods are also needed to model extremely rare processes that may not be captured in simulations. In this section, we do not attempt to accurately model backgrounds for LLP searches at Belle II, but instead we enumerate the main expected sources of background, and discuss how they can be characterized and mitigated in multi-track displaced vertex searches. As a general rule, the backgrounds to LLP searches are very rare, and so increasing either the multiplicity of vertices or the number of tracks per vertex can suppress them.
The backgrounds fall into three main categories. First, there are backgrounds from fake vertices and tracks, whether arising from collisions of primary particles with gas or material in the detector, accidental track crossings, or spurious tracks that arise from mis-reconstruction of tracker hits. Second, there are backgrounds from heavy-flavour decays such as b and c hadrons, which mimic the signal by having relatively large masses and high track multiplicities. Finally, there are backgrounds from hadrons with strange quarks, such as kaons and Λ baryons. While these decays typically give low track multiplicities and can often be fully reconstructed and rejected, they are produced in sufficient numbers that rare decays or mis-reconstruction of the vertices could populate the signal region. All masses, lifetimes, and branching fractions of SM particles are taken from the Particle Data Group world averages [94]. GeV. The white dashed contours correspond to a production cross section of 10 −2 fb, which is roughly equivalent to dark Higgs mixing θ ∼ 10 −4 for m φ ≤ 4.5 GeV.

Combinatoric and Detector Backgrounds
Displaced vertices can originate from particle interactions with the detector material, such as photon conversions into e + e − pairs, as well as random overlaps of high-impact-parameter tracks from the decays of kaons, pions, and strange baryons. Spurious vertices from material interactions can be removed by vetoing vertices that originate in material-dense regions of the detector. Additionally, photon conversions and material interactions should typically give rise to low multiplicities and small invariant masses of the tracks originating from the vertex, and so increasing the N tr requirement and applying a modest cut on the vertex track invariant mass can remove these vertices.
More problematic are random crossings or mis-reconstructed tracks that are accidentally combined into a displaced vertex. In the BABAR search for N tr = 2 vertices [49], this constituted the dominant background. Based on the total BABAR trigger rate and the number of background events, we calculate the probability of a fake single, two-track vertex to be approximately 10 −6 for that analysis, which would correspond to 10 5 such vertices with the full Belle II integrated luminosity. In reality, the fake-vertex background will be larger because of the larger beam backgrounds at Belle II. Nevertheless, random-crossing vertices become rarer for larger vertex multiplicities and track multiplicities per vertex, and so requiring N dv ≥ 2 and/or N tr ≥ 3 should render this background negligible. Fortunately, the fake vertex rate can be estimated from data: for example, ATLAS and CMS pair high-impact-parameter tracks and vertices from different events to efficiently estimate the magnitude of this background in LHC searches [95,96]. A similar approach at Belle II would likely be effective.

Heavy-Flavour Backgrounds
Heavy-flavour backgrounds, such as B mesons, D mesons, and τ leptons, are pernicious because they give rise to genuine high-mass, multi-track displaced vertices. Their typical decays cannot be fully reconstructed, which means they cannot be readily vetoed. Fortunately, the lab-frame decay lengths of these particles at Belle II are sub-mm: βγcτ is approximately 0.14 mm for B mesons (βγ ∼ 0.3) and 0.26 mm for τ leptons (βγ ∼ 3). D mesons have different kinematics depending on whether they originate from continuum cc production or B meson decay, but we expect a typical βγ ∼ 1 (with a maximum value of 3), which corresponds to a lab-frame decay length of 0.3 mm (0.9 mm).
With the illustrative minimum vertex displacement we used in Sec. 5 of 2 mm, we expect to reduce B backgrounds by 2 × 10 −7 , τ backgrounds by 4.5 × 10 −4 , and charm backgrounds by approximately 10 −3 (although this depends on how many charm mesons populate the highestmomentum tail). Increasing the minimum vertex displacement to 4 mm could reduce even charm backgrounds to the 10 −6 level, at which point requiring N dv ≥ 2 would effectively suppress heavy-flavour backgrounds. Tightening the minimum vertex displacement would eliminate some sensitivity to signals with the shortest values of cτ , but would otherwise have a negligible impact on signal acceptance. Because heavy-flavour backgrounds follow an exponential distribution in decay length and have kinematics that can be well-measured in other processes, they can be readily estimated by looking at sidebands that are just below the signal region in displacement.
While raising the minimum required vertex displacement can effectively suppress heavyflavour backgrounds, it may be possible to mitigate them using other means depending on the signal being searched for. As a result, sensitivity to smaller signal lifetimes could be possible at the expense of making the searches less inclusive. Examples of possible strategies could include requiring multiple leptons per displaced vertex, or (in the case of charm backgrounds) vetoing vertices with kaons. The leading multi-charged-pion decay modes of D mesons (such as D ± → π + π − π − ) occur at the 10 −2 level, compared to ≥ 20% for 3-track decays including kaons. In τ events, there should be a very low multiplicity of particles that are unassociated with the displaced vertex, since most taus decay to a single track plus missing momentum. B backgrounds can be suppressed using event-shape observables like thrust [97] and Fox-Wolfram moments [98], which can identify the characteristic kinematics and event topologies associated with B mesons; however, if this strategy were pursued, it would be important to test whether this would also limit sensitivity to LLPs produced in B decays, such as B → Kφ, φ → χχ.

Strange-Quark Backgrounds
Finally, we turn to strange-quark backgrounds such as K + , K 0 S , K 0 L , and Λ. These typically feature low-track-multiplicity decays, but have sufficiently long lifetimes that they frequently give displaced decays inside of the signal region.
Charged kaons. With cτ ≈ 0.4 m, most K ± are stable through the tracking system. Kaons are typically produced with an O(1) boost, which means that we expect 1%-10% of charged kaons to decay in the signal region. Given that we expect about one kaon per hadronic event [99], this translates to ∼ 10 11 kaons with the full Belle II integrated luminosity. The most common, dangerous decay mode is K + → π + π + π − with branching fraction of 5.5%. This is, however, a fully reconstructible decay and can be vetoed: if the veto has an inefficiency of 10 −3 for K ± , this should be sufficient to render the background negligible for N dv ≥ 2. Furthermore, tightening the track requirement to N tr ≥ 4 would eliminate these backgrounds while having little effect on neutral signal DVs. Other problematic decay modes include K + → π + π − e + ν e with a branching fraction of 4 × 10 −5 , K + → π + π + π − γ with a branching fraction of 7 × 10 −6 , and K + → π + π 0 e + e − with a branching fraction of 4 × 10 −6 . Other multi-lepton decay modes occur at the level of 10 −8 . These are sufficiently small that the background should be negligible with the requirement of N dv ≥ 2. If some kaon backgrounds persist in the signal region after these selections, they can be efficiently removed with a cut on the vertex track invariant mass of 0.5 GeV. In Fig. 13, we show the impact on the acceptance of such a track-mass selection when the LLP decays semi-visibly, finding that the signal acceptance is largely unchanged except at small LLP masses where the typical visible invariant mass is below the value of the cut. The effect of a track-mass cut is even less pronounced for fully visible LLP decays, since the invariant mass of the charged particles in the final state is larger.
Neutral kaons. The CP -even K 0 S state has cτ ≈ 3 cm, which is right in the middle of the signal region. K 0 S is therefore potentially the most dangerous background. It typically decays to two charged particles, and so a requirement of N tr ≥ 3 greatly suppresses this background. The dominant multi-track decay mode is K 0 S → π + π − e + e − with a branching fraction of 5 × 10 −5 . This decay is, however, fully reconstructible, and even a veto inefficiency of 10% for this decay would still render the background negligible with a requirement of N dv ≥ 2. An even more efficient veto, or a veto of vertices with this particular particle combination in the final state, could even make this background manageable with N dv ≥ 1 as the requirement. However, this depends on the detailed rate of track mis-reconstruction in K 0 S decays.
The CP -odd K 0 L state has cτ ≈ 1.5 m. Only about 1% of K 0 L decays will occur in the signal region for O(1) boosts. The most dangerous decay mode is K 0 L → π + e − e + e − ν e with branching fraction 10 −5 . Other high-track-multiplicity modes, such as K 0 L → π + π − e + e − , have branching fractions < 10 −6 and are fully reconstructible, and so a veto on the K 0 L peak should efficiently N dv ≥ 2, N tr ≥ 3 no mass cut w/ mass cut χ →νedu, m A = 8.0 GeV, cτ χ = 1.2 cm Figure 13: Comparison of signal acceptance with and without a vertex track-mass cut > 0.5 GeV to reduce kaon backgrounds. reject them. Together with the lower acceptance of K 0 L decays, the requirement of N dv ≥ 2 should suppress the K 0 L backgrounds to negligible level. As with K + decays, however, if some backgrounds remain then a cut on the vertex track invariant mass of 0.5 GeV should eliminate the remaining kaons in the signal region. As we discussed above and show in Fig. 13, this cut should not significantly affect the acceptance for most signal benchmarks.
Strange baryons. Baryon production rates are subdominant to the meson production rates. For example, the off-resonance inclusive Λ cross section is about 0.3 pb at √ s = 10.52 GeV [100], which is smaller than the total continuum quark cross section of 2.3 pb by almost an order of magnitude. Similarly, the inclusive B meson branching fraction to Λ c baryons (of which 40% decay to Λ baryons) is at the percent level. However, with cτ Λ ≈ 8 cm, Λ baryons still constitute a major background. Since Λ always decay to at most two tracks, requiring N tr ≥ 3 should greatly suppress this background. Additionally, the dominant 2-track decay Λ → pπ − is fully reconstructible, and these vertices can be efficiently rejected. The dominant non-fully-reconstructible 2-track decay is Λ → pe −ν e with a branching fraction of 8 × 10 −4 . A more targeted search with particle identification selections on the vertex tracks could render a N dv ≥ 2, N tr ≥ 2 track search background-free with respect to Λ decays.

Discussion and Conclusions
We have proposed several strategies for inclusive, long-lived particle searches at the Belle II experiment that should greatly expand sensitivity to GeV-scale hidden sectors. To motivate such new searches and characterize their reach, we have developed an EFT framework to classify LLP decay modes in a gauge-invariant manner, finding that multi-track displaced vertex analyses can have sensitivity to a wide range of models. We have also considered multiple well-motivated LLP production modes, including dark photon and dark Higgs portals, finding good sensitivity to couplings of interest from dark matter and muon g − 2 models. Crucially, many of these extended hidden-sector models would not have been discovered in existing searches, which tend to target specific exclusive production and decay modes.
While we have focused on the Belle II experiment, whose substantial integrated luminosity and clean environment should allow for the low-background reconstruction of multi-track displaced vertices, the hidden-sector models we presented here should also lead to large production rates at LHCb. Because of the large boosts characteristic of forward experiments at hadron colliders, LLP lifetimes must be shorter than at Belle II to give displaced decays in the tracking system. This gives a relatively large increase in B backgrounds, making inclusive searches somewhat more challenging; indeed, most inclusive LLP searches at LHCb only have sensitivity for LLP masses above approximately 20 GeV [101,102] because the large vertex mass can be used to efficiently reject B backgrounds, although the reconstruction of displaced resonances is possible at lower masses [51]. While it is possible for LHCb to improve on low-mass exclusive and inclusive signatures [103,104] and this merits further study, we nevertheless expect Belle II to have a unique role in covering intermediate-lifetime hidden sectors in an inclusive manner.

Acknowledgments
We are grateful to Victoria Lloyd for collaboration at an early stage of this project. We thank Torben Ferber and Abi Soffer for helpful discussions of LLPs at Belle II, and Olivier Mattelaer and Qiang Li for useful correspondence regarding MadGraph and the ISR plugin. This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. This work is supported by the U.S. National Science Foundation under Grant PHY-1820770.

A Effective Field Theory Framework of LLP Decays
In this appendix we describe in more detail how we classify LLP decay operators in an electroweakgauge-invariant EFT. This task is made simple by the fact that dark sector LLPs are neutral under SM gauge symmetries, which allows us to repurpose established results from SMEFT (the collection of higher-dimensional operators constructed from SM fields alone) [61][62][63][64] and SM + sterile neutrino EFT (νSMEFT) [65][66][67]. The references provided here contain comprehensive lists of operators which we do not reproduce; instead we describe their general structure at different operator dimensions and provide explicit examples that give rise to the specific decays studied in the main part of the paper.

A.1 Spin-0 LLP
If the LLP χ is a scalar, then the decay operator has the generic form χO SM where O SM is also a Lorentz scalar and therefore belongs to SMEFT. From Eq. (3) we see that we want to focus on operators with dim O SM ≤ 7.
At dimension 4 (corresponding to dim χO SM = 5), we can multiply any term in the renormalizable SM Lagrangian by χ: the total derivative terms G a µν G a µν , B µν B µν and W a µν W a µν , and pseudo-scalar couplings to fermions. These lead to φ →f i f i (for all kinematically accessible fermions other than ν) and φ → γγ, gg, with other decays (e.g. to EW gauge bosons) being kinematically forbidden for LLPs of interest to B-factories; these final states are shown in the first row of Tab. 1.
At dimension 5 (corresponding to dim χO SM = 6) there is only a single operator χ(L i H)(L j H), which leads to φ → ν i ν j at low energies after EW symmetry-breaking. This decay is shown in the second line of Tab. 1.
At dimension 6 (corresponding to dim χO SM = 7), there are 59 baryon number B-preserving and 4 B-violating operators [61,64]. These interactions are listed in Tab. 1 of Ref. [64] and Eq. (1) of Ref. [62], respectively. The LLP decay final states enabled by these operators are listed in the third line of Tab. 1. The leptonic, semi-leptonic and hadronic decay operators for scalar LLPs used in the body of the paper are generated by the following operators at this dimension: where lower (upper) case fields are right-(left-) handed Weyl fermions, and the Roman indices are flavour labels.
All dimension-7 operators (corresponding to dim χO SM = 8) violate B − L, and some break B or L individually [63,105], making them particularly interesting. However, we leave these for a future study as symmetry-violating operators deserve a careful treatment as they are likely to be strongly constrained by other observables.
The dim χO SM = 6 operators are listed in Tab. 1 of Ref. [67]. We use (d p γ µ u r )(χγ µ s ) (15) to generate semileptonic decays of a spin-1/2 LLP as described in the main text. We show this and other possible final states produced at this dimension in the first line of Tab. 2.
Many operators at dim χO SM = 7 (collected in Eqs. (14)- (17), (22) and (25) of Ref. [67]) are obtained by adding covariant derivative to a dimension-6-like operator. Thus these operators generate the same final states as dimension six, possibly with an addition gauge boson (with γ and g being the most relevant in our mass range). The corresponding decay channels are listed in the second line of Tab. 2.

B Hadronic Final States
In considering LLP decays to hadronic final states, we need to translate parton-level operators for LLP decays into operators involving mesons and baryons. For large LLP masses, it is appropriate to decay the LLP to quarks and then rely on a parton shower/hadronization model to obtain the hadronic final state. Conversely, for small LLP masses, chiral perturbation theory (ChiPT) and the Hidden Local Symmetry (HLS) frameworks provide a way of estimating the mapping of quark operators to meson operators. The boundary between these two treatments is not clearly defined, although we expect that the parton model should give reasonable results for q 2 GeV 2 , where q is the 4-momentum flowing into the composite quark operator. For each operator involving quarks, we perform calculations with both chiral perturbation theory and a parton shower + hadronization, and compare the results in the region of intermediate hadronic momentum where neither treatment is fully valid. We find that in all cases we do get charged track multiplicities of three or more at appreciable rates, although in certain regimes the two treatments differ in the degree to which this is true.
We emphasize that, for the purposes of the experimental search, the theoretical uncertainty on the hadronization model is irrelevant: multi-track searches are well motivated by both approaches to hadronization, and as long as the search results are transparent in terms of the precise model used in setting limits (or with efficiencies provided for displaced object reconstruction), the results can be reinterpreted and applied to other, more refined hadronization models, particularly as theories of hadronic GeV-scale decays continue to be refined.

B.1ū R γ µ d R Operator
This operator closely resembles the hadronic current in leptonic meson decays. For example, for exclusive decays into hadrons we can use the hadronic matrix elements where µ is the ρ meson polarization vector, f π = 0.130 GeV, and g ρ = 0.119 GeV 2 [110]. For the operator 1 Λ 2 (ū R γ µ d R ) ( R γ µ χ), we can then calculate the decay rates for χ → e − π + and χ → e − ρ + , ρ + → π + π 0 (neglecting the pion mass): where V ij is the Cabibbo-Kobayashi-Maskawa (CKM) matrix. These results are essentially the same as the exclusive semileptonic decay rates of heavy neutral leptons [78]. For m χ m ρ , we find that the ratio of decays to ρ + relative to π + approaches two. Thus, while both of these decay modes produce a single charged pion, the majority of LLP decays actually give an additional π 0 , and so an exclusive reconstruction of the LLP mass will cut out most of the signal events. In our analysis, we neglect the decay into K + , which is both Cabibbo-and phase-space-suppressed.
Higher charged pion multiplicities can occur once we take into account heavier vector mesons. It is well established that the 3-pion τ decays predominantly arise from a 1 intermediate states [111]. Similarly, the decay χ → e − a + decay χ → e − π + π − π + . The calculation of this process is complicated by the large, energydependent a + 1 width, as well as nearby broad resonances such as the a 1 . There are several models of 3-pion τ decays. We choose the model of Ref. [112] due to its simplicity and the readiness with which it can be implemented in the UFO format. There are some aspects of the model that cannot be simply accounted for in MadGraph (for example, the renormalization group running of the a 1 mass and the energy-dependent a 1 width); however, we have found that our somewhat crude implementation of the model predicts the τ − → π − π + π − ν τ decay width in agreement with the observed value within 15%, and that the kinematic acceptance for our proposed displaced vertex search in Sec. 5 is insensitive to model details. Away from the τ mass, the variation in the 3-pion decay rate is of O(1) when we change model parameters.
We implement an exclusive hadronic decay model in the UFO format and simulate events using MadGraph for m χ < 2 GeV. At higher masses, we directly implement the decay χ → e − ud in MadGraph, and perform showering and hadronization with Pythia. In the region 1.5 GeV < m χ < 2 GeV, we use both methods and compare their acceptance and efficiency in our analysis.

Meson Model Details & Validation
Our model consists of the LLP, χ, as well as the ρ ± , ρ 0 , π ± , π 0 , a ± 1 , and a ± 1 mesons, and is based on Ref. [112]. The charged pions and vector mesons are coupled to χ by replacing the hadronic currents with the meson fields to give the matrix elements in Eqs. (16)- (17). The a 1 hadronic matrix element is the same as for the ρ, but using g a = 0.250 GeV 2 as the decay constant.
The a 1 interaction Lagrangian is taken to be where sin θ ≈ 0.463 and being the field corresponding to a 1 . The coupling g a 1 ρπ is fixed by obtaining the observed a 1 width, and we obtain g a 1 ρπ = 54 GeV −1 .
The Lagrangian describing the ρ coupling to pions is where g ρππ ≈ 6 gives the correct ρ width [113]. The authors of Ref. [112] found that the inclusion of the a 1 , or a 1 (1640), is important to get good agreement with 3-pion τ decays. The hadronic matrix element is the same for the a 1 as for the a 1 , but with an additional multiplicative complex constant α = −0.31 ± 0.32i. The a 1 is taken to have the same hadronic couplings as the a 1 . In our model, we use m a 1 = 1.232 GeV, Γ a 1 = 0.431 GeV, m a 1 = 1.655 GeV, and Γ a 1 = 0.254 GeV.
The model of Ref. [112] implements a form-factor suppression of fully hadronic vertices; however, for the kinematic range of interest to us, this is not a significant effect and would greatly add to the complexity of the model and so we neglect it. The model of Ref. [112] also includes the σ, or f (500), scalar meson. Its contribution to the 3-pion decays of the a 1 meson are only at the 20% level, and as the σ meson is extremely broad (its width is comparable to  Figure 14: Ratio of acceptances of a N dv ≥ 1, N tr ≥ 3 search with a m A = 8 GeV, cτ χ = 1 cm signal for the original χ → π + π − π + e − decay model compared to the alternate model where we have adjusted some of the intermediate resonance masses. This ratio is within 5% of unity over the mass range for which we use the exclusive meson decay model. its mass), it is not straightforward to model. Therefore, we do not include the σ meson in our model. We implemented this model in the UFO format using FeynRules, and we use MadGraph to simulate events. To validate the model, we simulated the process χ → π + π − π + e − for m χ = m τ and compared the results to CLEO data of τ − → π − π + π − ν τ [112,114]. Our model predicts a value for the 3-pion partial width that is 87% of the PDG 3-pion width of the τ ; this is remarkably good agreement, and the discrepancy is at the level of the σ meson contribution that we neglect in our model. We can then calculate the χ → π + π − π + e − partial width at other values of the χ mass to determine the branching fraction to 4 charged tracks, applying the same correction factor of 1/0.87 to account for decay modes not included in our model.
We also used our MadGraph model to compute the m πππ invariant mass spectrum from χ → π + π − π + e − and compare with CLEO data from τ − → π − π + π − ν τ decays. Our spectrum is slightly broader than theirs and with the peak shifted 0.2 GeV to the right. We therefore found an alternative model benchmark point where we obtained good agreement with the data if we shift the masses to m a 1 = 1 GeV and m a 1 = 1.35 GeV. Clearly this is unphysical, but as a phenomenological model it allowed us to obtain a 3-pion spectrum that is consistent with the data from τ decays and we can use this to estimate a "systematic" on the signal acceptance by varying these parameters. We repeated our calculation of the kinematic and geometric acceptance of a search for a single displaced vertex with at least three tracks on the sample with 3-pion decays; we show the ratio of acceptances for the two models in Fig. 14. There is no appreciable difference in the acceptance of the two models, which gives us confidence in the robustness of our results in spite of the ad hoc nature of some aspects of our hadronic model.
We can also use our alternative model benchmark to calculate the 3-pion partial width of χ, normalizing to the τ -mass value as before. We find that the two models agree for m χ 1.5 GeV, but unsurprisingly the alternative model gives a larger 3-pion width at lower mass due to the lower a 1 /a 1 masses, which allows for on-shell contributions for a wider range of χ masses. The discrepancy between the two models can be as much as a factor of 3, although we expect the original model to give a better estimate of the 3-pion partial width. This does suggest that the uncertainty on the 3-pion width from our method is of O(1), but given that we see extremely good sensitivity of our multi-track vertex search, this does not alter the motivation  for the search or our general conclusions.
Finally, we compute the acceptances of our proposed displaced vertex searches from Sec. 5 twice: using the meson decay model, as well as with decays of χ to quarks using Pythia for showering and hadronization. We expect that Pythia should give the more correct description at larger LLP masses, while exclusive meson decays provide a better description for small LLP masses. We compare the acceptances of the two approaches for m χ = 2 GeV in Fig. 15, finding that the two descriptions agree very well for the N dv = 2, N tr = 2 analysis, and agree within a factor of 2 for the N dv = 1, N tr = 3 analysis (the acceptance of the N dv = 2, N tr = 3 analysis is approximately the square of the N dv = 1, N tr = 3 analysis). Given the challenge of correctly modelling hadronization in this intermediate mass range, this gives us confidence in the robustness of our results. We find a similar comparison between the two hadronization methods throughout the mass range 1.5-2 GeV.

B.2d L u R Operator
The mapping of quark-level operators into ChiPT is discussed in, e.g., Ref. [115] in the context of GUT-scale B-violating interactions leading to proton decay. First, we need to identify the chiral transformation properties of the operator whose EW-symmetric form is given in Eq. (14) and the quarks q L(R) transform as a 3 of SU (3) L(R) ; we also introduced a 3 × 3 flavour-space matrix C lequ . The choice generates χ(ν L e R )(d L u R ) used in the main text. We treat C lequ as a spurion of the chiral SU (3) L × SU (3) R symmetry, which transforms as where L (R) are elements of SU (3) L (R) . This transformation rule makes it easy to identify the leading operator corresponding to O lequ in chiral perturbation theory, since C lequ transforms in the same way as the quark mass matrix. Therefore O lequ maps onto where Σ = exp(2iΠ/f π ) and Π is the meson field 8 The exponentiated field transforms as The Wilson coefficient B can be estimated by analogy with the chiral mass term which has the same structure as the quark part of O lequ and gives rise to meson masses; this comparison yields B ∼ m 2 π f 2 π /(2(m u + m d )) ∼ (300 MeV) 3 .
In the last step of Eq. (25) we only showed the leading term in 1/f π ; higher-order terms contain more mesons, and at O(1/f 3 π ) terms with multiple charged mesons appear.
We are interested in the decays of LLPs with mass of O(GeV) which is close to the cut-off of ChiPT of ∼ 4πf π , implying the possible importance of vector mesons in multi-track decays. We construct operators containing vector mesons next.

Including Vector Mesons
Multiple charged tracks may be generated through the production and decay of vector mesons. We implement vector meson couplings using the Hidden Local Symmetry (HLS) formalism [116], in which the building blocks are ξ L,R : These fields are related to usual ChiPT field Σ via Unlike Σ, however, they transform non-trivially under the entire (global + hidden) symmetry group: where h(x) ∈ SU (3) V is an element of the hidden gauge group, and L (R) ∈ SU (3) L (R) as before. In order to find an operator that can lead to vector meson production, we need to use a modified spurion analysis which ensures that the result is invariant under the hidden gauge symmetry, and it is not degenerate with using Σ (i.e., does not just reduce to Eq. (25)). One way to achieve this is to start from the (hidden) gauge-covariant derivatives where g V ≈ 5.9 is the hidden gauge coupling, and V µ is the vector meson matrix The covariant derivatives of ξ have the usual simple transformation property under the hidden gauge symmetry: The simplest operator that preserves HLS and the global chiral symmetry is where in the last line we are showing only a single term from the 1/f π expansion of ξ L,R . We estimated the Wilson coefficient by considering the analogous term without the leptonic factor, which contains vector meson masses m 2 V ∼ g 2 V f 2 π if C lequ = 1. Since the decays ρ + → π + π 0 and ρ 0 → π + π − have O(1) branching fractions, these operators can easily give rise to multi-track signals, while avoiding the phase space and chiral (1/f π ) suppression of the multi-body decays involving only mesons (if the vectors can be produced on shell).

Comparison with Pythia
We have implemented the hadronization of the scalar semileptonic operator in Eq. (22) using Pythia, and by incorporating the chiral-perturbation theory and HLS results described above in a FeynRules model. In the latter approach we include terms in the 1/f π expansion that enable χ → νe + (π − , π − η, π − ππ, πρ) and allowing the heavier mesons to decay as η → π 0 π + π − and ρ → ππ; depending on the charges of ρ and π these decay channels can produce up to four tracks per vertex. Note that there is an ambiguity in the relative branching fractions to final states involving vector mesons, since the Wilson coefficients in Eqs. (25) and (34) are only fixed by inexact arguments. We have generated full event samples using both hadronization approaches and computed the signal acceptance using the selections discussed in Sec. 5. In Fig. 16, we compare acceptances for a single point in the LLP parameter space with m A = 8 GeV, and cτ χ ≈ 1 cm. In the left panel we show the N dv ≥ 2 and N tr ≥ 2 analysis which illustrates remarkable agreement between the two approaches. In this case, the dominant decay channel is χ → νe + π − which is already captured by the leading term in the chiral expansion. Significant differences emerge in the N dv ≥ 2 and N tr ≥ 3 analysis shown in the right panel of Fig. 16; we see that the chiral perturbation/HLS treatment of hadronization leads to a smaller acceptance for m χ 3 GeV compared to Pythia. We can understand this discrepancy as follows. In decays mediated by Eq. (22), about half of the rest mass energy of χ flows into the hadronic part of the operator; demanding ≥ 3 tracks further partitions this energy into a number of particles with total mass of at least 3m π , suggesting that the phase space suppression from finite meson masses becomes important. We have confirmed that this an important effect by evaluating partial widths to different hadronic final states at unphysically low meson masses in the chiral perturbation/HLS approach. 9 Additionally, the invariant mass distribution of the ud pair is broad, and even in the m χ ∼ 3 GeV case half of the events have √ s ud 1.2 GeV. These observations suggest that the Pythia hadronization model most likely overestimates the number of high-multiplicity events, leading to a higher signal acceptance for GeV-scale LLP masses. The chiral perturbation theory/HLS approach, however, is also deficient because we do not know how to precisely fix the relative sizes of different operators that contribute to the decay. It is reassuring that both methods agree for larger LLP masses m χ 3 GeV.
B.3ū R u RūR u R Operator Fully hadronic decays can be modelled in ChiPT/HLS similarly to the semi-leptonic case. We begin with the operator O uuuu = χ(q R C uu γ µ q R )(q R C uu γ µ q R ), whose EW-symmetric form is given in Eq. (14); C uu is a spurion of SU (3) L × SU (3) R which must have the expectation value to generate the χ(ū R u R )(ū R u R ) operator. Chiral invariance demands that this matrix transforms as C uu → RC uu R † .
Since there are no SU (3) L transformations involved, we can build our operators out of objects charged only under SU (3) R , such as ξ † R D µ ξ R which transforms as These objects are also invariant under the HLS. Some invariant combinations are As before, we only highlighted a few representative terms, but there are many more that can be important. Of particular interest to us are those containing η, ω 0 and other mesons that can decay to multiple charged particles. Unfortunately we do not know a way to fix the relative sizes of these operators, so the branching fractions of the different hadronic states are uncertain. For concreteness, we take the coefficient of the operator leading to χ → π + π − to be four times larger than the others, since this leads to a better agreement with Pythia in the N tr ≥ 2 analysis as we show below.

Comparison with Pythia
We have implemented the hadronization of the scalar hadronic operator in Eq. (36) using Pythia, and by incorporating the ChiPT/HLS results described above in a FeynRules model. In the latter approach we include terms with up to three meson fields and up to O(1/f 3 π ); the leading decays are into ππ, π 0 η, πππ, ρπ, ω 0 π 0 and others. We note again that the branching fractions are uncertain since multiple operators with unknown Wilson coefficients contribute to different channels; this is already apparent for χ → π 0 π 0 versus χ → π + π − , which are generated by operators in Eq. (40a) and (40b).
In Fig. 17 we show a comparison of the signal acceptance for χ →ūuūu hadronized in Pythia and using the ChiPT/HLS approach outlined above; we focus on a single benchmark point with m A = 8 GeV and cτ ≈ 1 cm. The results for the N dv ≥ 2, N tr ≥ 2 (N dv ≥ 2, N tr ≥ 3) analysis are given in the left (right) panel. In the N tr ≥ 2 case the signal is dominated by χ → π + π − decays and we see that the acceptances agree between the two hadronization methods within a factor of ∼ 2. This agreement can be improved by tuning the Wilson coefficients of the operators contributing to different branching fractions; we leave this exercise to future work. As for semi-leptonic decays, the N tr ≥ 3 exhibits significant differences between the Pythia and ChiPT/HLS approaches, with Pythia predicting a significantly better acceptance at lower LLP masses. We hypothesize again that this is due to the string fragmentation model employed by Pythia not assigning a large enough phase space suppression for producing multiple mesons. Thus, high-multiplicity events are much more common in the Pythia event samples compared to ChiPT/HLS, leading to a significantly larger acceptance at small LLP masses. At the same time, our ChiPT treatment doesn't include broad vector resonances beyond the lightest multiplet which (as in Sec. B.1), which could lead to an under-prediction of higher multiplicity final states in the ChiPT model.