Strongly interacting dark sectors in the early Universe and at the LHC through a simplified portal

We study the cosmology and LHC phenomenology of a consistent strongly interacting dark sector coupled to Standard Model particles through a generic vector mediator. We lay out the requirements for the model to be cosmologically viable, identify annihilations into dark vector mesons as the dominant dark matter freeze-out process and discuss bounds from direct detection. At the LHC the model predicts dark showers, which can give rise to semi-visible jets or displaced vertices. Existing searches for di-jet resonances and for missing energy mostly probe the parameter regions where prompt decays are expected and constrain our model despite not being optimised for dark showers. We also estimate the sensitivity of dedicated analyses for semi-visible jets and emphasize the complementarity of different search strategies.


Introduction
A rapidly growing effort is being dedicated to the exploration of dark matter (DM) scenarios where the DM particle does not appear in isolation but as part of a richer dark sector, which may in particular feature new strong interactions [1][2][3][4][5][6][7][8]. Such dark sectors provide a variety of mechanisms to set the DM relic abundance [9][10][11][12][13] and lead to novel signatures at collider and fixed-target experiments [14][15][16][17][18][19]. Furthermore, dark sectors with significant self-interactions may explain astrophysical small scale structure observations that are in tension with predictions of collisionless cold DM simulations [20].
The phenomenology of the model depends decisively on the internal properties of the dark sector, in particular the number of dark quark flavours and their respective charge assignments. These properties determine which mesons are stable and which mesons can decay into Standard Model (SM) particles, as well as the processes that can contribute to the freeze-out of the dark sector. Based on these considerations, one obtains rather strong constraints on the structure of the dark sector if cosmological constraints are to be satisfied and the observed DM relic abundance is to be reproduced. In this paper we focus on a QCD-like SU(3) dark sector, in which the dark quarks form various bound states, in particular dark pseudoscalars π and vector mesons ρ. We identify the case of 2 dark quark flavours as particularly interesting since this guarantees that all dark pions JHEP01(2020)162 are stable, avoiding strong constraints from the decay of dark sector particles in the late and early Universe. 1 The dark pion stability allows us to study the phenomenology of strongly interacting dark sectors at ground based experiments while simultaneously being consistent with all cosmological constraints. In our set up we find that the DM relic density is set by the annihilation process ππ → ρρ, which is kinematically forbidden in the present Universe [21].
At the same time, there is substantial freedom in the form of the portal interaction that determines how the dark sector as a whole couples to the SM. The main effect of such interactions is that they induce decays of dark ρ 0 mesons into SM particles. Relatively weak interactions are sufficient for the ρ 0 decays to proceed sufficiently fast to keep the dark sector in thermal equilibrium with the SM during freeze-out, while somewhat stronger interactions allow for the dark sector to be probed with direct detection experiments and at the LHC. Rho mesons that decay promptly on collider scales as well as long-lived mesons are a generic possibility of these models, leading to a wide range of novel collider signatures, such as semi-visible jets [22][23][24][25] and emerging jets [17,26].
In the present work we therefore adopt a hybrid approach, in which the dark sector is modelled in as much detail as possible (given the inherent limitations from non-perturbative physics), while the interactions with the SM are described using a simplified model similar to the ones commonly used in the context of LHC DM searches [27,28]. Specifically, we consider a spin-1 mediator Z with vector couplings to both SM and dark quarks and no other interactions. For this coupling structure direct detection constraints are known to be quite strong if the DM mass is sufficiently above the GeV scale [29,30]. However, the typical mass scale of the dark sector is largely unconstrained and can easily be of the order of a few GeV, such that direct detection constraints are weakened. If furthermore the Z mediator has a mass at the TeV scale, direct detection constraints are additionally suppressed, while LHC constraints become relevant, leading to an interesting interplay between the two different search strategies.
While simplified models of strongly interacting dark sectors have been proposed previously in the literature (see e.g. ref. [31]), our approach differs in that the dark sector is constructed in a way that guarantees a consistent cosmology, independent of the details of the portal interaction. Conversely, our approach differs from most models of Strongly Interacting Massive Particles (SIMPs) with a Z or dark photon mediator in the literature [16,32] in that we allow for a more general coupling structure, in which the interactions of the Z are dominated by its direct couplings to SM fermions (rather than interactions induced by mixing).
We find that large parts of the interesting parameter space are presently unconstrained by direct detection experiments and by LHC searches for di-jet resonances or missing transverse energy. However, the fact that a typical dark shower contains both stable and unstable dark mesons means that existing LHC searches are not optimised for the characteristic signatures of strongly interacting dark sectors. We propose dedicated searches for dark JHEP01(2020)162 showers, which offer great potential to explore the model in more detail, and point out the relevance of searches for displaced vertices [33]. At the same time, significant efforts are being made to develop new direct detection strategies for low-mass DM, which will substantially improve the sensitivity to light dark sectors [34].
This paper is structured as follows. In section 2 we present the model that we consider and derive the corresponding chiral Lagrangian. A particular emphasis is placed on the discussion of meson stability. A broad overview of the phenomenology of the model is then given in section 3. We calculate the lifetime of the unstable particles, the relic density and direct detection constraints. Finally, section 4 takes a closer look at various LHC searches. We first consider existing constraints from searches with visible and invisible final states and then discuss the potential of dedicated searches for the specific signatures of our model. Our conclusions are presented in section 5.

Dark sector model set-up
We consider a dark sector consisting of N f flavours of dark quarks q d = (q d,i ) with i = 1 . . . N f , which are in the fundamental representation of a dark SU(N d ) gauge group. The corresponding Lagrangian is given by where M q denotes the dark quark mass matrix. In order to have a theory that resembles QCD, we assume N d = 3. For reasons that will become clear below, we furthermore restrict ourselves to the case N f = 2. The dark sector described by eq. (2.1) is completely secluded from the SM. Such secluded dark sectors can have a viable cosmology, for example in models with asymmetric reheating [35,36], but make few testable predictions. We therefore focus on the case that there is an additional interaction between the two sectors, which establishes thermal equilibrium in the early Universe and allows for the exchange of entropy. We assume that these interactions arise from an additional U(1) gauge group under which both dark quarks and SM quarks are charged. The U(1) is broken such that the corresponding Z gauge boson acquires a mass m Z . The two dark quarks are taken to have opposite charges under the U(1) such that the interactions with the Z can be written as where e d denotes the product of the U(1) gauge coupling and the charge of the dark quarks and Q = diag(1, −1) is the dark quark charge matrix. The assignment of the U(1) charge Q is of relevance to the stability of dark mesons, as discussed in detail below. On the SM side, we consider a universal coupling of the Z to all SM quarks:

JHEP01(2020)162
Couplings to leptons, as well as mixing between the Z and SM gauge bosons, are assumed to be sufficiently suppressed to be irrelevant for phenomenology. Our set-up hence resembles the simplified model of a vector mediator that is frequently used for the interpretation of DM searches at the LHC [27,28]. At some scale Λ d the dark sector confines and the dark quarks form bounds states in the form of dark mesons and dark baryons. In the present work we will focus on the dark mesons, assuming that the dark baryons are sufficiently heavy that they are cosmologically irrelevant and do not contribute to the present-day DM density. More specifically, we will be interested in the pseudoscalar mesons π and the vector mesons ρ. 2 The former can be understood as the Goldstone bosons associated with chiral symmetry breaking It follows that the number of pions equals the number of broken generators T a , with a = 1 . . . N 2 f −1. We use the normalization Tr(T a T b ) = δ ab /2. For N f = 2, the pion matrix thus reads where π 0 , π + , π − denote the U(1) charge eigenstates, i.e. the π ± have charge ±2e d and the π 0 are uncharged. As we will discuss in more detail below, in our set-up all three pions are stable and hence they contribute equally to the DM relic abundance. Below the confinement scale, the interactions of the pseudoscalar mesons can be described by a chiral effective field theory (EFT), written in terms of where f π denotes the dark pion decay constant. For example, the mass term for the pseudoscalars is given by [16] L mass = f 2 where B 0 is a dimensionless constant. If both dark quark masses are the same, i.e. M q = diag(m q , m q ), one finds with m 2 π = 2B 0 m q . The detailed expressions for the interactions of dark pions with each other as well as the interactions between dark pions and the Z gauge boson are provided in appendix A.
The vector mesons can be parametrised as [16,37] Throughout this work π and ρ always refer to dark mesons, not to their SM counterparts. 3 In general, Vµ also contains the vector ω, and π contains the pseudoscalar η. However, we assume that these particles are sufficiently heavy that they play no important role in the phenomenology of the model and can therefore be omitted in this study. Table 1. Overview of independent parameters in the model.

JHEP01(2020)162
It is furthermore helpful to define the vector meson field strength where g is the pion-vector-meson coupling strength. The so-called KSRF relation [38,39], which relates properties of the ρ meson to the pion decay constant, implies [16] g ≈ m ρ √ 2f π , (2.10) which we use to express f π in terms of m ρ and g. The detailed expressions for the interactions between the pseudoscalar and vector mesons can again be found in appendix A. The interactions between the Z and the charged vector mesons arise from the term which induces in particular the ρ + ρ − Z vertex [37]. For this vertex to have the correct normalisation consistent with the charge ±2e d of ρ ± requires that g Z V = g, i.e. the pionvector-meson coupling strength defined above. Furthermore, the term (2.11) gives rise to mixing between the Z and the ρ 0 . This mixing is of central importance for the phenomenology of our model, as it induces small couplings between the ρ 0 and SM quarks, which render the ρ 0 unstable (see appendix A). A detailed calculation of the resulting ρ 0 lifetime will be provided in section 3.1.
To summarise, in the perturbative regime our model can be completely characterised by the five parameters m q , Λ d , m Z , e d and g q . In the confinement regime, the first two parameters are replaced by the three effective parameters m π , m ρ and g. An overview of all particles in the dark sector and the corresponding parameters is given in table 1.
Before exploring the phenomenology of our model in detail, let us briefly discuss how our model differs from similar scenarios with three flavours, as considered for example in [40]. First of all, for N f = 3, the Chiral EFT Lagrangian includes the Wess-Zumino-Witten (WZW) term [41,42] 2N d 15π 2 f 5 π µνρσ Tr (π∂ µ π∂ ν π∂ ρ π∂ σ π) , (2.12) which induces the 3π → 2π annihilations crucial to the SIMP mechanism [40]. For N f = 2 this term vanishes due to its antisymmetry under pion exchange, so that we need to consider alternative mechanisms for obtaining the DM relic abundance (see section 3.2). Further interaction terms involving π and ρ arise when the WZW term (2.12) is gauged [37].

JHEP01(2020)162
However, for N f = 2 all anomalous π-ρ interactions vanish, because Tr σ a {σ b , σ c } = 0 and therefore the SU(2) is anomaly-free. The most important difference however concerns the pion stability. While the charged dark pions are guaranteed to be stable if there are no lighter states carrying U(1) charge, the π 0 can in principle decay into qq or qqqq. Such decays result for example from the triangle anomaly, which gives rise to a π 0 Z Z vertex in complete analogy to the π 0 SM γγ vertex in the SM. Even if the anomaly vanishes (i.e. if the dark quark charge matrix satisfies Q 2 ∝ 1), neutral pions can still decay via higher-order terms of the WZW type in the chiral Lagrangian [43]. Such decays turn out to be extremely dangerous for the viability of the model. If π 0 decays are fast, they will keep the pions in equilibrium with the SM and lead to an exponential suppression of the DM relic abundance. Slow decays, on the other hand, are subject to strong constraints from nucleosynthesis and recombination.
To evade these constraints, one can impose a dark G-parity [16] G where C is the U(1) charge conjugation operator, and we have introduced a Z 2 transformation that takes Z → −Z , as well as an SU(N f ) transformation U q that satisfies In combination with the requirement Q 2 ∝ 1, eq. (2.14) can only be satisfied if the number N f of dark quarks is even. Specifically, for N f = 2 and Q = diag(1, −1) one finds such that π 0 is indeed odd under G d . For N f = 3, on the other hand, it is not possible to define an analogous G-parity and hence π 0 decays cannot be forbidden.

General phenomenology
For the discussion above there was no need to specify the mass scale of the dark sector or the mass of the Z mediator. From now on we will be more specific and consider a GeV-scale dark sector that interacts with the SM via a TeV-scale Z . As discussed in detail below, GeV-scale dark sectors provide various interesting mechanisms to produce the thermal relic density, and are less constrained by direct detection experiments than heavier DM particles. The phenomenology of dark sectors with light Z mediators has been explored in the literature, see e.g. [16]. Here we focus on heavy TeV-scale mediators which lead to new avenues for dark sector searches at the LHC. Since the Z is heavy compared to the confinement scale of the dark sector, we can calculate its production and decay in terms of free (SM and dark) quarks -this will be the topic of section 4. At low energies and in the present Universe, however, the appropriate degrees of freedom are the pseudoscalar and vector mesons that appear in the confinement phase. For the reasons outlined above, all dark pions are exactly stable in the model we JHEP01(2020)162 consider and can therefore potentially account for the DM in the Universe. In the present section, we will study the mechanisms that determine the relic abundance of dark pions, as well as constraints from direct detection experiments. Before doing so, we however need to consider the properties of the rho mesons, which turn out to be of central importance for the phenomenology of the model.

ρ 0 lifetime
In the presence of the interaction term in eq. (2.11) the ρ 0 meson mixes with the Z boson. This mixing gives rise to interactions between the ρ 0 and SM quarks, which to leading order in m ρ /m Z can be written as For large Z masses the induced couplings can be extremely small and hence the ρ 0 can potentially be rather long-lived. Calculating the ρ 0 lifetime for masses in the GeV range is a notoriously difficult problem. For m ρ 2 GeV and away from any spin-1 QCD resonances, one can employ the perturbative spectator model to estimate and the sum includes all quarks that are kinematically allowed, i.e. with m q SM < m ρ /2. For m ρ < 2 GeV, decays of the ρ 0 are more accurately described by calculating the mixing with QCD resonances, in particular with the SM ρ meson. However, since the Z in our model, and hence the ρ 0 , couples to all quarks with equal strength, tree-level interactions with SM mesons are absent. Decays into mesons can therefore only proceed via baryon loops, so that we expect the decay width of the ρ 0 to become very small in this mass range. In the present work, we will not consider this problem further and focus on m ρ 2 GeV. Finally, we assume that, contrary to QCD, m ρ < 2 m π such that decays into the dark sector are not possible kinematically. Such a small mass difference can arise if the masses of the dark quarks are comparable to the confinement scale, such that the explicit breaking of chiral symmetry is stronger than for QCD.
For concreteness let us introduce a benchmark scenario that we will study in detail in the following. We set e d = 0.4, g = 1, m ρ = 5 GeV and m Z = 1 TeV. We then find which corresponds to As we shall discuss in section 4, small couplings g q lead to displaced ρ 0 decays and thus interesting LHC phenomenology. The ρ ± can in principle decay into π ± qq via an off-shell Z coupled to π and ρ through an anomalous vertex of the gauged WZW type. However, for m Z m ρ , we find this channel to be extremely suppressed by the three-body nature of the decay as well as powers of momentum in the vertex factor. In addition, if m ρ − m π 2 GeV, the appropriate final states would again be SM mesons rather than quarks. As discussed above, however, tree-level couplings of the Z to SM mesons are absent. We can therefore treat ρ ± as effectively stable during dark sector freeze-out. As we will see below, the abundance of ρ ± is exponentially suppressed relative to the dark pion abundance, so that these particles do not contribute to the DM and their late-time decays do not lead to any observable effects.
Apart from the interesting phenomenological implications to be discussed below, the interactions between the ρ 0 and SM quarks play an important role in the cosmological evolution. Indeed, the (inverse) decays of the ρ 0 are the primary way in which the dark sector can maintain thermal equilibrium with the SM bath. For these decays to be efficient, one requires that where H ∼ 14.4 T 2 /M P is the Hubble rate before the QCD phase transition. For our benchmark point this corresponds to If this condition is satisfied, the temperature of the ρ 0 is equal to the SM temperature and the number density n ρ 0 is given by an equilibrium distribution, n ρ 0 = n eq ρ . The strong interactions between the ρ 0 and the ρ ± will then ensure that the same also holds for the charged rho mesons. Initially, interactions between the dark pions and rhos will also keep the dark pions in equilibrium. DM freeze-out happens when these interactions become inefficient and the dark pions decouple from the rho mesons, which will be discussed next.

Relic density from forbidden annihilations
The main process that keeps dark pions in thermal equilibrium with the dark rhos (and hence with the SM) is the pair annihilation ππ → ρρ (see figure 1). This process is JHEP01(2020)162 fully efficient as long as the temperature is large compared to m ρ − m π , but becomes exponentially suppressed for smaller temperatures [16,21,44]. Provided that the dark ρ mesons are always in thermal equilibrium with the SM, the Boltzmann equation for the dark pion number density n π readṡ Since the right-hand side has to vanish if both the dark pions and the dark rhos are in equilibrium, the cross section for forbidden annihilations can be expressed as (3.8) Here, σ ρρ→ππ v is unsuppressed at low temperatures and scales proportional to g 4 /m 2 π . The exponential suppression arises from the factor with x = m π /T and ∆ = (m ρ − m π )/m π . Inserting eq. (3.8), the Boltzmann equation (3.7) takes the formṅ From (3.10) we can read off that ππ → ρρ interactions decouple when For example, for the benchmark point above and m π = 4 GeV, such that ∆ = 0.25, one finds that decoupling happens for x ≈ 26. The resulting abundance of dark pions is found to be close to the observed value Ω DM h 2 = 0.12 [45]. Due to the exponential sensitivity of Ω DM to the mass splitting ∆, it is always possible to reproduce the observed value by varying m π (or alternatively m ρ ) slightly. We compute the relic density using MicrOMEGAs 5.0.6 [46], which allows for up to two dark species in its freeze-out calculation. Conversions between π ± and π 0 remain efficient until long after annihilations of dark pions into dark rho mesons have frozen out. Therefore, we can treat all dark pions as one species in MicrOMEGAs. The ρ ± are assigned as a second dark species, whose abundance is strongly suppressed at low temperatures. 4 Furthermore, in the freeze-out scenario considered here, ρ 0 is in equilibrium with the SM throughout dark-pion freeze-out, and can hence be treated like an SM particle by MicrOMEGAs.

JHEP01(2020)162
Let us briefly revisit our assumption that the dark rho mesons are in thermal equilibrium with the SM during dark pion freeze-out. In order for the dark rhos to remain in equilibrium, (inverse) decays should be at least as efficient as the conversion between dark rhos and pions until the point when the dark pions freeze out. This requires that before dark pion freeze-out, and hence Hence, requiring that rho decays not be a limiting factor for pion freeze-out yields the lower bound (3.13) on Γ ρ 0 , which is more stringent than the simple requirement of thermalisation in eq. (3.5).

Constraints from direct detection experiments
At low energies dark rho exchange induces an effective coupling of π ± to SM nucleons N = p, n, given by which turns out to depend on m Z rather than m ρ because of the way in which the interactions arise from ρ-Z mixing (see appendix A). This effective interaction gives rise to spin-independent scattering with cross section given by where µ πN = m π m N /(m π + m N ) is the reduced mass.
Since we have assumed universal quark couplings, the DM-nucleus cross section receives a coherent enhancement proportional to the square of the mass number A. We can therefore directly compare our model predictions for σ SI N to published exclusion limits and obtain a bound on the effective coupling e d g q /m 2 Z as a function of m π . However, it is important to account for the fact that neutral pions do not couple to SM quarks at tree-level. Thus, for the purpose of direct detection the effective local DM density is reduced by a factor 2/3, which can be captured by an appropriate rescaling of published exclusion limits.
For the mass range 1 GeV m π 10 GeV that we will be interested in, relevant constraints arise from a number of different direct detection experiments, namely CRESST-III [47], CDMSLite [48], DarkSide-50 [49], PICO-60 [50], PandaX [51] and XENON1T [52]. Rather than simply considering each experimental result separately, we use DDCalc 2.0 [53] to perform a statistical combination of all experimental results. However, as explained in detail in appendix B, we do not include the DarkSide-50 analysis, which relies on an overly optimistic extrapolation of the ionisation yield to low energies. In addition, we separately consider sensitivity projections for the SuperCDMS experiment [54], which should provide substantial improvements for small DM masses. The resulting constraints on the parameter space are shown in figure 2. While the direct detection constraints depend only on the parameters that are being varied explicitly, we can include additional information in the figure by fixing the π-ρ coupling g. Doing so enables us to calculate m ρ as a function of m π from the relic density requirement, as indicated by the second x-axis at the top. For given m ρ we can then determine the ρ 0 decay length and indicate the parameter region where one can expect LHC signatures with displaced vertices. Furthermore, we show the parameter region where the requirement of thermal equilibrium during dark pion freeze-out is violated (see eq. (3.5)), as well as the parameter region where ρ 0 decays are inefficient (see eq. (3.13)).

Astrophysical constraints
Let us finally consider potential constraints from astrophysical observations. Indirect detection experiments typically give strong constraints on models of thermal DM with a particle mass below 10 GeV. In our model, however, the dominant annihilation channel for dark pions is ππ → ρρ. Since m ρ > m π , this channel is only kinematically open in the early Universe and becomes strongly suppressed for small relative velocities. Thus, there are no relevant indirect detection constraints for our model.
One of the attractive features of strongly interacting dark sectors is their potential to resolve the so-called small-scale crisis of ΛCDM with DM self-interactions. In our model, the self-scattering cross section for π + π − → π + π − divided by the dark pion mass is given by .

(3.16)
Similar results hold for the scattering channels involving neutral pions.

JHEP01(2020)162
Hence, for GeV-scale dark pion masses the self-scattering cross section is too small to induce observable effects in small astrophysical systems, which would require σ self /m π 1 cm 2 /g. For the same reason, the parameter space that we consider in this work easily satisfies all constraints on the self-interaction cross section from merging galaxy clusters.

LHC constraints
The LHC phenomenology of our model is dominated by the on-shell production of a Z boson and its subsequent decay. Since the mass of the Z is assumed to be large compared to the confinement scale of the dark sector, the partial decay widths can be calculated in terms of free quarks: where the sum runs over all SM quarks and m q denotes the corresponding quark mass. In the absence of any other decay modes, the branching ratio into dark sector states is then given by It is worth noting that Γ(Z → q d q d ) is proportional to N f × N d and hence the branching ratio into dark sector states can be sizeable without the need for a large hierarchy of couplings between e d and g q . For example, for e d = 0.4 and g q = 0.1, the branching ratio into dark sector states is found to be 84%. The case of the Z decaying into SM quarks is constrained by LHC searches for di-jet resonances. A variety of search strategies have been developed to search for such resonances in different mass ranges. The strongest bound for high-mass resonances stems from a recent ATLAS search based on an integrated luminosity of L = 139 fb −1 [55]. A compilation of different constraints on low-mass resonances was recently presented in ref. [56]. We apply these constraints to our model by rescaling them with the appropriate branching ratios, which is a good approximation because the total width of the Z is sufficiently small.
If the Z decays into dark quarks, fragmentation and hadronisation proceeds within the dark sector (see figure 3). We simulate these dark showers using the Hidden Valley model in Pythia [57,58], which calculates the final yield and distribution of dark mesons. The number of dark mesons produced within each shower depends on the initial energy of the dark quark and varies from event to event. On average, 10-12 dark mesons are produced from a dark quark with an energy of 1 TeV, but events with more than 20 mesons per shower are not uncommon.
Accordingly, we find the typical boost factor for dark mesons to be of the order of γ ≡ E ρ /m ρ ≈ 10, where E ρ denotes the average energy of ρ mesons. Because of this boost factor, the actual decay length of the ρ 0 mesons in the laboratory frame, γcτ , is substantially larger than the estimate given in eq. (3.4). In the following, we will assume that for γcτ < 1 mm the rho meson decays can be treated as prompt, such that conventional experimental strategies apply. However, since both the boost factor and the actual distance travelled before the decay are subject to large fluctuations, displaced vertices may be observable even for smaller decay lengths. The average relative multiplicity of the different mesons depends on their respective number of degrees of freedom. Spin-1 ρ mesons are three times as abundant as spin-0 π mesons and charged ρ ± and π ± mesons are twice as abundant as their neutral partners. It follows that we expect on average 25% of a dark shower to consist of ρ 0 mesons, which subsequently decay into SM hadrons, while the remaining 75% are stable mesons that escape from the detector unseen. A dark shower will hence lead to a semi-visible jet [22,23] with an average fraction of invisible energy of r inv = 0.75.
Such semi-visible jets give rise to a range of interesting experimental signatures. If the Z is produced in isolation, i.e. without additional energetic SM particles from initial state radiation (ISR), the two semi-visible jets will be back-to-back. Defining the minimum angular separation in the azimuthal plane between the missing energy vector / E T and any of the leading jets ∆φ = min such events are expected to have small ∆φ, as the missing energy is aligned with one of the dark showers. Ordinary "mono-jet" searches (i.e. searches for energetic jets in association with missing energy) will reject such events because of prohibitive QCD backgrounds from misreconstructed jets [59,60]. Traditional searches for di-jet resonances are also expected to be insensitive to these kinds of events, since the visible jets only carry a fraction of the energy of the dark shower and hence their invariant mass does not peak at the mass of the Z . However, given the relatively large value of r inv , there is a non-negligible probability for a dark shower to remain entirely invisible. In this case, the Z decay would lead to a JHEP01(2020)162 missing energy vector pointing away from the visible jets, such that ∆φ ≈ π and the event selection cuts applied in mono-jet searches can be satisfied. While the missing energy that can be obtained in this way is limited to / E T < m Z /2, larger amounts of missing energy are possible if the Z recoils against an ISR jet (see figure 3). In this case, also events with two partially visible jets can contribute, as the missing energy vector will in general not be aligned with any of the jets.
These considerations are illustrated in figure 4, which shows the double differential cross section with respect to ∆φ and N j after imposing the requirement / E T > 250 GeV. In the left panel, we consider the case of dark showers with r inv = 0.75, whereas the right panel corresponds to fully invisible Z decays, such that jets can only arise from ISR. In the first case, one can clearly see the two contributions discussed above, where either N j = 1 and ∆φ ≈ π or N j > 1 and ∆φ ≈ 0. For the case of fully invisible decays, on the other hand, the second contribution is absent.
In the following, we will first consider existing constraints on the parameter space from searches for di-jet resonances and searches for jets in association with missing transverse energy. Afterwards, we discuss a novel search strategy that targets semi-visible jets more specifically and estimate the corresponding sensitivity. As we will see, all these searches probe the parameter regions where the ρ 0 is short-lived and decays promptly on collider scales. Nevertheless, even for smaller couplings the Z production cross section can be sizeable, leading to events with displaced vertices, which would look similar to the emerging jet signature discussed in ref. [26]. A detailed analysis of such signatures will be left to future work.

Constraints from missing energy searches
Mono-jet signatures can arise in our model if the dark showers remain mostly invisible and recoil against an ISR jet, or if some part of the dark showers produces an energetic jet while the rest remains invisible. Moreover, since mono-jet searches typically allow for more than JHEP01(2020)162 one jet, events with multiple jets from ISR or the dark shower can also contribute. On the other hand, events where both dark showers are partially visible can have similar kinematics as squark or gluino pair production. In both cases, events consist of two hemispheres with a jet and missing energy in each hemisphere. Hence, supersymmetry searches for jets and missing energy can be sensitive to dark shower production.
Simulation details. Having implemented our model with the FeynRules package [61], we generate parton-level events at leading order for the dark quark production process pp → q d q d with MadGraph5 aMC@NLO [62] using the NN23LO1 PDF set [63]. We perform MLM matching of samples with up to one hard jet setting xqcut = 100 GeV. The Z width is calculated self-consistently by MadGraph at each parameter point. The parton-level events are then showered and hadronised, both in QCD and in the dark sector, with Pythia 8 [64].
The simulation of the dark shower and of dark meson production relies on Pythia's Hidden Valley module [57]. The hidden valley meson states pivUp, pivDn, pivDiag, rhovUp, rhovDn and rhovDiag provided by Pythia8 are equivalent to the dark mesons π + , π − , π 0 , ρ + , ρ − and ρ 0 , respectively. Out of these, only ρ 0 (rhoDiag) decays into SM quarks while the others are stable. We adjust the Hidden Valley variable probVector appropriately to implement the expected ratio of invisible mesons r inv = 0.75. Furthermore, we turn on the running of the dark coupling α d determined by the confinement scale Λ d .
We scan over a range of Z masses between 500 GeV and 5000 GeV, generating 10 5 events for each Z mass. The other relevant parameters are set to the benchmark values m q = 500 MeV, m π = 4 GeV, and m ρ = Λ d = 5 GeV. Larger meson masses would reduce the average meson multiplicity, while different values of m q and Λ d would not significantly change our results.
Recasting and analysis details. We recast existing experimental analyses with Check-MATE 2 [65,66], and MadAnalysis 5 [67,68] in conjunction with its Public Analysis Database (PAD) [69]. Both these codes first pass the hadron-level events from Pythia8 to Delphes3 [70], which simulates the appropriate detector for the respective analysis.
Delphes3 internally calls the FastJet [71,72] package for jet clustering. Subsequently, the analysis cuts are applied to the signal events and 95%-exclusion limits S 95 on the number of signal events are calculated based on the number of observed background events. To incorporate the signal uncertainty without the time-consuming computation of CLs values [73], CheckMATE calculates the ratio r defined as  four jets with p T > 30 GeV. For the angular distance between the missing energy and the leading jets the search requires that ∆φ > 0.4, with ∆φ as in eq. (4.4). A range of inclusive and exclusive signal regions are defined in terms of the amount of missing transverse energy. Among SUSY searches available for recasting we find the most recent CMS squark and gluino search [76] to have the highest sensitivity to dark shower signal in the Z mass range we consider. This search uses 35.9 fb −1 of data at √ s = 13 TeV and is implemented in the MadAnalysis5 PAD [77]. In contrast to the mono-jet analysis, here at least 2 jets are required. Events need to fulfil / E T > 300 GeV and H T > 300 GeV, where H T denotes the scalar sum of transverse momenta of all jets. Moreover, ∆φ(j, / E T ) > 0.5 is imposed on the two leading jets, and ∆φ(j, / E T ) > 0.3 on the third and fourth jet. Each signal region is defined by the combination of an / E T range and an H T range. For each Z mass, we conservatively approximate the exclusion bound from a given search by the limit from the most sensitive signal region, based on observed event numbers. Since Γ Z /m Z , as given in eqs. (4.1) and (4.2), is well below 10% in the relevant parameter space, we translate bounds on the number of signal events S into bounds on couplings and masses using the narrow width approximation, i.e.
with the branching ratio given in (4.3). We present our results in figure 5, which shows the parameter regions excluded by searches for di-jet resonances and searches for missing energy in the g q -m Z parameter plane for different values of e d . Larger values of e d suppress the branching ratio of the Z into SM quarks and therefore the impact of di-jet constraints while enhancing the sensitivity of searches for missing energy. For m Z 2 TeV the combined constraints imply g q 0.1, while for m Z > 4 TeV values of g q as large as 0.3 are allowed.
We also indicate the combination of parameters corresponding to γcτ ρ 0 = 1 mm for γ ≈ 10. As discussed above, ρ 0 decays can be treated as prompt above this line, so that JHEP01(2020)162 the analyses discussed above can be safely applied. The search for di-jet resonances, on the other hand, does not involve any dark mesons and is therefore not affected by the ρ 0 lifetime. Below the green dashed line, one would generically expect displaced vertices. We emphasise that the number of such events can potentially be quite large. For concreteness, we show the combinations of g q and m Z that correspond to σ pp→Z = 10 fb, which would correspond to more than 10 3 such events having been produced to date at ATLAS and CMS.
To conclude this discussion, we emphasise that the analyses considered above have not been optimised to search for dark showers. The mono-jet search, for example, is intended to search for Z bosons that decay fully invisibly. In this case, signal events are expected to have large ∆φ (see figure 4), so that a cut on ∆φ substantially improves the signal-tobackground ratio. In our case the situation is quite different, as events with N j > 1 and small ∆φ contribute substantially to the signal cross section. This makes it important to reassess whether the requirement ∆φ > 0.4 is strictly necessary to suppress backgrounds. Conversely, it may be interesting to reject events with N j > 1 and ∆φ > 2, which are more likely to arise from background than from our model.

Prospective searches for semi-visible jets
Having reinterpreted existing generic searches for jets and missing energy, let us now discuss dedicated searches for the dark showers produced in the decay of the Z boson. If the dark showers were completely visible, one would obtain a peak in the total invariant mass of the resulting jets, which would be centred at the Z mass and could be used to distinguish signal from background. However, for semi-visible jets the peak in M 2 jj = (p j 1 + p j 2 ) 2 is washed out to a considerable extent by the fact that the fraction of visible energy in the jets differs from event to event.
We therefore consider an analysis proposed in refs. [22,23], which is based on the transverse jet mass where p T jj is the vector sum of the p T of the two jets. It was shown in ref. [23] that this variable maintains the ability to distinguish signal from background up to rather large values of r inv and can be competitive with ordinary missing energy searches for r inv ≈ 0.75.
To apply this analysis to our model we pass the events generated by the tool chain described in section 4.1 to Delphes3 with CMS settings and use the anti-k T algorithm with R = 0.5 for jet clustering. We then apply the same cuts as in ref. [23]. At the preselection level, these require / E T > 200 GeV and at least two jets with p T > 100 GeV and |η| < 2.4. Subsequently, jets are re-clustered with the Cambridge-Aachen algorithm and R = 1.1. We then compute the transverse mass of the two leading re-clustered jets and require / E T /M T > 0.15. Moreover, we impose |η j 1 −η j 2 | < 1.1. Electrons with p T > 10 GeV and |η| < 2.4 as well as muons with p T > 20 GeV and |η| < 2.4 are vetoed. Finally, the analysis imposes an angular separation cut that is inverted compared to the mono-jet case and requires ∆φ < 0.4 between the missing energy vector and at least one jet.
We show in figure 6 the M T spectrum expected for a Z boson with m Z = 2 TeV decaying into two dark showers with r inv = 0.75. We have furthermore set e d = 0.6 and JHEP01(2020)162 Figure 6. Differential cross section with respect to the transverse jet mass M T for semi-visible jets obtained from the dark showers produced by the decays of a Z boson with m Z = 2 TeV and for the SM backgrounds simulated in ref. [23]. g q = 0.1, compatible with current di-jet constraints (see figure 5). For comparison, we also show the background estimate obtained in ref. [23] rescaled to L = 300 fb −1 . One can clearly see the difference in shape between signal and background arising from the fact that M T ≤ m Z for the decay of an on-shell Z .
We estimate the sensitivity of this search by calculating the likelihood ratio L of signal+background and background-only and requiring −2 log L > 3.84 in order for a given parameter point to be testable. The resulting sensitivity estimate is shown in figure 7 in comparison to the combined bounds from existing searches (see figure 5). We find that the M T search has the power to probe larger Z masses than the missing energy searches discussed above and can be competitive with searches for resonances in visible decays for sufficiently large e d . 5 The M T search that we have implemented is not optimised for m Z < 1 TeV and hence we do not show sensitivity projections in this mass range. Searches for low-mass resonances decaying into jets are notoriously difficult to observe because of large QCD backgrounds. While the missing energy requirement is expected to improve the situation, reliable estimates of backgrounds with fake missing energy due to misreconstructed jets are difficult to obtain. To suppress these backgrounds it may be necessary to require additional particles in the event, for example a high-p T photon from ISR (in analogy to the recent ATLAS search presented in ref. [78]). We leave a detailed study of this mass range to future work.
We finally note that, following ref. [23], the sensitivity estimate shown in figure 7 assumes no systematic uncertainties. However, as shown in figure 6, backgrounds are 5 The projected sensitivity of the MT analysis extends into the parameter region where not all ρ 0 decays can be treated as prompt. Jets originating from displaced vertices would likely be rejected in such an analysis, which would reduce the sensitivity slightly compared to what is shown in figure 7. large and the signal is rather broad. Hence, even small uncertainties in the shape of the background may render a potential signal unobservable. We find that, in order to achieve a sensitivity close to the one shown in figure 7, systematic uncertainties must be smaller than 1%. For the more realistic assumption of 5% systematic uncertainties, the sensitivity of the M T offers no improvement over the searches discussed in section 4.1. 6 Detailed studies of systematic uncertainties will therefore be essential to obtain realistic sensitivity estimates.

Conclusions
In this work we have investigated how the phenomenology of strongly interacting dark sectors depends on their internal structure and on their interactions with the Standard Model. For simplicity, we have considered portal interactions arising from a Z mediator with vector couplings to both dark and visible quarks, but many of our results are insensitive to the details of this interaction. Below the confinement scale the dark quarks form π and ρ mesons, such that the dark sector can be characterised in terms of the masses m π and m ρ and the coupling strength g of their interactions.
The case of two dark quarks with opposite U(1) charges turns out to be particularly interesting as all dark pions are then stable, allowing for a non-zero DM relic density as well as avoiding constraints from late decays of DM particles. The ρ 0 on the other hand has a sizeable decay width into SM particles, which is sufficient to establish thermal equilibrium between the dark sector and the visible sector in the early Universe. Dark sector freeze-out then proceeds via the forbidden annihilation process ππ → ρρ. The observed relic abundance can hence be reproduced independently of the strength of the JHEP01(2020)162 portal interactions, provided m π and m ρ are sufficiently close. For example, for m π = 4 GeV and g = 1 we require m ρ ≈ 5 GeV. The observation that the cosmological properties of the dark sector and the relic abundance of dark mesons are largely independent of the interactions between the dark sector and the Standard Model justifies our simplified description. Additionally, the details of our dark sector have been chosen in a way to allow for a fully consistent cosmological picture to emerge, allowing us to consider LHC phenomenology from a dark sector that can realistically account for the observed DM in our universe.
Experimental predictions, on the other hand, depend more sensitively on the assumed portal interactions. For the Z mediator we consider, charged dark pions interact with SM nuclei via vector boson exchange, which gives rise to spin-independent scattering in direct detection experiments. By combining results from a number of different direct detection experiments, we obtain strong constraints on the coupling strength of the Z . Based on these constraints we have identified dark pion masses of the order of a few GeV as particularly interesting. Even for dark pions in this mass range, direct detection constraints require the Z mediator to be either heavy or weakly coupled. In the present work, we have explored the former possibility, which leads to interesting implications for LHC physics.
LHC phenomenology for this model is dominated by the on-shell production of the mediator (possibly in association with SM particles) and its subsequent decays into either visible or dark quarks. While the former case leads to di-jet resonances, that can be easily reconstructed, the latter case gives rise to dark showers. We emphasise that for a strongly-interacting dark sector the branching ratio into dark quarks is enhanced by colour and flavour factors, making the latter signature particularly important. Since most of the mesons in the dark shower are stable, one obtains so-called semi-visible jets, in which only a small fraction of the initial energy of the dark quarks can be detected.
We find that various existing LHC searches for missing energy can be sensitive to such a scenario. Nevertheless, these searches are not optimised for the case of dark showers, where the missing energy tends to be aligned with one or more visible jets, which is difficult to disentangle from QCD backgrounds from mis-reconstructed jets. We therefore consider an alternative approach, in which the information from all visible final states is combined to calculate the transverse mass M T to achieve better discrimination between signal and background. The expected sensitivity of such a search depends however sensitively on the assumed systematic uncertainties for the background distribution. We find that systematic uncertainties as low as 1% will be necessary in order to improve upon existing searches.
At the same time we identify a number of exciting directions for future research. First of all, a generic prediction of the model that we study is that the ρ 0 can be long-lived, which would give rise to displaced vertices at the LHC. The corresponding production cross sections can be quite large, and it is conceivable that thousands of such events have already gone unnoticed at ATLAS and CMS. Ongoing detector upgrades as well as new analysis strategies make these signatures a promising target for future LHC runs. But even for prompt decays there is room for substantial improvements. The jets obtained from our model differ quite substantially from ordinary QCD jets, for example because heavy quarks JHEP01(2020)162 are absent in the shower. It will therefore be very interesting to study whether a neural network trained to identify dark showers can help to suppress QCD backgrounds.

Acknowledgments
We thank Alexander Mück, Pedro Schwaller and Susanne Westhoff for discussions. This work is funded by the Deutsche Forschungsgemeinschaft (DFG) through the Collaborative Research Center TRR 257 "Particle Physics Phenomenology after the Higgs Discovery", the Emmy Noether Grant No. KA 4662/1-1 and the Research Unit FOR 2239 "New Physics at the Large Hadron Collider".

A Full Lagrangian
The Lagrangian of the underlying dark SU(3) gauge theory coupled to a Z vector mediator reads where q d = (q d,1 , q d,2 ) and M q = diag(m q , m q ). The dark quark covariant derivative has the form where A µ denotes the dark gluon field and Q = diag(1, −1) is the dark quark charge matrix. The chiral EFT Lagrangian (up to fourth order in the pion fields) is given by Here π denotes the pion matrix with the SU(2) generators T a = σ a 2 . In a similar fashion, we have introduced The pion covariant derivative is given by Eq. (A.1) gives rise to kinetic mixing between the Z and the ρ 0 . To recover canonical kinetic terms, while at the same time keeping the mass term diagonal, the interaction eigenstatesZ andρ 0 can be expressed in terms of the physical fields Z µ and ρ 0 µ as where = arcsin(2 e d /g) and we have only kept terms up to second order in m ρ /m Z . Note that the diagonalisation can only be performed if the matrix of kinetic terms is positive definite, which requires e d /g < 1/2. Because of the mixing, the ρ 0 obtains couplings to SM quarks, which can be written as At the same time, the mixing modifies the interactions between the Z , the ρ 0 and dark pions. For example, the trilinear interactions become Interestingly, the mixing does not modify the interactions between the ρ 0 and π ± but strongly suppresses the interactions between the Z and dark pions. Thus, for m ρ m Z the ρ 0 replaces the Z as the primary mediator between the dark and the visible sector at low energies.

B Sensitivity of DarkSide-50 to low-mass dark matter
In this appendix we take a closer look at the sensitivity of the DarkSide-50 experiment to low DM masses. Ref. [49] claims to provide the world-leading bound on the spinindependent DM-nucleon scattering cross section of DM particles with masses between 1.8 GeV and 5 GeV. To derive this bound ref. [49] adopts the Bezrukov model [79] for the ionisation yield Q y in liquid noble gases, which uses Lindhard theory [80] to predict the electronic stopping power of liquid argon and a model of recombination proposed by Thomas and Imel [81]. This model is then fitted to calibration data to obtain the detector response for low recoil energies. However, as discussed in ref. [79] the use of Lindhard theory for very low energies is highly doubtful and deviations from the simple model are JHEP01(2020)162 expected. Ref. [79] proposes to vary the functional form of the electronic stopping power by multiplying it with a correction factor F (v/v 0 ). It has been shown in ref. [82] that such variations can substantially affect the sensitivity of direct detection experiments to light DM particles. Indeed, we will demonstrate that the sensitivity of DarkSide-50 relies heavily on the assumed extrapolation for Q y .
We implement the DarkSide-50 analysis as follows. The expected number of electrons produced at the interaction point is calculated via N e = Q y (E R ) E R . We assume the fluctuations in N e to follow a Poisson distribution convoluted with a Gaussian distribution that accounts for the detector resolution. The width of the Gaussian distribution is determined as σ Ne = 0.33 by fitting the detector response for N e = 2. 7 An overall detector acceptance of 0.43 is already accounted in the total exposure of 6786 kg days.
We then calculate the predicted number of events in three different N e -bins: [4,7], [7,10], [10,22]. In the first two bins, one finds a clear discrepancy between the background model (predicting ∼ 390 and ∼ 550 events, respectively) and observation (∼ 680 and ∼ 630 events), which suggests an unknown background component. In these bins we therefore do not perform any background subtraction and only include them in the total likelihood if the predicted number of DM events exceeds the number of observed events. In the third bin, on the other hand, background prediction (∼ 4080 events) and observation (∼ 4130 events) are in good agreement, so that one can use the standard Poisson likelihood to perform a background subtraction.
In spite of the rather crude implementation our analysis approximately reproduces the exclusion limit published by DarkSide-50 when making the same assumption on the ionisation yield Q y . We can therefore use our implementation to study the impact of varying Q y on the exclusion limit. The impact of such variations are shown in figure 8. The solid black line corresponds to the official DarkSide-50 bound, the various blue lines correspond to the bounds obtained when adopting different functional forms for the ionisation yield Q y as illustrated in the inset. We emphasise that all the variations that we consider differ only for E R 5 keV, i.e. in an energy range for which no direct measurements of the ionisation yield are available. 8 Figure 8 clearly illustrates the strong dependence of the claimed exclusion limit on the assumed functional form of the ionisation yield used for the extrapolation to lower recoil energies. Indeed, the bound from DarkSide-50 for DM masses of a few GeV can be suppressed by more than an order of magnitude, such that the strongest constraints in this mass range arise from CDMSLite [48] and PICO-60 [83] (also shown in figure 8). We conclude that the DarkSide-50 exclusion limit is not robust unless theoretical calculations or experimental measurements can provide further evidence in support of the assumed functional form of Q y . We therefore do not include the exclusion limit from DarkSide-50 in our analysis. 7 This approach neglects a possible dependence of σN e on Ne, which cannot be inferred from the available information. 8 The ionisation yield in this energy range may still be constrained using calibration data from DarkSide-50. However, to obtain robust constraints at very low recoil energies, the fit to calibration data should include variations in the functional form of Qy like the ones considered here. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.