A Theory of Dark Pions

We present a complete model of a dark QCD sector with light dark pions, broadly motivated by hidden naturalness arguments. The dark quarks couple to the Standard Model via irrelevant $Z$- and Higgs-portal operators, which encode the low-energy effects of TeV-scale fermions interacting through Yukawa couplings with the Higgs field. The dark pions, depending on their $CP$ properties, behave as either composite axion-like particles (ALPs) mixing with the $Z$ or scalars mixing with the Higgs. The dark pion lifetimes fall naturally in the most interesting region for present and proposed searches for long-lived particles, at the LHC and beyond. This is demonstrated by studying in detail three benchmark scenarios for the symmetries and structure of the theory. Within a coherent framework, we analyze and compare the GeV-scale signatures of flavor-changing meson decays to dark pions, the weak-scale decays of $Z$ and Higgs bosons to hidden hadrons, and the TeV-scale signals of the ultraviolet theory. New constraints are derived from $B$ decays at CMS and from $Z$-initiated dark showers at LHCb, focusing on the displaced dimuon signature. We also emphasize the strong potential sensitivity of ATLAS and CMS to dark shower signals with large multiplicities and long lifetimes of the dark pions. As a key part of our phenomenological study, we perform a new data-driven calculation of the decays of a light ALP to exclusive hadronic Standard Model final states. The results are provided in a general form, applicable to any model with arbitrary flavor-diagonal couplings of the ALP to fermions.


Contents
A Decays of a light ALP coupled to Standard Model fermions 34 A.1 a → γγ 36 A.2 a → π + π − γ 37 A.3 a → π + π − π 0 38 A.4 a → 3π 0 40 A.5 a → π 0 π 0 η, π + π − η 40 A.6 a → π 0 π 0 η , π + π − η 41 A.7 a → ηηπ 0 41 A.8 a → K 0 K 0 π 0 41 A.9 a → K + K − π 0 42 A.10 a → K + K 0 π − , K − K 0 π + 43 A.11 a → ωω, φφ, K * + K * − , K * 0 K 1 Introduction and the model A light, confining hidden sector coupled feebly to the Standard Model (SM) is in general an interesting possibility for new physics, often referred to as a hidden valley (HV) [1]. More sharply, it can be part of the answers to outstanding questions of the SM. The (little) hierarchy problem may be solved by models of neutral naturalness [2][3][4], where the partners of the top quark are not charged under SM color but a dark color symmetry, and dark confinement around the GeV scale is a generic prediction [5]. If the hierarchy problem is solved by cosmological relaxation, a confining hidden sector may be the origin of the backreaction potential that stops the relaxion [6]. Dark strong dynamics can also provide attractive scenarios for dark matter, with several plausible candidates found among the dark hadrons.
If at least some of the dark hadrons decay to SM particles, the feeble coupling connecting the hidden and visible sectors generally implies macroscopic lifetimes. Thus, hallmark signatures of HV models are given by long-lived particles (LLPs), which have been a topic of rapidly increasing interest at the Large Hadron Collider (LHC) [7] and beyond [8]. At the LHC, searches for LLPs hold a strong discovery potential, provided that dedicated and innovative strategies can be implemented at the level of event selection and analysis. This is especially true for "dark jet" or "dark shower" topologies, where the decay of a heavy particle (such as a Z or Higgs boson) to the hidden sector produces jets made of light dark hadrons. The associated phenomenology deserves further attention. 1 Areas where important progress is needed include maximizing the dark shower coverage of existing detectors (primarily ATLAS, CMS and LHCb), understanding the interplay with low-energy production processes such as flavor-changing neutral current (FCNC) meson decays, as well as comparing to the sensitivity of proposed LLP-specific experiments. In this paper we study a model of dark QCD with light pseudo Nambu-Goldstone bosons (pNGBs), namely dark pions, coupled to the SM through irrelevant Z and Higgs portals. This theory of dark pions provides a new coherent framework to address the above questions.
The low-energy spectrum of the hidden (or dark) sector depends on the number of light quark flavors, N , charged under the dark QCD (assumed to have SU (N d ) as the color group) and having masses below the strong scale Λ. If N = 0, dark glueballs are at the bottom of the spectrum [10]. An example is the Fraternal Twin Higgs model [5], where the lightest dark glueball is expected to mix with the Higgs boson, giving rise to phenomenology that has been extensively studied [5,11]. For N = 1, the low-energy spectrum contains several mesons with masses around Λ, the lightest being a (η -like) pseudoscalar, a vector, and a scalar [12]. For instance, this scenario has been thoroughly analyzed in a realization of the tripled top framework for supersymmetric neutral naturalness [13] that features electroweak-charged top partners, where the dark mesons mix dominantly with the Z boson [14]. When N ≥ 2 (but still below the conformal window, N 4N d ), one expects chiral symmetry breaking and N 2 −1 associated pNGBs, which in a slight abuse of notation we call dark pions,π, for any N . As familiar from the SM, the dark pions can be much lighter than the rest of the hadrons, whose masses are at or above the dark QCD scale: mπ Λ. Here we focus on the multi-flavor case. The lifetimes of the dark pions depend on the amount and pattern of explicit isospin breaking, yielding a larger parameter space to explore compared to the one-flavor theory. With respect to the latter, notable differences are that dark meson production and decay are less tightly connected, and dark vector mesons dominantly decay to dark pions if the phase space is open, whereas for N = 1 they decay to SM particles.
The light dark quarks must be singlets under the SM if their masses are below O(100) GeV, to satisfy collider bounds. The interactions connecting them to the visible sector dictate the phenomenology. Here we focus on the interesting possibility of irrelevant portals obtained by integrating out heavy states [14][15][16][17][18][19][20], which have been less studied compared to renormalizable ones, even though they have solid theoretical motivations. A concrete example is given by the scenario of Ref. [14], where the mediation is provided by heavy fermions charged under both the SM electroweak (EW) and hidden color gauge symmetries, allowing for renormalizable Yukawa interactions between dark-colored quarks and the SM Higgs doublet. The supersymmetric partners of the heavy fermions play the role of scalar top partners, hence the mediation scale is naturally around TeV.
The model. Drawing from the above discussion, the theory of dark pions considered in this work contains N > 1 flavors of Dirac fermions ψ i , transforming in the fundamental representation of the dark color SU (N d ), but singlets under all SM gauge symmetries. In addition, N EW-doublet Dirac fermions Q i = (Q u Q d ) T i with hypercharge 1/2 are included, which also transform in the fundamental representation of SU (N d ). This field content allows for Yukawa couplings involving the SM Higgs doublet, where Y , Y , M , and ω are N × N matrices in flavor space. The mass matrices M and ω can be diagonalized with real and positive diagonal elements by separate unitary transformations on the Q L,R and ψ L,R fields, respectively, so we assume this form without loss of generality. The coupling matrices Y , Y can be complex in general, with N 2 +(N −1) 2 independent phases, decreasing to (N − 1) 2 if one of Y , Y , or ω vanishes. In addition there is always the strong CP phase of dark QCD, which will be consistently set to zero in this work. The masses M of the heavy dark quarks are taken to be larger than Λ (and around TeV). The masses of the light dark quarks, which receive independent contributions of order ω and Y Y v 2 /M , where v is the Higgs vacuum expectation value (VEV), are assumed to be much smaller than Λ. Hence, the dark QCD has N light flavors. If N 4N d [21], at low energies the light quarks are confined and form a condensate. The SU (N ) L × SU (N ) R chiral symmetry is spontaneously broken to the diagonal SU (N ) V , resulting in N 2 − 1 pNGB dark pions.
The Y = y t 1, Y = 0 limit of Eq. (1.1), together with appropriate TeV-scale supersymmetry breaking, embodies a solution to the little hierarchy problemà la tripled top [13,14]. 2 This sets a well-motivated target for the chiral structure and mediation strength, with Y ∼ y t ∼ 1 and M ∼ TeV. On the other hand, the same Lagrangian (1.1) was employed in the relaxion solution to the hierarchy problem [6], 3 to generate a backreaction potential without running into difficulties with the SM strong CP problem. In both cases N = 1 was originally chosen for the sake of minimality, but not necessity. The fact that this setup emerges naturally in very different approaches to the Higgs naturalness puzzle makes it a compelling choice for a benchmark theory of dark pions.
The properties of the individual dark pions depend on the symmetries and structure of the dark sector. If the CP symmetry is preserved, the dark pions are classified into odd and even states: for example, with N = 2 theπ 1 andπ 3 are CP -odd (J P C = 0 −+ ) while theπ 2 is CP -even (0 −− ), where the index corresponds to SU (2) generators. Therefore, this theory provides a coherent framework to study both CP -odd and -even light scalars feebly coupled to the SM. The CP -odd dark pions decay to SM particles through the Z portal, i.e., by mixing with the longitudinal component of the Z boson through dimension-6 operators. They behave as ALPs with an effective decay constant parametrically given by where fπ is the dark pion decay constant, defined in analogy with the SM pion decay constant f π ≈ 93 MeV. The precise form of Eq. (1.2) is derived later, but one can already see that for Y ∼ 1, M ∼ TeV and fπ ∼ GeV, the CP -odd dark pions have f a ∼ PeV. This highlights how the ALP decay constant does not necessarily correspond to a physical scale (no threshold exists near the PeV in our model), but is a combination of parameters of the underlying theory if the ALP is composite. The CP -even states decay to SM particles through the Higgs portal, i.e., by mixing with the Higgs boson through dimension-5 operators. As we show in detail later, the mixing angle is parametrically where m h is the Higgs boson mass. Since small dark pion masses are well motivated in our setup, and the dark pions couple to all SM fermions including quarks, the phenomenological analysis presented here requires a detailed description of ALP decays for an ALP mass m a 3 GeV, where exclusive hadronic SM final states must be considered. We obtain this by means of a novel calculation that extends the data-driven methods proposed in Ref. [24]. We emphasize that the results, reported in Appendix A, apply to any ALP with arbitrary flavor-diagonal couplings to SM fermions. For the decays of light scalars, we make direct use of previous calculations [25]. As can be gleaned from Eqs. (1.2) and (1.3), dark pion theories with Y ∼ Y and with Y Y are very different, because in the former the dimension-5 Higgs portal dominates, whereas in the latter the dimension-6 Z portal is most important. In addition, in general the dark sector contains CP violation, which leads to mixing of different states and induces couplings of all dark pions to both Z and Higgs portals, with relative strengths determined 2 Effective theory for dark quarks Starting from the Lagrangian (1.1) and assuming GeV, we can integrate out the Q fields at tree level and obtain the EFT 4 In a broader perspective, see also studies of heavier dark pions with masses above the EW scale [26,27].
where we have retained operators up to dimension 6. In general, this effective Lagrangian contains the same number of complex phases that appear in the UV, except if either Y or Y vanish, in which case the counting is reduced to (N − 1)(N − 2)/2 [the "apparently missing" phases then appear in additional operators that were not included in Eq. (2.1)]. The first terms in square brackets in the first two lines of Eq. (2.1) renormalize the dark quark kinetic terms after inserting the Higgs VEV. These small corrections are neglected in the following, unless otherwise noted. The second terms in square brackets in the first two lines generate interactions of the ψ with the Z boson. The third line gives rise to the mass matrix, where the last term is induced by the seesaw mechanism. For general Y and Y the mass eigenstates ψ are obtained via unitary transformations ψ L,R = U L,R ψ L,R , and their diagonal mass matrix is Barring cancellations, the ψ i are light if both terms in Eq. (2. 2) are small compared to Λ. This occurs most naturally if there is an (approximate) chiral symmetry acting on ψ L (or ψ R ) to suppress both ω and Y (or Y ). For example, that is the case in the tripled top model, where Y = 0 [13,14]. The third line of Eq. (2.1) also generates the leading couplings of the dark quarks to the Higgs.

Constraints from Z and Higgs invisible decays
The first, important constraints on the parameter space are obtained by assuming that the dark hadrons mostly go undetected at colliders, so that the bounds on the Z invisible width from LEP and on the Higgs invisible width from LHC apply. The EFT in Eq. (2.1) induces Z decays to dark quarks via dimension-6 operators, where the small dark quark masses were neglected. For M = M 1, this gives a branching ratio The LEP measurement of the Z invisible width requires ∆Γ inv Z < 2 MeV at 95% CL [38], and from Eq. (2.4) we obtain If Y ∼ Y parametrically, the leading interaction of the dark sector with the Higgs boson is the dimension-5 operator in the third line of Eq. (2.1), yielding For M = M 1 the associated branching ratio is where we have taken Y ∼ Y ∼ 0.1 as reference value for the Yukawas. Satisfying the current invisible Higgs width constraint, BR(h → inv) < 0.13 at 95% CL [39], requires Note that for Y ∼ Y ∼ 1 the bound is M 40 TeV. The above Z, h → invisible bounds are applied widely in the rest of the paper, as we focus mainly on GeV-scale dark pions, for which assuming invisible dark jets is a reasonable first approximation. Nonetheless, it should be kept in mind that these bounds may be weakened or lifted in regions of parameter space where most dark pions are short lived.
A quick glance at Eqs. (2.6) and (2.9) indicates that the product Y Y is much more severely constrained than Y 2 or Y 2 . Given a new physics scale M , scenarios where Y ∼ Y parametrically are subject to a coupling constraint about one order of magnitude stronger than scenarios with Y or Y ∼ 0. This will have an important impact on the phenomenology, as the dark pion lifetimes scale with the fourth power of the Yukawa couplings. In addition, if Y ∼ 0 [Y ∼ 0] the dominant decay of the Higgs to the dark sector is either toĝĝ via the one-loop operator where c Q = 1 arises from integrating out the Q at one loop, 5 or to ψ ψ via the first [second] line of Eq. (2.1). These processes are too suppressed to lead to current constraints from h → invisible, but are discussed here for completeness. The one-loop Higgs decay to dark gluons gives resulting in a branching ratio for M = M 1, (2.12) If Y = 0 we obtain, using the dark quark equation of motion (EOM), and for ω ≈ m ψ 1 and M = M 1, the branching ratio is (2.14) In the opposite case Y = 0, one replaces Y → Y and ω → ω † in Eq. (2.13). 5 More generally, cQ may receive contributions from additional states, e.g. scalars in the tripled top model [13,14].

Indirect constraints
At one loop, integrating out the heavy fermions Q in Eq. (1.1) also generates higherdimensional operators built only of SM fields, which can be subject to relevant constraints.
The most important one is (H † ↔ D µ H) 2 , encoding a contribution to the EW T parameter [40,41]. In fact T is most easily calculated in the UV theory, by applying, e.g., the results of Ref. [42]. The derivation of a general analytical expression is rather cumbersome, but the calculation simplifies if the dark Yukawas are diagonal, Y = diag i y i and Y = diag iỹi : at leading order in the large -M i expansion and taking real couplings for simplicity. The general case including flavor mixing can be treated numerically in a straightforward manner. It is useful to compare the T parameter and Z → invisible constraints, in the simple scenario where for the former we have used the rough estimate T 10 −3 and the latter follows from Eq. (2.6). Since the two are comparable for Y ∼ O(1), and additional beyond-SM contributions can a priori alter the interpretation of the T constraint, in most of our discussion we stick to the more robust invisible Z width bound. When Y ∼ Y , both are subleading to the invisible h branching ratio constraint.
The operators |H| 2 B µν B µν and |H| 2 W i µν W µν i are also generated at one loop. However, since the Q d are electrically neutral and the Q u are charged but do not couple to the Higgs, we expect the operators to come in the linear combination |H| 2 (g 2 W i µν W µν i − g 2 B µν B µν ) which gives a vanishing contribution to the hγγ coupling.
CP violation in the dark sector could feed into the visible sector, inducing electric dipole moments (EDMs) for SM particles. The strongest limit comes from the electron EDM [43]. Corrections to the electron EDM arise through the loop-suppressed operator O B B = |H| 2 B µν B µν , which in turn contributes at one loop to the EDM (similar considerations apply to |H| 2 W i µν W µν i ). Inspection of the relevant diagrams shows that O B B does not arise at one loop. Furthermore, if Y = 0 or Y = 0 the two-loop contributions turn out to be strongly suppressed by an extra ∼ ω 2 /M 2 factor. If both Y and Y are non-vanishing [44]. This is much weaker than the Higgs invisible branching ratio bound, M 40 TeV Y Y from Eq. (2.9). In summary, we find that EDMs do not provide additional constraints in this model.

Effective theory for dark hadrons
At energies below Λ, the SU (N d ) gauge group confines and the dark quarks and gluons form hadrons. For N ≥ 2, the lightest hadrons are pNGBs of the SU (N ) L ×SU (N ) R → SU (N ) V symmetry breaking and belong to the adjoint representation of SU (N ) V . As the simplest example and representative case for phenomenological studies, in this work we focus on N = 2. The three dark pionsπ a are defined in the basis where the light quark mass matrix is diagonal,π where σ a are the Pauli matrices. Importantly, theπ 2 has J P C = 0 −− whereasπ 1,3 have 0 −+ , as can be derived from Eq. (3.1) using ψ 2 P L,R ψ 1 C → ψ 1 P L,R ψ 2 . Note that in the absence of a U (1) flavor symmetry, i.e., if Y , Y are not diagonal,π 1 andπ 2 are distinct states. Their degeneracy will be lifted by Y , Y interactions, as demonstrated later by explicit examples. We do not discuss in detail the dark flavor-singletη , which at small N d receives a large mass from the dark U (1) A anomaly.
The couplings of the dark pions to the Z boson can be derived from the interactions in the dark quark EFT of Eq. (2.1), where g Z = g 2 + g 2 . We rewrite this as 3) where the dimensionless matrices A and A are defined as and σ 0 ≡ 1 2 . In addition, The pions are excited by the axial vector current. We define their decay constant from with normalization corresponding to f π ≈ 93 MeV in the SM. Thus the last term on the right-hand side of Eq. (3.3) yields a tree-levelπ a -Z mixing, and the partial width for the decay to a pair of SM fermions f is where a f = T 3 Lf and N f c = 3 (1) for quarks (leptons). It is important to note that, in the absence of CP -violating phases,π 1,3 decay through the single Z exchange butπ 2 does not, because where we have used the hermiticity of A, A. For light dark pions it is in fact convenient to integrate out the Z boson, obtaining Due to the conservation of the vector current, the interaction relevant to describe dark pion decays is a is the effective decay constant ofπ b . Equation (3.10) enables us to apply the new calculations presented in Appendix A, where for arbitrary (flavor-diagonal) ALP-SM fermion couplings we perform the matching to the SM chiral Lagrangian, augmented with exchange of scalar, vector, and tensor resonances above 1 GeV, and by extending datadriven methods pioneered in Ref. [24] we evaluate the ALP decay widths to an extensive set of exclusive hadronic SM final states. The results are reported in Fig. 1, which is one of the main novelties of this work. The lifetime is also shown in Fig. 8, see Appendix A.
To gain some initial insight on the scales we take, e.g., Y = 0, giving the parametric scaling where CP conservation was assumed for simplicity. As the constraint from Z → invisible gives roughly M/Y TeV, see Eq. (2.6), for fπ ∼ GeV the CP -odd dark pions can be regarded as light ALPs with effective decay constants PeV.
In fact, from the quark-level EFT in Eq. (2.1) we can directly derive that the dark pions couple to the Higgs current, namely a . 6 Because the dark pions are appropriate degrees of freedom only at energies below Λ and the latter is smaller than the EW scale in most of our parameter space, the use of the broken EW phase is warranted and such an effective description is not fully justified. Nonetheless it affords us a first brief discussion of FCNC meson decays [45][46][47], by applying the leading-log results of Ref. [48]. The flavor-changing couplings of the dark pions to quarks arise at one loop, Note that the appropriate mass scale to cut off the logarithm is M ∼ TeV -the largest physical threshold here -and not f a , which is a combination of parameters with dimension of a VEV and does not correspond to the mass of any particle. In addition, owing to the modest separation between M and m t , finite pieces are expected to be important. Both expectations are confirmed by the explicit calculation in Section 5. There, we show that current meson FCNC constraints are at the level f a ∼ 10 3 TeV, as obtained from B → X sπ decays (where X s denotes a strange hadron state) with long-livedπ → µ + µ − at CHARM, LHCb and CMS for mπ 2m µ , and from searches for K + → π +π with invisibleπ at E949 and NA62 for smaller dark pion masses.
The dark pions can also decay through tree-level Higgs exchange. To derive the decay width, the starting point are the following interactions in Eq. (2.1), where we have already rotated to the quark mass eigenstate basis and the coupling matrix B is dimensionless. The piece of Eq. (3.13) containing γ 5 is relevant for dark pion decay, and we rewrite it as Finally, recalling Eq.
which allows us to calculate the decay width mediated by a single Higgs exchange, where we have employed the relation m 2 πa =B 0 Tr(m ψ ), valid at leading order in the dark sector chiral perturbation theory (ChPT), see Appendix B. It is immediate to see that if CP is conserved, the trace in Eq. (3.16) can be non-vanishing only for a = 2, since (iσ 2 ) * = iσ 2 whereas (iσ 1,3 ) * = −iσ 1,3 . Note that the interference between the Z-and h-mediated amplitudes vanishes in theπ a → ff process.
Comparing Eqs. (3.7) and (3. TeV andB 0 ∼ 10 GeV (dimensionally, we expectB 0 ∼ 4πfπ). However, for hierarchical Yukawas the ratio is suppressed by Y 2 /Y 2 or viceversa, and the pions decaying via the Z mediation can have much shorter lifetimes than those decaying via the Higgs exchange.
For GeV-scale dark pions, Higgs-mediated decays to exclusive hadronic SM final states become important. We do not attempt to reassess them here, but account for them following the results of Ref. [25] (see also Ref. [49] for a recent reappraisal), by matching to their definition of the couplings As in Ref. [25] we take m s = 95 MeV, however we include the running of m c,b in the perturbative spectator model and consider the decay to photons, The matching constant parametrizing φ → 4π, ηη, . . . [25] is fixed to C ≈ 4.8 × 10 −9 GeV −2 . The resulting decay widths and branching ratios are shown in Fig. 2.
Parametrically, the dark pion-Higgs mixing angle takes on the scaling, where, as in Eq. (3.11), CP conservation was assumed for simplicity. Since the bound from h → invisible reads roughly Y Y /M 0.03 TeV −1 (see Eq. (2.9)), for fπ ∼ GeV the CP -even dark pions can be viewed as Higgs-mixed scalars with mixing angles 10 −6 . When a single dark pion decays through both the Z and h portals in the presence of CP violation, we neglect the interference between the two amplitudes, which vanishes forπ a → ff as already noted, but can a priori be nonzero for more complex final states. In this work we focus on the mass range mπ 2m b ∼ 10 GeV, where a wider range of experiments are relevant and our results are expected to have the most impact.
We end this section with some brief comments on the heavier dark hadrons, including non-pNGB mesons and baryons. The dark vector (and axial-vector) mesons may be relevant to intensity frontier phenomenology, where they can be produced first and subsequently decay to dark pions, if kinematically allowed. In the N = 2, N d = 3 theory considered here, lattice QCD calculations at pion masses larger than their physical values can be repurposed [50] to parametrize the hidden sector, at least for moderately heavy pNGBs with 0.1 m 2 π /m 2 V 0.7, whereV denotes the dark vector resonance. As for the baryons, the lightest among them is stable due to dark U (1) B , but its relic density can easily be very suppressed unless a dark baryon asymmetry is present. In this work we focus on the properties and phenomenology of the dark pions, neglecting the heavier hadrons.

Benchmark scenarios for dark pions
In this section we discuss the range of possibilities for the dark pion properties, beginning with general arguments. If the theory respects the SU (2) V isospin symmetry, i.e., ω, M , Y , Y ∝ 1 2 , then A, A, B ∝ 1 2 and all dark pions are stable. It is also possible that the SU (2) V is explicitly broken to its U (1) subgroup, i.e., ω, M , Y , Y and hence A, A, B are diagonal. In this caseπ ± ≡ (π 1 ∓ iπ 2 )/ √ 2 is charged under the U (1) and therefore stable, whileπ 0 ≡π 3 can decay (to avoid any confusion, we remark that the subscripts do not indicate SM electric charge).
Furthermore, specific models can give rise to distinctive patterns for the masses and couplings. For example, in a setup inspired by the tripled top [14], but where the two Scenario Symmetries possessed Decay portals Table 1: Summary of the benchmark scenarios considered in this work. A chiral symmetry can be the origin of Y = 0. In the second scenario,π ± = (π 1 ∓ iπ 2 )/ √ 2 is stable because it is the lightest particle charged under a dark U (1). hidden sectors share a common dark color gauge group, we expect M = M 1 2 , Y y t 1 2 and Y 0, which implies A, B 0. In this case the main source of isospin breaking in the hidden sector comes from the diagonal m ψ ω. As a result, the U (1) subgroup is approximately preserved and theπ 1,2 have much longer lifetimes thanπ 3 .
The above considerations make it clear that, even for the minimal dark pion theory with N = 2, the parameter space is too vast to be covered systematically in this first study. Therefore we choose to discuss a few benchmark scenarios that give rise to distinct phenomenology. With these, our aim is to be illustrative rather than exhaustive, and we expect that other interesting patterns may be found in future work. We begin with a few comments on the case of stable dark pions and their possible role as dark matter, and then turn to the study of three benchmark scenarios where at least some of the pions are unstable and decay to SM particles. The key features of these three are summarized in Table 1. For each scenario, in the phenomenological analysis we fix generic textures for the Yukawa and mass matrices, paying attention to avoiding enhanced symmetry points. This reduces the number of independent parameters to a manageable handful.
Scenario 0: Isospin-symmetric limit and dark pion dark matter As already mentioned, for ω, M , Y , Y ∝ 1 2 the dark pions form a stable triplet of SU (2) V , which is a dark matter candidate. 7 However, in this limit the dark pions do not couple to the Z (see Appendix B), hence reducing their cosmological abundance to a viable level requires adding extra ingredients to the theory. For N ≥ 3 the number density can be depleted via 3 → 2 processes mediated by the Wess-Zumino-Witten action, potentially realizing Strongly Interacting Massive Particle (SIMP) dark matter [31], although an additional mediation between the dark pions and the SM should still be introduced to transfer the dark matter entropy to the SM. If the mediator is a dark photon that mixes kinetically with the hypercharge [51][52][53], care must be taken to check dangerous decays of singlet pions, which can be made viable through appropriate mass splittings for odd N [53,54], or prevented by imposing suitable discrete symmetries for even N [54]. Such scenarios provide appealing origins for light thermal dark matter, but as they are rather tangential to the central aspects of this work, we do not discuss them further.

Scenario 1: Y = 0
In this case there is no constraint from the invisible decay branching ratio of the Higgs. In general Y contains 1 physical phase, which can be parametrized, e.g., as with real y ij . It is convenient to perform a further redefinition which renders and the same for ψ L1 so that the mass matrix remains real. 9 All three dark pions are unstable. As anticipated,π 1 andπ 3 have unsuppressed decay to SM particles via the Z portal, so their lifetimes and branching ratios can be directly obtained from Appendix A. Since Y = 0, instead of the Higgs portalπ 2 decays through CP -violating mixing with the other pions. To estimate its lifetime, we need to take into account several corrections to the leading-order pion Lagrangian: • The pion mass splitting generated by O(p 4 ) ChPT operators with insertions of the quark mass matrix, e.g., where c 7 is a coefficient expected to be of O(1) by naive dimensional analysis. For generic dark isospin breaking |ω is the leading correction to the pion masses. Therefore, to estimate the CP -violating decay ofπ 2 we can focus only on its mixing withπ 1 . 10 • The effects of tree-level Z exchange, which correct the kinetic term ofπ 1 (andπ 3 ), but notπ 2 . 8 We assume cos α y11y12/M 2 1 + y21y22/M 2 2 > 0 for definiteness. 9 The low-energy quark masses are given by as obtained after applying the leading-order EOM to the first term on the right-hand side of Eq. (2.1). This is diagonalized to m ψ by ψL,R = UL,Rψ L,R , but in practice we neglect the (Y v/M ) -suppressed corrections. 10 Notice that if ω1 = ω2 the theory actually preserves CP , because the phase α can be removed by a • The one-loop contributions from box diagrams, with parametric scaling These yield in particular a CP -violating mixing ofπ 1 andπ 2 , provided is nonvanishing. This is the case if α = 0, all y ij = 0, and Once the above effects are included, the kinetic terms forπ 1,2 are made canonical by the rotation .
(4.7) To understand quantitatively the dark pion properties we focus on the following pattern for the Y and M matrices, which is of generic nature. From Eqs. (3.10) and (4.7) we find the effective decay constants ofπ 1,3 and the CP -violatingπ 1 -π 2 mixing angle, respectively, The decay width of the physicalπ 2 is then Γπ 2 ≈ sin 2 θ 12 Γπ 1 . Dark pion decays are mainly controlled by the three parameters y/M, fπ, and mπ. The mediation strength is constrained by the Z invisible width: Eq. (2.4) gives y/M 1.1 TeV −1 (assuming N d = 3). The dark pion lifetimes are shown in the left panel of Fig. 3, choosing y/M that saturates the LEP bound and fixing fπ = 1 GeV; other results are obtained by rescaling τ ∝ f −2 π (y/M ) −4 , see Eq. (4.9). Remarkably, for y ∼ 1, M ∼ 1 TeV and fπ ∼ 1 GeV, i.e. parameter choices motivated by (neutral) naturalness [14], the lifetime ofπ 1 falls between 10 meters and 1 millimeter across the mass range 2m µ mπ 2m b . Therefore, this dark pion is a natural LLP target for present and future experiments. On the other hand,π 2 andπ 3 have much longer lifetimes. As θ 12 depends on M but not on y, for illustration we show the range of τπ 2 obtained by varying M ∈ [1.1, 5] TeV, where the lower edge corresponds to the current bound on the Q mass from direct searches at the LHC (see Section 7). We stress that for theπ 2 lifetime we have performed an estimate, rather than a precise calculation, as sufficient for our purpose. The right panel of Fig. 3 shows selected branching ratios, which are the same for the three dark pions as they all decay through the Z portal (if the CP -violating mixing with the other pions is very suppressed,π 2 may decay through Higgs mediation via a small Y = 0, but we do not study that possibility here).
This benchmark scenario provides a theoretically motivated and remarkably simple target for current and future experimental probes. The constraints from and future opportunities in FCNC meson decays are discussed in Section 5, whereas the prospects for discovery at the LHC via Z decays to dark showers are presented in Section 6.

Scenario 2: exact U (1)
If the U (1) symmetry {ψ, Q} → e ixσ 3 {ψ, Q} is preserved, the Yukawa matrices are diagonal. Parametrizing the two physical phases as Y = diag (y 1 , y 2 ), Y = diag (ỹ 1 e iα 1 ,ỹ 2 e iα 2 ), the EFT quark mass matrix (2.2) is diagonal but complex, and is transformed into a real and positive m ψ by the rephasings 11 which also leave the (real) Zψψ coupling matrix Y † M −2 Y unaffected. While the charged pionπ ± is stable,π 0 decays through the Z portal and, in the presence of CP violation, the Higgs portal. The cosmological history can be easily safe. The mass splitting of charged and neutral pions is controlled by the operator in Eq. (4.3), which yields mπ 0 < mπ + if c 7 < 0. Thenπ +π− →π 0π0 conversions followed by decays ofπ 0 to the SM with lifetime τπ 0 To simplify the analysis of the parameter space we assume ω i = y iỹi v 2 M i cos α i , which can be regarded as a particularly simple case of the scenario where ω and YỸ v 2 /M are of the same order. This choice leads toα i = α i and m ψ i = y iỹi v 2 /(2M i ).
In addition, we choose the generic patterns giving the bounds from invisible Higgs and Z decays where we have setB 0 = 4πfπ.
Owing to the CP violation, the decays ofπ 0 are an intricate combination of Z-and h-mediated processes. For the former, the effective decay constant is found from Eq. (3.10), (4.14) For the latter, the coupling of the dark quarks to the Higgs reads Tr .
(4.15) The first equality holds in general and shows that s (0) θ would vanish for ω = 0, as a consequence of sin(α i +α i ) = 0. At this stage we can take {yỹ/M, M, r, mπ} as the four independent free parameters, with fπ fixed via Eq. (4.13). In Fig. 4, we set yỹ/M to its upper bound from h → invisible and explore the remaining three-dimensional parameter space. As expected, theπ 0 lifetime depends strongly on r: if M ∼ O(TeV), for r 1 or 1 the Z portal dominates (with branching ratios that are well described by Fig. 1). Note that the dependence of the dark pion lifetime on its mass is very different from scenario 1, as can be observed by comparing with Fig. 3. The reason is that, while in scenario 1 fπ is independent from mπ, here Eq. (4.13) dictates the scaling fπ ∝ m 2 π , resulting in a much shorter lifetime as the dark pion mass increases.
Conversely, for r ∼ 1 the Higgs portal plays an important role, dominating the total width for mπ 2 -3 GeV. For this reason, in the bottom panels of Fig. 4 we show the branching ratios at r = 1, which best illustrate the complexity of the decay pattern. If mπ 2 GeV the branching ratio to the CP -even KK final state is of several percent, which could rise up to ∼ 15% in the region below the cc threshold, although the description adopted here [25], based on the ss final state in a perturbative spectator model, does not permit a more accurate prediction. For mπ above the cc threshold Higgs exchange completely dominates the width, yielding the interesting prediction that a heavier (and therefore shorter-lived)π 0 mainly decays to CP -even final states if the dark Yukawa interactions contain sizable CP violation.

Scenario 3: exact CP
The third and last scenario we consider is one where CP is exactly preserved by the dark Yukawa interactions. As in scenarios 1 and 2, to reduce the number of independent parameters the Yukawa and mass matrices are set to definite patterns. These are chosen to be of generic nature, avoiding points of enhanced symmetry. We take .

(4.19)
In the top left panel of Fig. 5 we show theπ 2 lifetime as a function of mπ, for several values of κ. τπ 2 becomes very long for κ 1, because in the limit κ → 0 we have B ∝ m ψ which is diagonal, hence d(κ) → 0 . τπ 2 also increases for κ 1, due to the larger c(κ) which for fixed mπ requires a smaller value of fπ, thereby suppressing s (2) θ . The shortest lifetime for a given mπ is thus obtained for κ ∼ 1, i.e., when the mass scales ω and Y Y v 2 /M are close.
The CP -odd dark pionsπ 1,3 decay only via the Z portal, with decay constants that depend strongly on r = y/ỹ , where p (b) , q (b) are dimensionless functions. It is instructive to compare the lifetimes of all three pions. To do so we focus on κ = 1, showing in the top right panel of Fig. 5 the lifetimes for illustrative values of r. 12 At small masses the CP -even pionπ 2 has the longest lifetime irrespective of r, but for mπ 4 (6) GeV it becomes the shortest-lived for r ∼ 1 (10). Recalling the branching ratio patterns shown in Figs. 1 and 2 for CP -odd and -even pions, we conclude that the expected signatures from dark shower events display a striking dependence on mπ. For simplicity, here we have considered r ≥ 1; in the opposite regime r ≤ 1 the behavior is very similar, but with the roles ofπ 1 andπ 3 reversed: in particular, for r 1 it isπ 3 that has the shortest lifetime. Thus far, we have fixed yỹ/M to the upper bound from h → invisible. If this parameter is decreased by a factor n > 1, fπ must be correspondingly increased by n in order to keep the same dark pion mass, as dictated by Eq. (4.18). These two effects exactly compensate (for fixed M ) in the decay constants f we read that the net effect on s (2) θ is an n-fold increase, and therefore theπ 2 lifetime becomes n 2 times shorter. We illustrate this somewhat counter-intuitive effect in the bottom panel of Fig. 5, which shows that even for mπ 2 m c,τ theπ 2 lifetime can be as short as O(1 -10) m, provided fπ ∼ 10 GeV. Furthermore,π 2 can easily have the smallest lifetime among the dark pions. These results are especially interesting in view of a proposed LHCb search for LLPs decaying to K + K − [55], which may have sensitivity to ourπ 2 since in this mass region its BR to KK is sizable, see Fig. 2. As for the largest plausible value of fπ, the neutral naturalness framework suggests Λ ∼ 4πfπ 100 GeV, corresponding to fπ 10 GeV, while the ultimate limit is Λ M , otherwise, the Q's cannot be treated as heavy dark quarks anymore, the global symmetry pattern is modified and the EFT breaks down.

FCNC meson decays
Light dark pions may be produced in FCNC meson decays if kinematically allowed. To describe these decay rates, we calculate the four-fermion effective operators of the form d Lα d Lβ ψ ψ with α < β. In our theory they arise through two classes of one-loop diagrams: Z exchange with insertion of thed Lα d Lβ Z coupling, and box diagrams containing W and Q u internal lines. The amplitudes can be fully obtained from the classic results for ds → νν in Ref. [56], leading to where x q ≡ m 2 q /m 2 W and y k ≡ M 2 k /m 2 W (recall that the mass of Q k u is simply M k ). For our purposes we can safely take the large-y k limit of D, where only the first few terms of the dominant top loop were retained. The meson decay amplitude is then, assuming factorization of the hadronic matrix elements into a SM factor and a hidden factor, 4) where we have focused on B → Xπ a decays with X = K, K * , and applied Eq. (3.6). For the decay widths we find . The log-enhanced contribution to Γ(B → Kπ a ) is in agreement with what one finds [48] from Eq. (3.12), but the finite terms have an important quantitative impact: for M = 1 TeV, retaining only the logarithmic piece overestimates the rate by a factor ≈ 3. The definitions and numerical values of the form factors f 0 , A 0 are taken from the light-cone QCD sum rules analysis of Ref. [57], with f B→K 0 (0) ≈ 0.27 and A B→K 0 (0) ≈ 0.31. An expression analogous to the first line in Eq. (5.5) applies to K → ππ a , with the appropriate replacements of masses, CKM elements, and the form factors available from lattice QCD with f K→π 0 (0) ≈ 0.97 [58]. 13 FCNC decays can also produce the CP -even dark pions, through Higgs mixing. The corresponding amplitudes are proportional to the Higgs penguin, resulting in [60,61] and for the decay widths. Evaluating Eqs. where in the CP -odd case we have set M = 1 TeV in the logarithm.

Constraints and projected sensitivity
We now highlight a few implications for our parameter space, focusing mainly on mπ > 2m µ . The theoretical predictions in Eq. (5.8) can be compared with the current BaBar [62] and Belle [63] 90% CL bounds on invisible decays, For CP -odd scalars, branching ratios at the 10 −5 level require f (b) a < 100 TeV, but in this regime the dark pion lifetimes become sufficiently short to ensure that decays to SM particles occur inside the detector (see Fig. 1 or 8), thus violating the search assumptions. Therefore more relevant are searches for B → K ( * ) (χ → µµ) with long-lived (scalar or pseudoscalar) χ at LHCb [64,65], as well as the re-interpretation in terms of these decays [66] of results from the CHARM beam dump experiment [67]. In addition, CMS has recently presented a novel search based on data scouting [68], setting limits on the inclusive branching ratio for B → X s (χ → µµ) [69]. In our setup this may be related to the exclusive branching ratios via as estimated from the observed values of BR(b → s ) and BR(B → K ( * ) ) [70]. The sizable uncertainty reflects the still-unsettled experimental status of these measurements. The relation (5.10) enables a direct comparison of the CMS and LHCb/CHARM bounds. In Fig. 6 we show such comparison for four representative ALP masses in the range m a 2m c , where searches for a → µ + µ − are relevant, as seen from the branching ratios in Fig. 1. The LHCb and CHARM constraints are taken from Ref. [66], whereas we apply here for the first time the CMS bound [68] with the help of Eqs. (5.10) and (5.5). For each value of m a , CMS provides limits for τ = 1, 10, 100 mm, corresponding to the red points in the (f a , BR) plane of Fig. 6; we simply interpolate between those points and include the uncertainty band arising from the relation between inclusive and exclusive branching ratios. 14 Figure 6 shows a clear pattern: for low (high) mass the strongest constraint comes from CHARM (CMS), while in the intermediate region LHCb has the best sensitivity. To better estimate the bounds as functions of the dark pion mass, we combine the above results at fixed m a with the findings of the Expression of Interest for CODEX-b [71], where constraints on f a for an ALP coupled universally to SM fermions were reported following the analysis of Ref. [66], but with updated lifetime and branching ratio calculations employing data-driven methods [24]. The main differences between our setup (where the ALP couples to weak isospin) and the universal coupling scenario [60,66,71] are the ALP total width and the treatment of finite terms in the B → Ka calculation. For the former, a detailed comparison in Fig. 8 (right panel) shows qualitative agreement, although important quantitative differences are present; for the latter, in Refs. [60,66] only the leading-log term was retained and the cutoff was set to 1 TeV, which combined with slightly different values for the form factors gives a rate ∼ 4 times larger than here. In light of these considerations we apply the f a bounds for universal couplings [71] to our setup, after weakening them by a factor ∼ 2 to account for the smaller production rate. Where relevant, the resulting estimates agree with Fig. 6.
For 2m µ m a 0.6 GeV, the re-interpretation [66,71] of CHARM results gives the strongest constraint. 15 In this region we estimate f a 1.3 -1.9 PeV, 16 translating in 14 In the top right panel of Fig. 6 we actually use the CMS bound for ma = 610 MeV, as 600 MeV is masked in the analysis [68]. We neglect the impact of this small difference. 15 In the universal coupling scenario [71] the CHARM bound was found to extend up to ma ∼ 1 GeV, but a direct comparison shows that in our setup it is limited to ∼ 600 MeV, see Fig. 6. 16 Here we quote the lower limit on fa from CHARM, but note that a small "wedge" of allowed fa may remain between the LHCb, CHARM and CMS exclusions for 0.3 ma/GeV 0.6, see the top panels of Fig. 6. when applied toπ 1 using Eq. (4.9), and taking conservatively the weakest bound in the given mass range. Thus, for fπ 2 GeV the CHARM sensitivity surpasses Z → invisible. Considering inclusive decays would likely strengthen the CHARM bounds compared to those used here, which were derived from B → K ( * ) a only [66]. Here again we have been conservative, adopting the weakest bound in this mass range; close to the upper end, the constraint can actually be about twice as strong, as seen in the bottom right panel of Fig. 6. It should be emphasized that a theoretical uncertainty affects this bound, stemming from Eq. (5.10). Looking ahead, several proposed LLP experiments at the LHC have the potential to improve the sensitivity on ALPs coupled to fermions in the mass range 2m µ m a 2m c , including CODEX-b [71], FASER 2 [72], and MATHUSLA [73]. Importantly, in contrast to current bounds that rely on a → µ + µ − , these experiments would be sensitive to any decays to ≥ 2 charged tracks and therefore to a → π + π − π 0 , which in our model dominates between 1 and 3 GeV (see Fig. 1). As already discussed above for LHCb and CHARM constraints, we can roughly estimate the projected sensitivities from the results for universal ALP-fermion couplings in Ref. [71]. Caveats concern the total ALP width, as shown in Fig. 8, and the production rate, which is assumed to arise dominantly from FCNC B meson decays. For example, in our setup mixing with π 0 , η, η may enhance the production, owing to the non-trivial U (3) transformation properties of the ALP. With these disclaimers, we obtain for m a = 1 GeV the projections f a 10 PeV at FASER 2, f a 20 PeV at CODEXb, and f a 80 PeV at MATHUSLA200. In addition, for SHiP with 10 20 protons on target we find f a 14 PeV [66], based on the a → µ + µ − signature.
The decays B → K a, a → hadrons with m a in the GeV range have also been studied as probes of a heavy QCD axion, where the dominant coupling to the SM is aG G. Both prompt a → π + π − π 0 , ηπ + π − , KKπ, φφ [74] and displaced a → π + π − π 0 [75] have been considered and projections for Belle II obtained. Our branching ratio calculations in Appendix A can serve as the basis to extend those results to the class of models where the ALP couples dominantly to SM fermions.
Finally, for decays to CP -even dark pions we find from Eq. (5.7) hence NA62 [76] is probing mixing angles of O(10 −4 ). Comparing this with the expectation in benchmark scenario 3, where we have set mπ ∼ 250 MeV and κ = 1 (recall that the product fπ × yỹ/M is then fixed by Eq. (4.18)), suggests that FCNC meson decays to CP -even dark pions are out of experimental reach, unless one is willing to consider an extreme hierarchy between fπ and mπ, with the former exceeding the TeV.

Z -initiated, muon-rich dark showers at the LHC
In the previous section we have discussed processes at energies well below the weak scale, where the dark pion properties can be fully described through the low-energy parameter combinations f a and s θ , for CP -odd and -even states respectively. Here we take a step up in energy and consider production of dark pions via Z and Higgs decays to dark partons, followed by showering and hadronization. As we are going to show, these processes access new directions in parameter space compared to FCNC meson decays. The LHC inclusive production cross sections for Z and Higgs bosons are (see, e.g., Ref. [14]) σ(pp → Z) ≈ 54.5 (58.9) nb, σ(pp → h) ≈ 48.6 (54.7) pb, (6.1) at 13 (14) TeV. The coupling structure of our model implies that Z decays dominate in scenarios with Y or Y ∼ 0, whereas h decays are most important if Y ∼ Y , as quantified by the branching ratios to dark quarks in Eqs. (2.5) and (2.8). Here we focus on Z decays to the dark sector, which are largely unexplored but hold a strong LHC discovery potential, as the forthcoming discussion illustrates. The Z → ψ ψ decay results in two dark jets, dominantly composed of dark pions with high multiplicity. GeV-scale dark pions eventually decay to a variety of SM final states, as seen in Figs. 1 and 2. For mπ 2m c , the FCNC meson decays discussed in Section 5.1 set a lower bound f a O(PeV), implying in turn a lower bound on the dark pion lifetimes. Concretely, in scenario 1 with mπ ∼ 1 GeV we obtain from Eq. (5.12) a constraint τπ 1 > O(1 -10) cm, which sets the target for dark shower searches in this mass range. Differently from scenarios with t-channel mediation such as emerging jets [16,81], here the signal is not automatically accompanied by hard SM jet activity, hence the trigger strategy is a central issue. For this reason in the first exploration we focus onπ → µ + µ − decays, which result in striking displaced vertices (DVs) at the LHC and a narrow resonance peak that can be exploited to suppress the combinatorial and misidentification backgrounds [82]. 18 The sensitivity of LHCb to dark shower signals is well established [14,82] (see also a recent overview [83]) and the most recent search for dimuon resonances [84] has already provided a HV interpretation. Building on these results, in Section 6.1 we perform a detailed recast to set current bounds and estimate the future reach of LHCb on our Zinitiated, muon-rich dark shower signals. By contrast, for ATLAS and CMS we limit ourselves to some qualitative comments in Section 6.2, whereas a detailed study is deferred to a separate publication due to its more complex nature [85] (see also Refs. [9,[86][87][88][89] for discussions of other dark shower signals).

LHCb sensitivity
We base our reinterpretation on the latest LHCb search for displaced dimuons [84]. We generate pp → Z → ψ ψ at 13 TeV using the HV module of Pythia8 [90][91][92], with the production cross section in Eq. (6.1) as normalization. To set the dark pion parameters we focus on benchmark scenario 1 (Section 4.1), where all three dark pions decay through the Z portal, considering two mass points with the following characteristics: Here Nπ is the average number of dark pions per dark jet. As three different lifetimes cannot be accommodated by the HV module, we neglect the longest-livedπ 2 (which is also subject to larger uncertainties) and fix the ratio τπ 3 /τπ 1 ≈ 37 as expected from Eq. (4.9). This leaves τπ 1 and BR(Z → ψ ψ ) as free parameters of our analysis.
To derive the current constraint, we apply at truth level the displaced search cuts listed in Table 1 of Ref. [84] and compare to the cross section limits for promptly-produced X → µ + µ − [84] (this is the appropriate choice, as dark parton shower and hadronization are prompt in our model, and we require the reconstructed X to come from the primary vertex). We find the p X T ∈ [5, 10] GeV bin dominates the sensitivity, resulting in the solid black exclusion curves in Fig. 7. The right minimum of the exclusion contours corresponds to optimal sensitivity to theπ 1 signal with τπ 1 ∼ few mm, whereas the left minimum corresponds to optimal sensitivity to decays ofπ 3 , with τπ 3 ≈ 37τπ 1 ∼ few mm.
To estimate the future reach, we follow a slightly different strategy: we calculate the signal rate after cuts and parametrize remaining detector effects through a DV efficiency µµ that is varied in the range [0.4, 0.8]. This is compared to the background rate extracted from Fig. 2 in Ref. [84], which is found to be ≈ 1.6 (≈ 0.7) events per 5.1 fb −1 for the mπ = 650 MeV (1 GeV) hypothesis, by averaging over the m µ + µ − ∈ [600, 700] MeV ([0.9, 1.1] GeV) window and considering a bump-search interval |m µ + µ − − 650 MeV (1 GeV)| < 2σ with σ being the experimental resolution. When applied to the 18 Hadronicπ decays are alternative opportunities, especially when the final states are fully charged: for example,π → K * 0 K * 0 → (K + π − )(K − π + ) through the Z portal, orπ → K + K − through the h portal.
The phenomenology of these hadronic final states within dark showers deserves future study.  Figure 7: Projection of the 90% CL LHCb sensitivity [84] to Z-initiated, muon-rich dark showers for mπ = 650 MeV (top) and mπ = 1 GeV (bottom). The two minima of the exclusion contours correspond to optimal sensitivity to decays of two dark pion species with different lifetimes,π 1 andπ 3 , while decays of the longest-livedπ 2 are neglected. The current exclusion is shown by the black curve, while the widths of all other bands are obtained by varying the single-DV efficiency µµ ∈ [0.4, 0.8]. Brown lines indicate the relation between BR(Z → ψ ψ ) and τπ 1 obtained from benchmark scenario 1, for representative choices of fπ.
current luminosity, this procedure gives the dashed gray bands in Fig. 7. The reasonable agreement with the actual LHCb constraint gives us confidence in the method, which is then applied to Run 3 (23 fb −1 ) and High-Luminosity LHC (HL-LHC, 300 fb −1 ) scenarios to obtain the red and blue bands. For mπ = 650 MeV, LHCb will probe Z branching ratios down to ∼ 10 −7 in the high-luminosity phase, with further improvements possible either through optimization to the dark shower signal or future detector upgrades. The reach for mπ = 1 GeV is somewhat weaker, due to the lower dark pion multiplicity and smaller dimuon branching ratio.
The brown lines in Fig. 7 show the relation between BR(Z → ψ ψ ) and τπ 1 that is realized in benchmark scenario 1, as a function of fπ. The dependence on the underlying parameters should be contrasted with complementary bounds from other processes, namely Z → invisible, which probes y/M , and B decays, sensitive to f a ∝ M 2 /(y 2 fπ). We learn that for 1 fπ/GeV 20 the LHCb dark shower search has already probed new parameter space, highlighting the strongly complementary role of this type of analysis with current and upcoming data.
In addition to the single-DV analysis we consider requiring 2 DVs per event, assuming zero background in this case. The corresponding exclusions, shown by the orange and green bands in Fig. 7, turn out to be weaker than the single-DV ones. This is explained by the fact that the background is already very suppressed for 1 DV, hence removing it completely results in a limited gain, and by the additional efficiency cost.
The potential sensitivity of LHCb to heavier pseudoscalars, with masses above a few GeV, has also been discussed in several final states [93].

ATLAS and CMS prospects
In the light of the results shown in Fig. 7, and in particular the correlation observed in our framework between BR(Z → ψ ψ ) and the dark pion lifetimes, a priori ATLAS and CMS may lead to dramatic improvements in the region τπ ∼ 0.1 -1 m, thanks to their larger volumes (and integrated luminosities). However, owing to the soft nature of the signals considered here, progress requires targeted experimental strategies that enable efficient triggering on low-p T displaced muons.
A major step in this direction has recently been achieved by CMS with the search for dimuon DVs [68] in data collected with scouting triggers, which permit the unprecedented exploration of very low muon transverse momenta and thus DV masses, down to the m µµ ∼ 2m µ threshold. This approach is well suited to test theπ → µµ signals discussed here, as demonstrated by the new constraints on the parameter space we have derived in Section 5 from the CMS B → X s (χ → µµ) results [68]. Thus a recast to the dark shower signal is warranted, which will be presented elsewhere [85]. We note that the CMS analysis imposes a cut l xy < 11 cm on the transverse displacement of the dimuon DVs, due to the definition of the scouting trigger stream which requires hits in at least two pixel layers. Looking ahead to future upgrades, CMS-specific triggers targeting LLP dimuon signals have also been proposed [94].
At ATLAS, a search for two "dark photon jets" [95] targeted final states related to those of interest here: Higgs decays to two jet-like structures, each composed of an invisible particle and two GeV-mass LLPs decaying to µ + µ − . Events were selected by means of a trigger requiring ≥ 3 L1 muons with p T > 6 GeV, then confirmed at HLT using only muon spectrometer information. It results in optimal sensitivity for O(cm) lifetimes. Compared to the signal model used by ATLAS, our Z-initiated dark shower has larger multiplicity, lower transverse momenta, and for mπ 1 GeV larger branching fraction to muons.
Heavier LLPs have been searched for in a number of analyses by ATLAS and CMS, mainly focusing on rare Higgs decays to the hidden sector, see e.g., Refs. [96,97] for very recent results.

Probing the ultraviolet completion
Finally, we take another step up in energy and discuss the expected LHC signals from direct production of the heavy fermions Q. Since these carry SM EW charges, they undergo Drell-Yan (DY) pair production such as, for instance, ud → W + * → Q u Q d . The decay patterns can be read from the Yukawa interactions by means of the Goldstone equivalence theorem: Q u → W + ψ, whereas Q d decays to Zψ and hψ with ≈ 1/2 branching ratios for M m Z,h . Flavor indices have been suppressed for simplicity. The ψψ pair in the final state give rise to two dark jets, which characterize the signal.
Assuming the dark pions are sufficiently long-lived to escape the detector we obtain W Z/W h + MET, a typical signature of EWinos in supersymmetry. Similar considerations apply to the production of the electrically neutral pairs Q u Q u , Q d Q d . Consequently, bounds on M can be directly set by applying the results of searches for Higgsinos, which are assumed to decay directly to the lightest supersymmetric particle (LSP), taken to be the bino-like neutralino. Our signal matches this topology in the limit of very light neutralino LSP. The strongest sensitivity has been achieved, remarkably, in the all-hadronic + MET search by ATLAS [98], which outperforms analogous searches for 3 + MET and bb + MET. For degenerate Higgsinosχ and massless bino LSP, a bound mχ > 900 GeV (95% CL) was obtained. Our signal cross section reads at partonic level We have assumed the Q i are not mass-degenerate, which applies to all benchmark models considered in Section 4 (for two degenerate Q i , the constraint strengthens to 1.2 TeV). Given the current Higgsino expected bound [98], we rescale the cross section by L/L with L, L = 139, 3000 fb −1 and derive M 1.3 TeV as our estimate of the (13 TeV) HL-LHC sensitivity.
If the dark pions are heavy enough to decay inside the detector, the phenomenology becomes similar to the emerging jets scenario [16,100], albeit with EW rather than QCD production of the mediators. Evaluating the impact of the existing CMS search [81] on our signals is beyond the scope of this work, and left as an interesting avenue for future studies.
Beside DY production we consider single Q production mediated by off-shell Higgs, gg → h * → Q d ψ. This yields Z/h + ψψ final states, leading to mono-Z/h signatures if the dark pions escape undetected. The partonic cross section is found to be, neglecting the ψ mass, .4), and again we have neglected flavor indices. By folding in the gg parton luminosity we obtain the 13 TeV cross sections 9.7 ab} , (7.5) where the renormalization and factorization scales were set to M . These results show that single production cannot compete in rate with DY, though the sensitivity to the Yukawa couplings makes it a complementary probe of the UV completion.

Conclusions
In this paper we have formulated a theory and initiated the study of dark pions, coupled to the SM via irrelevant Z and Higgs portals. The corresponding operators are obtained by integrating out TeV-scale EW-doublet fermion mediators. This setup has strong UV motivations, appearing in various modern approaches to the hierarchy problem, such as neutral naturalness models and the relaxion scenario. It provides a concrete framework where the GeV-scale phenomenology of the dark pions, the EW-scale decays of Z and h bosons to the hidden sector, and the TeV-scale signals of the mediators are all coherently linked.
The decays of CP -odd and CP -even dark pions proceed via tree-level mixing with the Z and h, respectively, providing explicit realizations of light composite ALPs and scalars coupled feebly to the SM. For CP -odd dark pions, we have provided a new comprehensive calculation of the decay widths to exclusive hadronic SM final states, obtained by applying data-driven methods. The results are valid for any ALP with arbitrary flavor-diagonal couplings to SM fermions, and can therefore be widely used to study other models.
The dark pion phenomenology depends on the symmetries possessed by the model, including CP , dark isospin, and chiral symmetries. To illustrate the range of possibilities we have analyzed in detail three benchmark scenarios. We find that for masses and couplings of the mediators that can be related to the hierarchy problem while satisfying experimental constraints, and for dark pion decay constants around the GeV scale, dark pions with 2m µ mπ 2m b have lifetimes varying from a millimeter to 10 meters. Intriguingly, this is the most interesting range for LLP searches at the LHC (and beyond), making the dark pions a natural target. We have begun the exploration of the signatures with two applications, meson FCNC decays and Z-initiated dark shower searches, focusing primarily on the mass region 2m µ mπ 2m c where the striking dark pion decay to dimuons has a significant branching ratio.
Searches for flavor-changing b → sa decays, with long-lived a → µ + µ − , set important bounds on the effective decay constant of the CP -odd dark pions, f a PeV. In addition to well-known constraints from CHARM and LHCb, we have derived new ones from a recent CMS search leveraging the data scouting technique. Each of these experiments turns out to have the strongest sensitivity in a different mπ window. Proposed LLP detectors at the LHC interaction points, including FASER 2, CODEX-b, and MATHUSLA, have the potential to extend the sensitivity on f a by 1 -2 orders of magnitude. For mπ < 2m µ , there are lower bounds f a PeV from K → π + invisible searches at E949 and NA62. On the other hand, the CP -even dark pions remain out of reach due to their very small mixing with the Higgs.
Dark shower searches at the LHC access the additional structure that partially completes the theory at the EW scale. They probe decays of on-shell Z and h bosons to dark jets composed mainly of long-lived dark pions. Z decays to the dark sector, in particular, have been largely overlooked so far, but here we have shown that they probe new directions in the parameter space, supplying orthogonal information to meson FCNC decays. We have performed a thorough recast of the most recent LHCb search for displaced dimuons. The resulting constraints demonstrate that the sensitivity to Z-initiated dark showers has already reached new parameter space, surpassing competing bounds from meson FCNC and Z → invisible decays. ATLAS and CMS have strong potential to extend the reach to longer dark pion lifetimes, which are well motivated in our framework, by exploiting larger decay volumes and luminosities. Dedicated experimental strategies are increasingly being implemented, such as data scouting/trigger-level analysis, and a detailed assessment of their impact on our framework will appear elsewhere.
As for the direct LHC reach on the EW-charged mediators, a straightforward reinterpretation of Higgsino searches in all hadronic + MET final states gives the constraint M 1.1 TeV. The improvement expected in the high-luminosity phase is mild, leaving open the possibility that a dark pion discovery may take place at the LHC, while the direct production of the mediators would need to wait for a future collider.
Looking ahead, many paths deserve further exploration. Hadronic decays of GeVscale dark pions are shown to be important by our results, warranting new studies both for FCNC meson decays and dark shower searches at the LHC. Notable modes include:π → π + π − π 0 , which we find to dominate the width of light CP -odd dark pions,π → K * 0 K * 0 → (K + π − )(K − π + ) andπ → K + K − , which can be fully reconstructed and have sizable branching ratios in some parameter regions, and several others discussed in Sections 3 and 4. The sensitivity of Belle II to such modes requires detailed studies, as well. In addition, we have not touched upon the heavier mass range mπ 2m c , where hadronic decays dominate and lifetimes become significantly shorter. In particular, it would be interesting to understand if in this region there are any constraints on the EW pair-production of the heavy dark quark mediators from the existing CMS search for emerging jets.
The dark pion phenomenology at fixed-target experiments also remains to be investigated. We note that dark hadrons heavier than the dark pions may be relevant there, due to different production mechanisms which could be exploited to test specific regions of parameter space. For instance, bremsstrahlung production of dark vector mesons can be strongly enhanced if their mass is around 1 GeV, due to mixing with SM vector meson resonances.
Finally, the sensitivity of future colliders to the scenario presented here warrants further studies. In particular, an e + e − machine like FCC-ee would offer extraordinary possibilities to probe decays to the hidden sector at a Tera-Z phase, as it has already been demonstrated for one-flavor dark QCD models. We believe the present work sets a solid foundation to tackle all the above aspects, while providing several new results of general applicability in the study of light, feebly coupled hidden sectors. ES thanks the organizers of the MIAPP "Novel Hidden Sectors" and Portorož "Physics of the flavorful Universe" workshops for kind invitations to present preliminary results of this work.

A Decays of a light ALP coupled to Standard Model fermions
The starting point is the Lagrangian with f ∈ {q, , ν} for quarks, charged leptons and neutrinos. The width for decay to a pair of charged leptons is If the ALP is much heavier than the SM QCD scale, Λ SM m a , its hadronic decays can be analyzed perturbatively. The width for decay to two gluons is [101,102] where n q counts the quarks lighter than m a , while the loop function is 19 We have B 1 (x) ≈ 1 (− 1 3x ) for x 1 ( 1), implying that light quarks contribute ≈ c q /(32π 2 ) to the sum in Eq. (A.3) whereas heavy quarks rapidly decouple. For decay to heavy quarks Q = c, b, where m Q is the running quark mass in the MS scheme. We use two-loop running for both α s and m c,b , and set m c = 1.67 GeV, m b = 4.78 GeV.
For m a Λ SM we must consider decays to exclusive hadronic final states instead. To do so we match Eq. (A.1) to the low-energy effective Lagrangian [24,103,104], The pseudoscalar matrix is written as where f π ≈ 93 MeV. The hard U (1) A breaking due to the anomaly is parametrized by m 2 0 and the physical η, η are related to the octet and singlet fields by This approximate value of the mixing angle is sufficiently accurate for our purpose, while simplifying analytical expressions [24]. The relevant pieces of the Lagrangian describing the vector resonances are [105] The structure of B1 can be understood upon integrating by parts the interaction with quarks in Eq. (A.1) and using the expression of the divergence of the axial current, where V µν = 1 2 µνρσ V ρσ (with 0123 = 1) and g V V P = −N c g 2 /(8π 2 f π ) is determined by the anomaly. The coupling g is fixed by the Kawarabayashi-Suzuki-Riazuddin-Fayyazuddin (KSRF) relation [106,107] to g = g V ππ = m V /( √ 2f π ) ≈ 6.0, where for m V we have taken the ρ mass. 20 The vector meson matrix reads It is important to note that Eq. (A.9) realizes vector meson dominance (VMD) for π 0 → γγ but retains an (anomalous) γP 3 contact interaction, with coefficient equal to −1/2 of the one given by the WZW action. This choice was shown to provide a better fit to data compared to "complete VMD" [105], and will impact the calculation of a → π + π − γ.
We then assign to the ALP the U (3) representation a = 1 13) which are taken to be valid up to m a ≈ 3 GeV. Above this mass we switch to the perturbative description. The model studied in this paper has c f = T 3 Lf i.e. c u = −c d = −c s = 1/2, giving K aπ 0 = 1, K aη = 1/ √ 6 and K aη = −1/ √ 3. However, we stress that our results are general and also apply to other models with different patterns of ALP-fermion couplings, for example those in Refs. [19,108]. We are now in the position to calculate the decays of low-mass ALPs to exclusive final states.

A.1 a → γγ
We begin with the decay to two photons [24,101], (A.14) 20 This g should not be confused with the SU (2)L gauge coupling, which never appears in this appendix.
where C γ is defined by the effective operator Cγ α 8πfa a µνρσ F µν F ρσ . The individual contributions are where the ρ 0 , ω, φ matrices are implicitly defined by Eq. (A.10) and m * a is the scale where the VMD and pQCD terms are matched, which equals ≈ 2.9 GeV for our benchmark model. The form factor F ≡ F 4 accounts for the suppression of the V V P interaction at high mass, extracted in Ref. [24] by comparison to e + e − data, A basic cross-check of Eq. (A.14) is that, setting a → π 0 and f a → f π , it reproduces the classic result for Γ(π 0 → γγ), which in the VMD picture is mediated by π 0 ρ 0 ω = 1/2. In addition, the predicted widths for η, η → γγ match the experimental values within 20%.
A.2 a → π + π − γ The amplitude is described by 5 diagrams: two with ρ 0 exchange, two with ρ ± exchange, and one contact interaction. For the spin-summed squared matrix element we find and, adopting a convention we follow consistently, the final-state particles were ordered according to how we define the decay (i.e. in a → π + π − γ, 1 denotes the π + and so on). In addition, the four-momenta satisfy p a = i ∈ final p i . The width is ) 1/2 + m 2 a + 2m 2 π − m 2 12 and the symmetry factor S = 1. We cross-check this result by applying it to the η : setting a → η , η P = δ η P and f a → f π gives Γ(η → π + π − γ) ≈ 56 keV, in excellent agreement with the PDG value of 55 keV. The same procedure applied to the η yields Γ(η → π + π − γ) ≈ 90 eV, to be compared with the experimental value of 55 eV. Our η prediction would get significantly closer to the observed rate if we used a more precise value of θ ηη and accounted for the SU (3)-breaking differences among the pseudoscalar decay constants [109], which however go beyond the scope of this work. Nonetheless, we remark that the γP 3 contact term in the vector meson Lagrangian (A.9) is important to improve agreement with data: omitting this term (as done, e.g., in Ref. [24]) we obtain 154 eV (63 keV) for η ( ) → π + π − γ, so the η partial width is off by a factor ≈ 3 relative to the observed value. In the numerics we actually replace the quantity in square parentheses with its expression including isospin breaking up to O(δ I ), where δ I ≡ (m d − m u )/(m d + m u ), as provided in Eq. (S32) of Ref. [24]. A k factor equal to 2.7 is included, derived from comparison with η ( ) → 3π data [24]. where F V ≡ F 3 . The first two pieces arise from ρ ± exchange diagrams, while the third one originates from the ∂ 2 P 4 /f 2 π interactions in Eq. (A.9) and is essential to ensure that M VMD vanishes at low energy, as can be explicitly verified by taking BW ρ m −2 ρ and applying the KSRF relation. Exchange of the σ scalar meson yields M σ = −2γ 2 σππ f π f a aπ 0 p a · p 3 p 1 · p 2 BW σ (m 2 12 )Θ(4m 2 K − m 2 12 )F(m a ) , (A. 25) where γ σππ = 7.27 GeV −1 , as well as all the couplings of the scalar nonet mesons that appear in the following, are taken from the fit to data performed in Ref. [110] without assuming U (3) symmetry (we use the second set of couplings given in Ref. [110]). To 21 In Ref. [24] only the diagram containing an aρ0ρ0 vertex was included, which corresponds to retaining only the piece proportional to √ 6 aη + √ 3 aη BWρ 0 (m 2 12 ) = 6 aρ0ρ0 BWρ 0 (m 2 12 ) in the second line of Eq. (A.20). Upon integrating this partial amplitude over dm 2 23 we find agreement with Ref. [24].

A.13 Comparison with previous studies
The main predecessor in the study of light ALP hadronic decays is Ref. [24], with which our analysis shares several important aspects. In particular, we adopt their choice of vertex form factors in Eq. (A.19) to suppress the resonance exchange amplitudes at large m a . There are, however, some major differences that we wish to summarize here: • The key distinction is that, as we consider scenarios where the couplings to SM fermions dominate, in general the ALP has a non-trivial U (3) representation for all masses up to ≈ 3 GeV (where we match to perturbative QCD). By contrast, Ref. [24] focused on the case where the coupling to gluons dominates, therefore C u = C d = C s was assumed for m a 1 GeV. The nontrivial U (3) representation of ALPs with mass above 1 GeV implies that here the a → P (V → P P ) decays are in general unsuppressed and play a crucial role. This is clearly demonstrated by our benchmark model c f = T 3 Lf , where a → π ± (ρ ∓ → π ∓ π 0 ) dominates not only the a → π + π − π 0 amplitude, but also the total ALP width for m a 1 GeV, as shown in Fig. 1. This is a consequence of the sizable ALP mixing with π 0 and the strong coupling g V ππ ≈ 6. Other effects of the nontrivial U (3) charges include strong relative suppressions for certain channels, such as e.g., Γ 3π 0 Γ π + π − π 0 and Γ K * + K * − Γ K * 0 K * 0 (see Fig. 1).
• We do not assume U (3) invariance to determine the scalar nonet contributions to a → 3P decays, as the results of Ref. [110] show this to be a rather poor approximation. Instead, we make use of all the couplings fitted to data in Ref. [110], taking into account all relevant a -P mixings. As a result, our amplitudes for scalar mediation agree in kinematic structure with Ref. [24], but differ in the values of the couplings.
• For the tensor meson f 2 we assume U (3) invariance with f 2 = diag (1, 1, 0)/2 and determine the g f 2 ππ coupling from data, as in Ref. [24]. However, we differ from that reference in that we use the unitary gauge propagator for the massive spin-2 field, leading to corrections to the f 2 contribution to a → 3P amplitudes. In addition, we fix the coefficient of the g µν piece in the ∂Σ † ∂Σf 2 interaction (this piece does not enter the calculation of on-shell f 2 → ππ, so its coefficient has to be fixed from other considerations) to the value corresponding to f 2 coupled to the energy-momentum tensor [113], see Eq. (A.26). Finally, we turn off the f 2 exchange amplitudes for m 2 ij < (m f 2 −Γ f 2 ) 2 , to avoid unphysical contributions to the O(p 4 ) terms in the chiral Lagrangian. The impact of different prescriptions for the f 2 couplings and propagator is shown in the left panel of Fig. 8, considering for illustration the a → π + π − π 0 decay.
• Other differences compared to Ref. [24] are described above for each process. These include a different treatment -with several new contributions -for a → π + π − γ and the addition of further decay channels such as a → ηηπ 0 and a → (ρ 0 → π + π − ) ω.
The ALP lifetime for the scenario with universal couplings to fermions, derived from the methods of Ref. [24], has also appeared before in the literature [71]. In the right panel of Fig. 8 we compare it to our determination for c f = T 3 Lf . While the results are qualitatively compatible, important quantitative differences appear for m a ∼ m π and in the region m η m a 2m c . Lf , but different prescriptions for the f 2 couplings and propagator lead to different expressions for M f 2 . In solid orange, the choice made in this work: f 2 couples to the energy-momentum tensor and has a unitary gauge propagator, leading to Eq. (A.29). In dotted blue (dashed purple), alternative versions where the coefficient of the g µν piece in Eq. (A.26) is set to −1/4 (0), still with unitary gauge propagator. In dot-dashed green, the version used in Ref. [24] where the propagator has the Landau gauge expression, i.e. in Eq. (A.30) one replaces m 2 f 2 → k 2 . For this choice, which does not seem justified, the result is independent of the coefficient of the g µν piece in the coupling. (Right) total lifetime obtained from our calculation with ALPfermion couplings proportional to weak isospin, compared to the lifetime for universal couplings [71].

B Chiral perturbation theory for dark pions
At energies below the scale of resonances, the dark pions are described using ChPT. To lowest order for N = 2, where U is the pion matrix transforming as U → L U R † under SU (2) L × SU (2) R , m ψ is the generalized quark mass matrix containing also the interactions with the Higgs, andB 0 is a non-perturbative constant that determines the dark pion masses, where the form of m ψ follows from Eq. (3.13). According to Eq. (3.3), the covariant derivative of U takes the form The above equations allow us to derive, in particular, the linear mixing between the dark pions and the Z, L π ⊃ −g Z fπTr[σ a (A − A)]∂ µπa Z µ /4, and the linear mixing between the dark pions and the h, L π ⊃B 0 fπTr[iσ a (B − B † )]π a h/2, both of which are of course consistent with the current algebra results given in Section 3. If SU (2) V is exact and therefore A, A, B ∝ 1 2 , all interactions of the dark pions with the Z in Eq. (B.1) vanish. For the single -Z terms this is a consequence of Tr(U † ∂ µ U ) = 0, valid for any N (see e.g., Ref. [117]).