The Elusive Muonic WIMP

The Weakly Interacting Massive Particle (WIMP) paradigm is one of the most popular scenarios for Dark Matter (DM) theories that however is strongly constrained, in particular by direct detection experiments. We stick with the WIMP hypothesis and consider a Dirac fermion candidate for DM that interacts with the Standard Model (SM) via a spin-1 $Z'$, arising from the spontaneous breaking of an Abelian $U(1)'_{\mu}$ gauge symmetry, under which only second generation leptons and the DM are appropriately charged. Due to the charge assignment, the model is gauge anomalous and can only be interpreted as an effective field theory (EFT) at low energy. The $Z'$ couples at tree level only to the vector DM current, to the axial muon current and to left-handed muonic neutrinos, so the WIMP-nucleon cross section is beyond the experimental reach of spin-independent (SI) direct detection searches. We study the current bounds on this model coming from direct and indirect detection of DM, collider searches, contributions to $(g-2)_{\mu}$ and to neutrino trident production. We find that large regions of the parameter space remains to be explored. In the context of LHC searches, we study the impact of a muon-exclusive signal region for the $3\mu$ + ${E}^{{\rm miss}}_T$ channel with an invariant mass window around $m_{Z'}$. We show that this search can significantly improve the current collider bounds. Finally, from the anomalous nature of our EFT, there remain at low energy triboson anomalous interactions between the $Z'$ and the electroweak (EW) SM gauge bosons. We explore the possibilities of probing these interactions at the LHC and at a 100 TeV proton collider finding it extremely challenging. On the other hand, for a muon collider the resonant channel $\mu^{+}\mu^{-}\to Z'\to ZZ$ could be discovered in the most promising scenario with luminosity of $\mathcal{O}({\rm few}\; 10)$ ${\rm fb}^{-1}$.


Introduction
There exists an overwhelming gravitational evidence for the existence of a new type of matter which is stable on cosmological scales and neutral under Electromagnetism. Its origin and composition cannot be explained by the Standard Model (SM) of particle physics, so it represents one of the greatest puzzles of modern physics. Many theoretical ideas beyond the SM attempt to explain the nature of this dark matter (DM), and within the attempts to describe DM as a particle, one of the most popular ones is that of a weakly interacting massive particle (WIMP) that introduces a new particle associated with the DM with a mass and interactions to SM states that are of order as those expected for an electroweak (EW) state. In that case one naturally obtains annihilation cross sections for the DM that provide via the mechanism known as freeze-out a relic density Ω DM h 2 ∼ 0.1, in accordance with measurements done by several astrophysical probes.
Despite the great appealing of the WIMP paradigm, the physics community would have naturally expected its existence to be detected already by now. There are several searches that probe the non-gravitational interactions of the DM with the SM states. Among them, direct detection experiments that involve the scattering of DM particles against heavy nuclei impose some of the most restrictive constraints [1][2][3], in particular in the parameter space regions in which WIMPS would tend to naturally lie. In fact typical cross sections expected for WIMP scattering against nuclei are already ruled out by several orders of magnitude. As direct detection searches continue probing smaller scattering cross sections and no sign of DM is detected in them, the WIMP paradigm becomes further in tension with data although it still remains as a viable possibility [4][5][6][7]. Another probe of the non-gravitational interactions of the DM are what are known as indirect detection experiments that measure the signals of DM annihilation in high density regions of the Universe, such as the center of the Milky Way [8] or in DM dominated galaxies such as dwarf spheroidals (dSph) [9,10]. Finally it may also be possible to directly produce the DM in particle accelerators, which would then be indirectly detected as missing energy [11][12][13][14][15].
In regards to direct detection experiments and the WIMP paradigm, given that the former test DM-nucleon cross sections, one may think that one way to avoid them is to propose a model in which DM is coupled to a leptophilic mediator [16], leading to vanishing scattering tree-level diagrams with nuclei. Direct detection searches have become however so constraining in the allowed scattering cross sections that even interactions with nuclei at 1-loop order (with leptons running in the loop) are under pressure [17][18][19][20]. In this work we propose to take the WIMP paradigm all the way to the end and so we consider a simple effective WIMP model that provides such small contribution to the DM scattering cross sections probed by direct detection experiments that they would naturally be buried under what is known as the neutrino scattering floor. We show that under simple assumptions about the nature of the interactions of the DM and its mediator to SM particles but sticking to masses and couplings of order EW, a model emerges in which second generation leptons and a Dirac fermion DM are charged vector-axially and vectorially respectively under a new U (1) µ spontaneously broken gauge symmetry. Under these conditions, the WIMP miracle is still accomplished and at the same time we are able to avoid any signal that could be measured by direct detection experiments at present and in the future (all the way up to the neutrino scattering floor). We find that indirect detection signals from current searches done by the Fermi-LAT collaboration from signals at the center of our own galaxy and from dwarf spheroidal galaxies do not currently provide strong constraints on the model. There can be however interesting collider searches that probe certain regions of the parameter space from current and future measurements at the LHC. There are also potential contributions to neutrino trident production and the magnetic dipole moment of the muon that though they are unrelated to the DM interactions, provide strong constraints to the mediator mass and coupling strength to muons and muonic neutrinos.
Furthermore, our effective model in particular being gauge anomalous provides a window into the UV physics that is responsible for curing the anomalies in the form of triple gauge boson couplings, as it has been analyzed in theories with an anomalous U (1) within different contexts [21][22][23]. By studying these anomalous interactions, we show that though it appears to be hopeless probe them at the high luminosity LHC (even hard at a 100 TeV hadron collider), they may be able to be probed at a future muon collider providing a clear signal and a great opportunity as a window into the UV physics.
The paper is organized as follows. In Section 2 we introduce the effective U (1) µ DM model describing its interactions. In Section 3 we discuss the experimental constraints on the model coming from DM relic density calculations, direct and indirect detection DM experiments, collider searches as well as contributions to (g − 2) µ and neutrino trident production. In Section 4, we focus on the projections at the 14 TeV LHC with luminosities of 300 fb −1 and 3000 fb −1 for the 3µ + E miss T channel and show that large regions of parameter space can be probed if a window for the invariant mass of the muon pair around the Z mass is added, though there are still regions in which the model remains elusive. Finally in Section 5, we discuss a very interesting feature of our model which is the presence of anomaly induced triple gauge boson couplings and their possible collider signatures at the LHC, and at a hypothetical muon collider or a 100 TeV hadron collider. Our conclusions are given in Section 6.

The U (1) µ axial model
We decide to pursue the idea of obtaining a model in which the WIMP miracle happens in a simple way, for couplings and masses of order the EW scale (or not too far from it). The lack of signal in direct detection experiments leads us to consider a leptophilic DM model with a massive vector mediator between the visible and dark sector, that couples vector-axially (VA) to leptons and vectorially (V) to the DM. These kind of interactions (VA with leptons and V with DM) have been shown to provide vanishing contributions at all loop-levels to the mixing between the Z and photons which tend to dominate spin independent DM-nuclei scattering [18][19][20] and only provide contributions to spin independent DM scattering via mass mixing between the Z and the Z gauge boson, proportional to the lepton's Yukawa couplings, which would put interaction with third generation leptons in tension. There are also e + e − → e + e − processes known as compositeness bounds at the Large Electron-Positron (LEP) collider that strongly constrain interactions with first generation leptons implying m Z 3 TeV in that case [20]. Thus we decide to consider vectoraxial interactions only with second generation leptons and vector interactions with the DM, which furthermore we choose to be Dirac such that in the annihilation process the s-wave contribution in the non-relativistic limit is non-vanishing.
Specifically, we introduce a new Abelian gauge symmetry U (1) µ under which second generation leptons are charged, such that the interaction in the mass basis for muons are axial and we also introduce a Dirac fermion χ charged vectorially under U (1) µ , which we identify with the DM; the rest of the SM fields including the Standard Model Higgs remain neutral under U (1) µ 1 . This new gauge symmetry is assumed to be spontaneously broken at some larger scale (possiblyà la Higgs, via the vacuum expectation value of a scalar field) and we work in the effective theory of the SM matter and gauge field content, with the addition of the Dirac DM χ and a massive vector gauge where Z µν = ∂ µ Z ν − ∂ ν Z µ is the Z -field strength, g χ = Q χ g , g µ = Q µ g are the coupling strengths for the DM and charged muon interactions with respective charges Q χ and Q µ under U (1) µ , g the U (1) µ coupling, and m χ and m Z are the masses of the DM and the Z gauge boson. The Z − ν µ coupling is fixed to preserve the EW symmetry. We assume a vanishing tree-level kinetic mixing term between the Z and the SM EW gauge bosons. In addition, we work in the scenario with m Z > 2m χ in which the phenomenology at the LHC is enriched. When m Z < 2m χ , the only invisible decay channel is Z → ν µνµ so that the invisible branching ratio is 1/3 irrespective of the values of m χ or g χ and the phenomenology of the LHC searches is disconnected from the DM parameters of the model. From the interactions in Eq. (1) we obtain the following expressions for the partial decay widths of the Z boson at leading order: where z x = m 2 x /m 2 Z . If the total width is saturated by the channels in Eqs.
(2)-(4), then its branching ratio can be obtained as is the branching ratio into states that appear at the LHC as missing transverse energy and can be written as, where we have used the definition ξ = g χ /g µ . From Eqs. (2)-(5), the ratio between the total decay width, Γ Z , and the Z mass, can be parametrized in terms of the coupling g µ and BR inv as follows, In the analysis of the model we will restrict ourselves to scenarios in which Γ Z /m Z < 0.3. It follows then that this expression is useful to translate this bound, for each value of BR inv , into a bound for the coupling g µ . The requirement Γ Z /m Z < 0.3 allows to have a not too wide Z while keeping available regions of the parameter space that are interesting from a phenomenological standpoint 3 . For example, as can be seen from Fig. 1, for BR inv as large as 0.9 we can still probe g µ values up to ∼ 1. Given that in particular the electromagnetically charged muon is also charged under the new U (1) µ gauge symmetry, one cannot write the usual Yukawa interaction for the muon responsible Figure 1: Values of Γ/m Z in the plane BR inv -g µ according to Eq. (6). In our analysis, we restrict ourselves to the region Γ Z /m Z < 0.3. via EW symmetry breaking to generate the muon mass since it would explicitly break U (1) µ . This interaction can however arise from a higher-dimensional operator involving the usual Yukawa interaction combined with a SM singlet Higgs field responsible for the spontaneous breaking of the U (1) µ gauge symmetry, which we shall denote Φ, with charge Q Φ with respect to U (1) µ and whose vacuum expectation value (vev) we denote Φ ≡ f v EW = 246 GeV sets the scale at which U (1) µ is broken 4 , where λ µ is a dimensionless coupling of order one and M f is the scale of the new physics that generates this operator. It is then clear that once Φ acquires a vev the low energy muon Yukawa takes the value y µ = λ µ (f /M ) −2Qµ/Q Φ , that can naturally accommodate the measured muon mass due to its smallness. We assume that this or a similar mechanism is responsible for the generation of the muon Yukawa and simply work with the low energy muonic Yukawa interaction included in L SM in Eq. (1).
Notice also that the fact that we have chosen a VA gauge interaction for which only the second generation leptons are charged implies naively the appearance of gauge anomalies which are known to lead to inconsistencies at the quantum level either by the breaking of unitarity and/or gauge invariance. Of course, the quantum field theory which completes our effective theory cannot have these anomalies in the UV and therefore there must exist additional fermions which cancel the possible gauge anomalies present. The effect of these additional fermions in the low energy effective theory manifest itself via the presence of triple gauge bosons couplings involving U (1) µ and the SM gauge bosons. We shall comment more on this in Section 5 and its possible phenomenological collider signatures in Sections 5.2 and 5.3.

Experimental constraints
In this section we study the conditions that must be met in order to reproduce the correct relic density of DM in the Universe today. Furthermore, we look into the possible DM standard signatures in direct and indirect detection experiments, as well as more specific signals that are particular of our model such as neutrino trident production and the muon anomalous magnetic moment. These latter affect indirectly the DM predictions by constraining the Z interaction with the second generation leptons. We also comment on the current collider constraints in the form of final states with muons and missing energy.

Relic density
The first main prediction we should satisfy is the DM relic abundance as measured by the Planck collaboration, Ω DM h 2 = ρ DM /ρ c = 0.120 ± 0.001 [27], where ρ DM is the DM density and ρ c the critical density of the Universe. As it is well-known, the WIMP miracle occurs within the standard freeze out mechanism, where initial conditions are washed out once the DM candidate enters in thermal equilibrium with the early Universe plasma. Given that the coupling strengths are assumed to be of order EW but not much smaller, it is clear that at earlier times in the Universe, the DM candidate χ can easily reach equilibrium with the thermal bath at temperature T via the process, where = µ, ν µ . The DM density then follows the thermal equilibrium density and as the temperature in the Universe drops due to the expansion it ends up following a Maxwell-Boltzmann distribution which leads to an exponential decrease in the DM density once the energy available is not sufficient to produce DM pairs from the thermal bath. However, once the annihilation rate is smaller than the Hubble expansion of the Universe at some temperature T F ≈ m χ /25, the DM particles fall out of equilibrium and decouple from the thermal bath, freezing at a density that gets diluted only by the expansion of the Universe. We will not demand that χ makes up all of the DM (though we will certainly consider this case as well), but instead we demand that it does not lead to a closed Universe by requiring the following condition to hold, which allows more freedom for our parameter space by assuming that some other source completes the full DM content of the Universe. All of this is done by integrating Boltzmann equation numerically and computing Ω χ h 2 for a particular point in our parameter space. We perform this automatically by using the publicly available micrOMEGAs code [28].

Direct detection
Direct detection experiments aim to detect DM particles present in the gravitationally bound halo of our galaxy. Through scattering with subatomic particles, DM could leave a signal in underground detectors, in the form of recoil energy. The most stringent bounds on DM-nucleon scattering cross sections come from these experiments, in particular for spin independent (SI) interactions for which the scattering cross sections are enhanced by a coherent sum of the nucleons in the nucleus. This leads to strong constraints for the possible effective interactions of WIMPS with quarks as shown in [1][2][3]29]. In order to avoid these bounds we have chosen as mentioned before a leptophilic Z mediator with axial couplings to muons, which guarantees that not only tree level diagrams are not present, but also that Z -photon mixing vanishes at all loop orders 5 . There remains however an effective operator that provides a contribution to both spin independent and spin dependent (SD) interactions and involves Z-Z mass mixing [20] 6 , where g H is a loop-induced effective coupling proportional to the Yukawa coupling of the leptons that run in the loop and are charged under U (1) µ , so its contribution is naturally suppressed (however, notice that a model with the tau charged under the U (1) symmetry is significantly constrained [20]). One can calculate the corresponding quark-DM contact interactions that provide contributions to SI and SD interactions respectively, where the vector and axial coefficients C A can be obtained from the renormalization group equations. We apply to our model the results given in [20], where the running of the coupling constants is performed from the energy cutoff of the effective leptophilic model, Λ U V ∼ f , down to the energy scale µ N associated with nuclear processes. These coefficients read, where y µ is the Yukawa coupling for the muon, T and Q q are the weak isospin and charge of the quark q respectively, and s W is the sine of the Weinberg angle θ W .
Only the first contact term of Eq. (11) contributes to a spin-independent signal, and thus we focus on this term only. The amplitude for the scattering of χ by a nucleon N is then calculated as, When evaluating matrix elements of quarks in a nucleon N at zero momentum transfer, only valence quarks amount for the result [30], where the first form factor is approximated by the number of valence quarks n q in the nucleon, and the second term is dropped in the q 2 = 0 approximation. Also, in the non-relativistic regime we can write for both nucleon and DM parts, where ξ is a two-component spinor and s and s are the initial and final spin states. Squaring the amplitude, averaging over initial states, summing over final states and integrating in the phase space we obtain the total cross section, 5 The largest contributions to the scattering cross section comes from the Z -photon mixing. 6 Loop-induced Z − Z kinetic mixing is also possible, but it is suppressed with respect to mass mixing by (q/MZ ) 2 /y 2 µ ∼ (MeV/100 GeV) 2 /y 2 µ ∼ O(10 −4 ), where q is the momentum transfer in the scattering [18]. and by setting n u = 2, n d = 1 we can obtain the χ-proton cross section. Finally, we evaluate the C (q) V coefficients at the benchmark values g µ = g χ = 1, m Z = 200 GeV. The nuclear scale is set at µ N = 1 GeV where the hadronic matrix elements are evaluated, and since we expect heavier states to appear at few TeV, we set for practical purposes Λ U V = 10 TeV, though the result only depends mildly on the scale Λ U V . For the range of χ masses evaluated in this work, this results in a χ-proton cross section of, and thus it is unconstrained by current experimental research, whose more restrictive upper bounds are around 10 −47 cm 2 [31]. Moreover, this cross section value sits well below the neutrino floor ∼ 10 −50 cm 2 , inside a region of the parameter space which seems to be inaccessible by conventional SI direct detection experiments 7 .

Indirect detection
In our Universe today, dark matter annihilation is heavily suppressed by its low abundance. Thus it can only be observed in high density regions, such as galactic centers or DM-dominated dwarf galaxies. For our case, muon pairs are produced from χχ annihilation, which then lead to a continuum spectrum in the form of bremsstrahlung radiation, whose differential flux per unit energy per unit solid angle is given by, where σv χ is the thermally-averaged χ-annihilation cross section into muons, dN/dE is the radiation spectrum generated by the muons, and ρ χ is the χ-density which must be integrated along the line of sight. The Fermi-LAT collaboration is able to measure the γ-ray sky via its satellite experiment with unprecedented precision. In particular, it provides strong constraints on possible DM annihilation coming from the center of our Milky Way [8] and also from a series of dwarf spheroidal galaxies [9,10], which are expected to be DM dominated. Indirect detection experiments typically give constraints in the m χ , σv χ plane, assuming χ accounts for all the DM in the Universe. Since as mentioned before we relax such assumption and only demand that the relic density should not overclose the Universe, in order to translate the constraints from the Fermi-LAT results to our particular model, we rescale ρ χ as [34], where Ω DM h 2 = 0.120 and ρ DM is the total DM density used in the indirect detection analysis. If, for a given m χ , the cross section is bounded above by σv lim , then this corresponds to the following limit in our model, We use the values for σv lim as given in Fig. 10 of [35] at 95% Confidence Level (CL), taking the more restrictive Bayesian limit. The relic abundance and cross section for our candidate χ are calculated using micrOMEGAs. For the latter we assume it depends on the model parameters in a straightforward manner, as we expect for the χχ → µ + µ − tree level dominant process, where F encapsulates the portion of the thermally-averaged cross section that does not depend on the coupling parameters. It is also expected that the relic density is inversely proportional to the annihilation cross section σv ann . Since χ can only annihilate to muons or muonic neutrinos, which interact with the Z boson through the same coupling strength g µ , the annihilation cross section will exhibit the same dependence on the couplings as in Eq. (22). We assume then, where Ω 1 (m χ , m Z )h 2 is the relic density for given masses and g µ = g χ = 1. Substituting equations (22) and (23) into (21) and setting g χ = ξg µ , By calculating the factors F and Ω 1 h 2 for given masses and g χ = g µ = 1 inside micrOMEGAs, using the experimental limits extracted from [35], and setting Ω DM h 2 = 0.120, we can obtain the exclusion limits in the m Z − g µ plane as in the plots below.
We make a brief comment regarding the constraints coming from cosmic-ray positrons as measured in particular by the PAMELA [36] and AMS-02 [37] experiments. It was shown in [38,39] that DM models with annihilations into µ + µ − which reproduce the DM thermal relic density via freeze-out can be strongly constrained by these measurements. In particular, in the recent study of [39], it is stated that DM masses below ∼ 150 GeV would be excluded for this channel. There are however large uncertainties stemming from the local dark matter density and local magnetic field strength values used, as well as a full lack of knowledge of the characterization of the galactic magnetic field and its effects in the propagation of charged particles, all of which can lead to softer constrains on the DM masses that could potentially be excluded [40]. Due to these large uncertainties, we will only consider for indirect detection constraints the exclusion limits provided by the γ-ray observations from the annihilation of DM in dSph galaxies, though there are certainly points of our parameter space that could satisfy the cosmic positron constraints as well.

Muon anomalous magnetic moment
Recent measurements at Brookhaven National Laboratory and at Fermilab [41,42] of the muon magnetic moment result in a combined excess of 4.2 σ with respect to the SM predicted value (see [43] and references therein). Due to the axial nature of the Z coupling to charged muons, there is a negative one loop contribution to the muon magnetic moment [20] with respect to the SM contribution, which worsens the agreement with respect to the measured value. This thus implies a constraint on the possible Z mass and coupling to the visible sector in our model and indirectly a constraint on the DM as well. We demand that the new physics contribution does not exceed the combined theoretical and experimental uncertainties at the 2σ-level, where σ exp = 41 × 10 −11 and σ th = 43 × 10 −11 . Notice that with this requirement the NP contribution of Eq. (25) is hidden in the theoretical and experimental uncertainties, keeping the disagreement with the measured value below the 5σ level.
It is worth mentioning that lattice QCD calculations have been shown to modify the SM theoretical predictions [44], reducing the discrepancy with respect to the experimental measurements to a 1.6σ excess [45], which would also reduce the tension within our model.

Trident diagrams
The Z interaction with muonic neutrinos are constrained by what are known as neutrino trident experiments, which probe lepton production by neutrino scattering with nucleons. In our model, the process ν µ N → ν µ µ + µ − N takes place with the contribution of Z exchange diagrams. Measurements of the associated cross section at the Columbia-Chicago-Fermilab-Rochester (CCFR) neutrino experiment at the Fermilab Tevatron [46,47] imply the following constraint, where v 246 GeV is the EW breaking vacuum expectation value. Interestingly enough, for the case of only axial contributions as in our scenario, it is expected that future Deep Underground Neutrino Experiment (DUNE) will have similar sensitivity as CCFR and therefore no improvement of this particular constraint is expected.

Collider searches at the LHC
At the LHC the leptophilic Z boson is mainly produced through Drell-Yan processes with the Z radiated from one of the final leptons, as shown in Fig. 2. The specific final state depends on which is the exchanged EW gauge boson and the decay mode of the Z . For the muonic Z the exchange of γ and Z bosons lead to final states with 4µ or 2µ+ E miss T [48], while the exchange of W gives rise to signals with 3µ+E miss T or 1µ+E miss T . Note that the missing transverse energy in the modes with one or two muons can arise not only from the neutrinos but also from the DM particles. According to the analyses in [49][50][51] the constraints coming from the modes with at least three leptons are the most restrictive ones for the case of a leptophilic vector boson. Moreover, it is also shown in [50] that in general the 3µ+E miss T mode provides stronger bounds than the 4µ mode, except when the Z coupling to right-handed muons (g R ) is larger than to left-handed muons (g L ), with the limits still being similar for g R /g L as large as 2. In our case, since the Z is VA coupled to muons, we have g R = −g L , which translates to cross sections for the 3µ+E miss T channel that are ∼ 3 times those of the 4µ mode along the parameter space explored here.
In order to obtain an estimation of the limits imposed by the existing searches at the LHC, we confront the most promising channel of the model, 3µ+E miss T , with the collection of LHC analyses implemented in the Public Analysis Database (PAD) of MadAnalysis5 [52][53][54]. We also include the limits coming from the 2µ+E miss T channel since in the model we are considering its cross section increases with g χ and then it is interesting to test if there is some regime in which it gives better exclusion bounds than the 3µ+E miss T mode. The signals of both channels were simulated with MadGraph 3.1.1 [55] using an implementation of the model in UFO format obtained with Feynrules [56]. The parton shower and hadronization were carried out with Pythia 8.2 [57], while a fast detector simulation is performed inside MadAnalysis5 by Delphes 3 [58]. In Fig. 3 we show the exclusion limits in the mass range (200 − 500) GeV 8 arising from the most restrictive search of the PAD along with the bounds imposed by the CCFR neutrino experiment and the latest measurements of the muon magnetic moment at Fermilab. We show two representative cases: g χ = 0 (left panel) and g χ = 0 with BR inv = 0.6 and z χ = 0.01. For the latter we also include limits from indirect detection and relic density. The most restrictive limit for the 2µ + E miss T channel corresponds to the inclusive signal region SR SF 1J with m T 2 > 160 GeV of the ATLAS search for electroweakinos and sleptons in the 2 + E miss T channel [59,60]. For the 3µ+E miss T mode the most constraining bound is obtained with the signal region A44 of the CMS search for electroweakinos in multileptons final states [61,62]. As can be seen from Fig. 3, none of these bounds are competitive with those arising from the trident diagrams and the muon anomalous magnetic moment. Regarding the comparison between the 2µ+E miss T and 3µ+E miss T signatures, we see that the latter provides the stronger limits as expected. However, the opening of the DM decay mode of the Z slightly improves the limit provided by the 2µ+E miss T and at the same time worsens the 3µ+E miss T bound. This is also expected due to the behavior of the corresponding cross sections when the decay channel Z → χχ is available. Nevertheless, the gap between the exclusion limits is still significant and we have checked that they become comparable only when BR inv is as large as 0.9. For this reason we will focus exclusively on the 3µ+E miss T signature in what follows. The lack of sensitivity of the existing searches at the LHC to exclude and/or discover the model studied here is in part due to the fact that the signal regions mentioned above select events with both electrons and muons in the final state instead of only muons, which would be better suited to probe a muonic Z . Thus, in order to have a better sense of the LHC potential to test the U (1) µ axial model it is necessary to consider a dedicated search strategy. Here we apply the one proposed in [50], which is based on the ATLAS search [63]. In the following list we summarize the cuts defining the search strategy: 1. A set of basic cuts are applied to the muons and jets at generator and detector level: p µ T > 50 GeV, |η µ | < 2.4, ∆R jµ > 0.4, p j T > 20 GeV and |η j | < 2.5. 2. A b-jet veto is applied (N b = 0) and only events with exactly three muons (N µ = 3) with net charge ±1 are kept.
3. Two more vetoes are applied in terms of requirements on the invariant mass of the opposite sign muon pairs: m µ + µ − > 12 GeV (low resonance veto) and |m µ + µ − − m Z | > 10 GeV (Z veto).
4. The missing transverse energy is required to be above 100 GeV.
5. The transverse mass of the muon not belonging to the pair that better reconstructs the Z mass must be larger than 110 GeV.  For the latter we also include the regions excluded by the relic density and the dSph galaxies bound from Fermi-LAT. We add for reference the lines corresponding to Γ Z /m Z ratios 0.1 and 0.3.
6. Finally, at least one opposite sign muon pair must fulfill the requirement As it is pointed out in [50], the last cut is crucial to enhance the sensitivity and obtain significant limits. Notice that this cut would not be suitable for the 2µ + E miss T channel where part of the signal events arises from the Z µ + µ − production followed by the Z invisible decay, which makes impossible the reconstruction of the Z boson. This is another reason to prefer the 3µ + E miss T channel instead.
In order to obtain the prospects for discovery and exclusion derived from this search strategy, we simulated the dominant backgrounds using the same simulation setup already described for the simulation of the signal. We considered the diboson processes leading to ν and final states, where , = µ, τ , and, in addition, the triboson processes W W W, W W Z, W ZZ and ZZZ that produce final states with the same charged lepton content. Since the search strategy described above is based on the signal region SRnoZc of [63], we used the corresponding background rates to validate our simulation procedure, finding good agreement within the reported uncertainties. The results obtained by applying the dedicated search strategy to the most promising channel will be presented in Section 4 along with the experimental constraints discussed previously in this section.

LHC projections for Run 3 and high luminosity
In order to estimate the exclusion/discovery prospects for the U (1) µ axial model at the LHC with the search strategy described in Section 3.6, we will consider the four benchmarks listed in Table 1. Each benchmark is defined by its invisible branching ratio BR inv and squared mass ratio z χ = m 2 χ /m 2 Z values, which also fix the coupling ratio ξ = g χ /g µ in virtue of Eq. benchmark we generate pp → 3µ + E miss T signal events varying the coupling g µ and the mass m Z , while adjusting g χ and m χ to the corresponding values that fix BR inv and z χ . We use the same UFO model as the one implemented for the PAD recasting described in Section 3.6. Parton level events are generated using MadGraph, parton shower and hadronization are simulated by Pythia 8.2, and fast detector simulation is done by Delphes 3. Detector-level events are then passed through the analysis described in Section 3.6, which is implemented in MadAnalysis5, where we define a single signal region (SR) for the events that pass all cuts. We obtain the detector-level number of events that pass the SR cuts as where L is the total integrated luminosity, σ is the process cross section and A is the selection acceptance. The latter is calculated as the fraction of events that pass all the cuts over the total number of generated events. Eq. (28) is used for both signal and background processes. The 95% C.L. exclusion limits and the discovery prospects are obtained [64] by demanding and respectively, where s and b denote the number of signal and background events that pass the SR cuts. We do not include systematic uncertainties in our analysis. We show in Fig. 4 the exclusion and discovery limits in dashed contours, along with the other experimental constraints described in Section 3, for the Z mass range of (200 − 500) GeV. We use L = 300 fb −1 as the LHC is expected to reach this value in the near future, and L = 3000 fb −1 to account for the full LHC lifetime. We also plot Γ Z /m Z reference values, as we restrict our analysis to Γ Z /m Z < 0.3. Note that, for every benchmark considered, there is a portion of the parameter space that remains unconstrained by current experimental data. These available regions are bounded below by the relic density constraint, for which the gray area in the plots represent points that predict a χ relic density that would overclose the Universe (i.e. Ω χ > Ω DM ), and therefore cannot be allowed. On the other hand, upper bounds can be taken up to the more restrictive magenta lines, below which our model predicts a (g − 2) µ correction that lies within the 2σ experimental uncertainty, as described in Section 3.4. Taking into account that there could be other new physics sources for this quantity that cancel the Z exchange contribution, but more importantly that there are lattice results [44,45] that show a smaller discrepancy between the SM prediction and experimental values, we could relax the muon anomalous magnetic moment bound in the plots and consider the region below the trident constraint to be explored by our collider analysis.
For the benchmark point BMI there is a region of the parameter space above m Z ∼ 350 GeV that is unconstrained by current experimental bounds. The 95% C.L. exclusion limits we obtain for the collider analysis at 300 fb −1 would exclude our model up to 450 GeV for this benchmark.
Moreover, the 3000 fb −1 exclusion and discovery prospects cover the available parameter space in the range of masses considered. The current experimental bounds evaluated at the benchmark point BMII allow for m Z values along the whole span of masses considered. A portion of this region could be explored at 300 fb −1 , improving the (g − 2) µ limit in the (200 − 500) GeV range. At 3000 fb −1 for BMII it would be possible to completely exclude our model for m Z < 300 GeV, and most of the g µ range can be excluded for larger masses (for example, with m Z = 500 GeV, g µ 0.35 is excluded, which corresponds to Γ Z /m Z 8 × 10 −3 ). Discovery prospects at the LHC for this benchmark also cover a portion of the available parameter space in the range of masses considered. At benchmark point BMIII our 300 fb −1 exclusion limit is weaker than the (g − 2) µ constraint. As we mentioned before, this constraint can be relaxed, so our collider limit at 300 fb −1 shows an improvement with respect to the trident bound below 280 GeV, becoming no longer competitive for larger masses. The discovery prospect also sits between (g − 2) µ and trident limits for masses in the (230 − 320) GeV range, while being slightly below (g − 2) µ bound for lower masses. The 3000 fb −1 exclusion limit we obtain fully covers the available parameter space for m Z < 300 GeV, and improves the (g − 2) µ limit for masses up to 430 GeV, where it excludes g µ 0.7 (Γ Z /m Z 0.13). Collider limits in benchmark point BMIV show a similar behavior to BMIII. For example, at m Z = 200 GeV, g µ 0.17 (Γ Z /m Z 7.7 × 10 −3 ) is excluded at the HL-LHC, while for m Z = 430 GeV it excludes g µ 0.7 as before. However, since lower masses are allowed by the relic density bound in comparison with BMIII, there is a larger unconstrained region in this benchmark in the whole (200 − 500) GeV mass range, even at 3000 fb −1 .
Note that for benchmark points BMI and BMIII the relic density lower bound is higher than in BMII and BMIV. This is due to the fact that their z χ value is lower, further away from the resonant production, thus needing a larger coupling to reproduce the annihilation cross section value associated with the measured DM relic density. The same reasoning applies to the dSph constraints, as they are also related to the χχ annihilation to SM states via Z (see Section 3.3). Also, collider limits obtained for BMI and BMII are stronger since the Z branching ratio to muons, i.e. 1 − BR inv , is larger. In these scenarios it would be possible to improve the bounds for our model in the near future at the LHC. For a Z with larger BR inv (BMIII and BMIV) one would need the estimated luminosity of the full LHC lifetime to probe a significant portion of the available parameter space, while still leaving some areas inaccessible to the LHC. Finally, it is worth noting that collider limits shown for BMIII and BMIV become steeper for higher masses and couplings. This is due to the fact that Γ Z becomes larger for these values, as we can see in Eq. (6). Therefore, the mass window cut (cut 6 of the analysis in Section 3.6) becomes less effective as the width of the resonance increases and more signal events are left out of the signal region. We tried to substitute the mass window cut with other similar cuts that attempt to include more signal muons, but we only observed a small improvement in the limits for a "wider window", using the condition |m µ + µ − − m Z | < 0.3 m Z instead of cut 6. Nonetheless, the regions where the wider window limits are competitive with the ones in Fig. 4 are above trident and/or (g − 2) µ bounds, and thus are already excluded. It is also worth noting that in all of the cases described the limits in which we include a m µ + µ − window around the m Z value are stronger than the ones obtained without a mass window. We show some examples of cutflows for this analysis in the Appendix B.

Mixed gauge anomalies
In the model under consideration the axial (chiral) coupling of the muon (muonic neutrino) to the Z boson induce mixed anomalies between U (1) µ and the EW gauge symmetry group and gravity. Gauge symmetries of the classical fields cannot be simultaneously satisfied in the QFT and this is signaled by anomalous Ward identities which imply the loss of unitarity and/or Lorentz invariance in the first place and the non-renormalizabilty of the theory in second place [65]. By allowing the gauge field to acquire mass (as we are doing via spontaneous symmetry breaking) one naturally solves the problems of unitarity and Lorentz invariance but the issue of renormalizability remains and implies that the theory can only be regarded as an effective theory with a cutoff that cannot be made arbitrarily large without suffering a loss of calculability [65]. In a fully consistent QFT all anomalies must be canceled and in fact, anomalous theories can be regarded as effective models where the anomaly cancellation fully happens at a higher energy scale and in which part of the spectrum responsible for anomaly cancellation has been integrated out in the effective theory 9 .
In the case of our anomalous U (1) µ symmetry, the Ward identities involving the U (1) µ and EW gauge symmetry groups in the EW unbroken phase can be accommodated such that, where the SU (2) L ×U (1) Y gauge bosons are labeled collectively as A i and p 1 , p 2 are their momenta. The vertex functions Γ µνρ A a A b Z corresponding to triangle diagrams with Z , A a and A b in the external legs are not invariant under a shift of the momenta of the fermions running in the loop. A particular shift to restore the Ward identities for the SM gauge group can be always chosen leading to Eqs. (31), (32). However, there is no possible shift that satisfies simultaneously the three Ward identities and thus Eq. (33) remains non-vanishing. The coefficients of the mixed anomalies between the Z and the SM gauge bosons can be calculated as, with Q , T i being the U (1) µ and EW gauge group generators, respectively, and Tr R R(L) standing for the symmetric part of the trace evaluated in the right-handed (left-handed) chiral representation of the SM fermions running in the loop. The anomalous coefficients A Z Z B , A Z Z W a and A Z BW a are zero because the EW generators are traceless and factor out in Eq. (34). The non-vanishing anomalous coefficients are, The presence of 1-loop anomalous triple gauge couplings [66], in particular the ones involving the Z and the EW gauge bosons, renders the proposed model interesting from a phenomenological perspective. Indeed, we can consider these interactions as remnants at low energy of the UV physics which completes our effective theory and thus, if we were able to probe these couplings, we would be having access through a window into this UV physics. We shall discuss in Sections 5.2 and 5.3 the potential reach of searches for anomalous couplings at the LHC and at hypothetical 100 TeV proton and muon colliders.

Example of an anomaly-free UV completion
We briefly discuss now a possible anomaly-free UV completion of our effective model. With this aim we introduce new fermions in order to cancel all the gauge and gravity anomalies in the high

Hadron collider searches for anomalous triple gauge boson couplings
As it has been discussed throughout the text, our model represents an anomalous EFT which can be interpreted as a UV complete theory in which part of the chiral fermion spectrum charged under the U (1) µ has been integrated out. This implies in particular an enhancement with the energy of the process in the coupling between the longitudinal mode of the Z and two EW gauge bosons [21,23], as can be seen from the anomalous Ward identities, signaling the breaking of unitarity. Such behavior of course is tamed at high energies by the appearance of the spectator fermions and the restoration of a full unitary theory in the UV, though at intermediate energies the enhancement remains and could potentially soften the 1-loop suppression of the anomalous couplings, providing hope of probing these at current and future particle accelerators.
In Fig. 5 we show two situations in which the anomalous triple gauge coupling could enter in the production of a Z at hadron colliders. The first one corresponds to a vector boson fusion process and as such, due to the t-channel production of the Z , we expect no enhancement with the energy in the process. The second process on the other hand is an s-channel Z -strahlung production, where the Z is emitted from an EW gauge boson in the s-channel. In this case, due to the s-channel production of the intermediate EW gauge boson, the Z can potentially be sufficiently boosted such that the enhancement with energy in the production of the longitudinal mode of the Z leads to testable cross sections. Since we expect the longitudinal mode of the Z to dominate the production cross sections in the high energy regime and in order to avoid the complications from having to integrate over the anomalous momentum-dependent couplings convoluted with the proton's PDF, we can further simplify our calculations using the Goldstone equivalence theorem to estimate the amplitude of the process of interest, where L here denotes the longitudinal polarization of the Z and φ is the corresponding eaten Goldstone boson. We can obtain the anomalous coupling of the Goldstone boson to the EW gauge bosons by replacing in the tree-level Lagrangian Eq. (1) the interaction term between the Z and the fermionic currents, where f is the Goldstone decay constant which coincides with the vev of Φ. Integrating by parts, we can use that the divergence of the anomalous fermion current is related to the anomaly as [68], where in our case the field strengths correspond to the EW sector SU (2) L × U (1) Y . Replacing in Eq. (38) and using the values of Eq. (35) we obtain, where g 1 and g 2 are the U (1) Y and SU (2) L gauge couplings. Rewriting the last expression in the EW broken physical base (γ, W ± , Z) we get, where s W y c W are the sine and cosine of the Weinberg angle and for each vectorial field C µ we define C µν = ∂ µ C ν − ∂ ν C µ and its dualC µν = µνρσ C ρσ /2. Note that effective vertices between φ and two or three EW gauge bosons are generated, while vertices between φ and four gauge bosons vanish as a consequence of pentagonal diagrams being not anomalous [68]. The three and four boson anomalous couplings given in Eq. (41) are then implemented in a UFO model via Feynrules, which we use to calculate parton-level cross sections using MadGraph. We simulate the pp → φ V signals 10 as in diagram 5b, with V = γ, Z, W ± and the model parameters fixed at g µ = 0.46, g χ = 0.1, m Z = 200 GeV, m χ = 30 GeV (that correspond to ξ ≈ 0.22, z χ ≈ 0.06, BR inv ≈ 0.35), which satisfy Ω χ h 2 ≈ 0.1 and sit at the trident bound described in Section 3.5, aiming to maximize the anomalous couplings. We get the largest cross section for the pp → φ γ channel, resulting in σ φ γ = 5.9 × 10 −3 fb, leading to an approximate of 20 events produced at the LHC at √ s = 14 TeV and at L = 3000 fb −1 , which is estimated for its full lifetime.
The cross sections found for the Z and W channels are smaller by a factor of at least 3, resulting in less than 10 events. For a dedicated search of this signal it is necessary to impose cuts on the final state particles, further reducing the number of expected events, so probing these anomalous vertices is extremely challenging at the LHC. In a hypothetical √ s = 100 TeV hadron collider we obtain σ = 7.6 × 10 −2 fb, giving ∼ 230 events by assuming an integrated luminosity of 3000 fb −1 . This might be a more promising scenario for probing the anomalous vertices, though the backgrounds are also expected to grow, but a much more favorable experimental facility for our model would be a muon collider, in which we focus the next section.

Muon collider resonant searches for anomalous triple gauge boson couplings
Muon colliders provide a great potential to explore new physics in the sub-TeV to the multi-TeV energy range [69][70][71]. The large mass of the muon in comparison with the electron mass suppresses synchrotron radiation by roughly a factor of 10 9 for beams of the same energy, and therefore rings can be used to accelerate muon beams efficiently and bring them repeatedly into collision. Furthermore, the physics reach of a muon collider extends that of a proton-proton collider of the same energy since all of the beam energy is available for the hard collision, whereas a fraction of the proton-beam energy is carried by the colliding partons 11 . A dedicated muon collider can scan the Higgs resonance and precisely measure its mass and width [72][73][74]. In fact, a muon collider is ideal to search for new physics and for resolving narrow resonances both as a precision and/or as an exploratory machine. There are nonetheless challenges that arise from the short muon lifetime and the difficulty of producing large numbers of muons, which requires the development of demanding technologies and new concepts. The beam background from the muon's decay has also consequences on both machine and detector design. An ambitious research and development program is needed to assess the feasibility of a muon collider in the tens of TeV range [75,76]. Therefore it will be important to study the physics potential of smaller-scale machines in the sub-TeV range that may be built along the way as technology demonstrators. In this respect, a sub-TeV muon collider becomes the ideal machine not only to measure with starking precision by resonant production the main decays of the U (1) µ Z into muon pairs or into invisible particles (muonic neutrinos or DM), but moreover, it may even be possible to probe the triple gauge couplings stemming from the mixed anomalies.
Before exploring the possibility of probing the triple gauge couplings at a muon collider, we first provide an estimation at Monte Carlo truth level of the discovery prospects for our model given by processes at tree level. In particular, we study two channels which would give the most promising signals: the dimuon channel µ + µ − → µ + µ − and the monophoton channel µ + µ − → γ + E miss T . For each of these processes we generate signal and background events at parton level with MadGraph. As mentioned before, it is expected that a sub-TeV muon collider starts operating at center of mass energy values near the SM Higgs mass, so we fix √ s = 125 GeV in our simulations. This implies that, for the range of Z masses considered, the channels studied are mediated by an off-shell Z . We estimate the luminosity required to achieve a 5σ signal significance for discovery (see Eq. (30)).
In Fig. 6a we show the required luminosity for the discovery of the dimuon production process at different g µ and m Z values 12 . We impose the MadGraph default transverse momentum cut for leptons, p T > 10 GeV, which is consistent with pre-selection cuts applied on full simulations of muon colliders [77]. On the plot we only show points with (m Z , g µ ) values that are unconstrained by (g − 2) µ measurements. For integrated luminosities of the order of a few tens of fb −1 , discovery would be possible for g µ = 0.6 in the available region of the m Z range considered. Furthermore, at L ∼ 300 fb −1 , discovery can be achieved for: g µ = 0.5 in the available m Z range, g µ = 0.4 with m Z 440 GeV, or even g µ = 0.3 and m Z 325 GeV. Note that for L ∼ O(few) ab −1 one could even probe couplings as small as g µ = 0.2 up to masses of order m Z 250 GeV. For the case of the monophoton channel, we show in Fig. 6b the target luminosities for the benchmark points provided in Table 1. For each of these benchmarks we fix g µ in order to maximize the m Z range that is unconstrained by the relic density and (g − 2) µ bounds from Fig. 4. Since the dominant background of this channel is the on-shell γZ production with invisible Z decay, and therefore the photon energy is fixed at E γ ≈ 29 GeV, we impose a cut in the transverse photon momentum p γ T > 30 GeV in order to suppress this background contribution. We find that the only realistic scenario for discovery of this channel is BMIII, where the z χ value is low enough to keep the χχ production open at m Z 450 GeV (m χ 45 GeV), and at the same time the relic density and (g − 2) µ bounds allow for a wide m Z range to be probed. Even in this case, only Z masses below 350 GeV can be discovered with L < 10 ab −1 . Note that the BMI curve would show a similar behavior as BMIII, with a considerable decrease in the target luminosity at lower masses, but since (g − 2) µ excludes m Z 400 GeV at g µ = 0.65 these points are not shown. On the other hand, BMII and BMIV have m χ > 80 GeV, which does not allow DM in the final state at √ s = 125 GeV, and therefore the signal cross section is significantly reduced, depending only on m Z and g µ . In fact, this effect can be seen by comparing the BMII and BMIII curves, both at g µ = 0.4: when the DM channel opens for masses m Z 450 GeV in BMIII the required luminosity rapidly drops for lower masses in comparison to BMII, in which only neutrinos contribute to the invisible decay channel. The luminosities required for discovery at benchmarks I, II and IV are above 60 ab −1 , which is considered to be too large for a muon collider, in particular for one with the fixed energy value of √ s = 125 GeV.
Comparing both channels it is clear that the sensitivity for the monophoton channel is significantly lower than the dimuon channel. If a hypotetical sub-TeV muon collider is able to perform a collision energy scan, it could eventually probe energies that allow on-shell Z production, which would significantly improve the discovery prospects for both channels in the mass range considered above. A muon collider is also ideal to probe m Z values higher than 500 GeV, where the hadronic collider searches at the 14 TeV LHC lose sensitivity.
In the following, we work in a region of parameter space under the assumption that the Z gauge boson has been already discovered at any of the possible experiments in which it can be searched for and focus on probing the triple gauge couplings stemming from the mixed anomalies, showing that indeed this could be achievable by resonantly producing the Z at a future muon collider 13 .
Close to the resonance, the cross section corresponding to the s-channel production of a Z boson decaying into a two-body final state consisting of particles C and D, can be written as  However, since we are interested in the resonant production of the Z , by virtue of the Landau-Yang theorem, the width of the γγ channel is zero, and so the resonant cross section vanishes for this channel. On the other hand, we expect the prospects of the W + W − channel to be similar to those of the ZZ channel, so that we focus in what follows on the Zγ and ZZ final states. By approximating the energy spectrum of each beam by Gaussian shapes, the effective cross section at the muon collider,σ Z ( √ s), can be computed by convoluting σ Z ( √ŝ ) with the Gaussian distribution in √ŝ centered at √ s and with a standard deviation given by [73], where R is the resolution in the energy of the muon beams. For a proton driver muon facility this resolution is around 0.004% for c.m. energy of 126 GeV and increases up to 0.1% in the multi-TeV range [75]. Since we explore Z masses between 200 and 500 GeV, we will take the former as reference value. The ratio Γ Z /m Z is at least one order of magnitude larger than σ √ s=m Z /m Z except for couplings g µ as small as 0.08 15 . Therefore, it is reasonable to use the approximation Γ Z σ √ s=m Z to compute the resonant effective cross section at √ s = m Z . In addition, in the regime Γ Z σ √ s=m Z the effective cross section turns out to be independent of the resolution R. The corresponding expression is given by, Figure 7: s-channel diagram for the production of a Z boson at muon colliders. The black blob represents the triple gauge couplings arising from the mixed anomalies.
where CD = Zγ or ZZ, with m C = m Z , m D = 0 or m C = m D = m Z , respectively. From Eq. (44) we see that the effective cross section depends on the branching ratios of Z → µ + µ − and Z → Zγ/ZZ. The former can be approximated by 1−BR inv , with BR inv given in Eq. (5), due to the smallness of the anomalous decays, while the latter can be obtained by computing the partial width corresponding to the anomalous decay channel and approximating the total decay width by Eq. (6). Therefore, in contrast to the processes studied in the previous section for the hadron collider, in this case we do not need to rely on the Goldstone equivalence theorem to estimate the cross sections.
The decay rate for the process Z → Zγ is where p = (m 2 Z − m 2 Z )/2m Z is the momentum of the final state vectors in the center-of-mass frame and |M| 2 is the spin-averaged squared of the matrix element which, following the procedure described in [78][79][80][81], is given by where e = g 1 g 2 / g 2 1 + g 2 2 and g Z = g 2 1 + g 2 2 , with g 1 and g 2 the SM gauge coupling constants of U (1) and SU (2), respectively, and t Z Zγ µ = (1 − 4 sin 2 θ W )Q µ . The expressions of the integrals I iµ can be read from the general formulae provided in the Appendix A.
For the process Z → ZZ the momentum of the outgoing Z bosons is p = m 2 Z − 4m 2 Z /2 and the spin-averaged squared of the matrix element is in this case given by [78,79], where  By computing the effective cross sections of the ZZ and Zγ channel from Eq. (44), we found that the former is between one and two orders of magnitude larger within the range of Z masses considered here. For example, with BR inv = 0.34 we obtain cross sections in the range of (7.5 − 25.8) fb for the ZZ channel, while for the Zγ channel these values are within (0.07 − 0.59) fb. In addition,the cross section we obtain for the irreducible ZZ background is in the range of (0.41 − 1.31) pb for the √ s values considered, which is significantly lower than the Zγ cross section values of (10.6 − 79.0) pb. Then, the detection of the Zγ channel appears to be very challenging in comparison with the ZZ channel. For this reason in the following we concentrate only in this channel. In order to provide an estimation of the discovery prospects we compute the luminosity required for the significance s/ √ b to be 5, with the number of background events obtained from the cross section of the irreducible background simulated with MadGraph 16 . The results are shown in Fig. 8 for BR inv of 0.34, 0.6 and 0.7. It is clear that the required luminosity increases with BR inv . This behavior is expected: on the one hand, the effective cross section is proportional to BR(Z → µ + µ − ) = 1−BR inv , and on top of this the total width of the Z used to compute BR(Z → ZZ) increases with BR inv as can be seen by rewritting Eq. (6) as Γ Z = g 2 µ m Z /(12π(1 − BR inv )). Therefore, the most promising scenario corresponds to a invisible decay rate dominated by the neutrino final state.
For each value of BR inv , the minimum required luminosity is reached at m Z = 300 GeV which is consistent with the fact that the anomalous partial decay width Γ(Z → ZZ) exhibits a peak around that mass. For the most promising scenario, luminosities between ∼ 30 fb −1 and ∼ 60 fb −1 would be enough to detect the anomalous decay channel for Z masses above 250 GeV. Assuming a total integrated luminosity of 20 fb −1 per year [73], this would correspond to 1.5 − 3 years of data taking. For scenarios in which the DM channel contributes significantly to the invisible decay rate the situation worsens dramatically, with the required luminosity being pushed to values above 200 fb −1 (BR inv = 0.6) or even close to 1 ab −1 (BR inv = 0.7). It is important to emphasize that these are preliminary results and that the developement of a dedicated search strategy with a thorough treatment of the backgrounds is needed to obtain conclusive results. We leave the implementation of such search strategy for future work.

Conclusions
Guided by the current null findings of DM signals, beyond its gravitational influence, which exert strong constraints on WIMP DM models, we embarked on the construction of a DM theory in which second generation leptons and the DM are suitably charged under a new spontaneously broken Abelian U (1) µ gauge group such that the WIMP paradigm is attainable and at the same time all current direct detection, indirect detection and collider experimental constraints are satisfied. In fact, by construction due to the axial nature of the interaction between the DM and the muon, the theory predicts such small contributions to spin independent DM-nuclei scattering cross sections that they are buried under what is known as the neutrino floor. The strongest current constraints on the model are indirectly related to DM and, in those regions of parameter space where the relic density does not overclose the Universe, they come from contributions beyond the SM to the muon anomalous magnetic moment (g − 2) µ and to neutrino trident production.
We studied all these experimental constraints, focusing in particular on the current searches done at the LHC involving leptons and missing energy and showed that there remain large regions of parameter space that evade all of them and for which the DM content can satisfy the latest DM relic density measurements. We also showed that parts of these unconstrained regions could be probed by future collider measurements at the √ s = 14 TeV LHC for luminosities of 300 fb −1 and 3000 fb −1 (even reaching a discovery level in some regions in this case) with a muon-specific search strategy using an invariant mass window around the Z mass in the more sensitive 3µ+E miss T channel, though some unconstrained regions would remain leading still to an elusive WIMP.
Another very interesting feature of the model is that due to the new abelian charges of the SM second generation leptons, crucial to suppress the spin-independent direct detection cross sections, the model is gauge anomalous and can only be interpreted as a low energy effective theory of a non-anomalous UV theory in which part of the fermion spectrum responsible for anomaly cancellation has been integrated out. An implication of this anomalous nature are triple gauge couplings in the low energy effective theory between the anomalous Z and the EW gauge bosons of the SM, which due to their loop-nature tend to be suppressed, but if able to be probed lead to a window into the UV physics responsible for the anomaly cancellation. We showed that attempts to probe these anomalous couplings at the LHC even at very high luminosities L = 3000 fb −1 are extremely challenging. However, should a muon collider be built in the future, we have demonstrated that due to the large on-shell Z production cross section it would be feasible with relatively low luminosities to probe the anomalous couplings, in particular in the µ + µ − → Z → ZZ resonant search. Specifically, for this signal we showed that in scenarios where the invisible branching ratio of the Z is dominated by the decay into neutrinos, luminosities in the range (30 − 60) fb −1 would be enough to detect the anomalous decay channel, provided that the Z mass is between 250 and 500 GeV.

A Expressions of loop integrals appearing in anomalous triple gauge couplings
Here we provide general expressions of loop integrals necessary to compute the decay rates of Z into Zγ and ZZ [78]:

B Example of cutflows for the LHC analysis
By way of illustration, we provide some cutflows in Tables 2 and 3 for the dedicated search of the 3µ + E miss T signal described in Section 3.6 and whose results are presented in Section 4. The total integrated luminosity is set at 300 fb −1 . Selection cuts include points 1 and 2 of the search strategy in Section 3.6.  Table 2: Cutflow of signal events generated for a particular point of BMII, with m Z = 350 GeV and g µ = 0.46. The number of generated events is rescaled here to a total integrated luminosity of 300 fb −1 . The selection cuts are constituted by: p µ T > 50 GeV, |η µ | < 2.4, p j T > 20 GeV, |η j | < 2.4, ∆R jµ > 0.4, N b = 0, N µ = 3 with net charge = ±1. The m T > 110 GeV cut is applied to the muon that is not included in the pair which better reconstructs the Z mass. For more details regarding the cut definitions see Section 3.6.  Table 3: Cutflow of signal events generated for a particular point of BMIII, with m Z = 200 GeV and g µ = 0.35. The number of generated events is rescaled here to a total integrated luminosity of 300 fb −1 . The selection cuts are constituted by: p µ T > 50 GeV, |η µ | < 2.4, p j T > 20 GeV, |η j | < 2.4, ∆R jµ > 0.4, N b = 0, N µ = 3 with net charge = ±1. The m T > 110 GeV cut is applied to the muon that is not included in the pair which better reconstructs the Z mass. For more details regarding the cut definitions see Section 3.6.