Collider Detection of Dark Matter Electromagnetic Anapole Moments

Dark matter that interacts with the Standard Model by exchanging photons through higher multipole interactions occurs in a wide range of both strongly as well as weakly coupled hidden sector models. We study the collider detection prospects of these candidates, with a focus on Majorana dark matter that couples through the anapole moment. The study is conducted at the effective field theory level with the mono-$Z$ signature incorporating varying levels of systematic uncertainties at the high-luminosity LHC. The projected collider reach on the anapole moment is then compared to the reach coming from direct detection experiments like LZ. Finally, the analysis is applied to a weakly coupled completion with leptophilic dark matter.


Introduction
One of the defining features of dark matter (DM) models is the nature of the interaction between the DM candidate and Standard Model (SM) particles. This interaction is usually assumed to be mediated by heavy states (for example the SM Z or Higgs boson [1,2], or Z bosons [3,4] belonging to extensions of the SM), although in recent times there has been a surge of interest in hidden-sector models in which the interaction is mediated by new dark forces [5].
A class of models with a long history that lies somewhat between these two options is one in which DM interacts with the SM by exchanging photons through higher multipole interactions [6] - [16]. At dimension 5 and 6, the effective operators for multipole interactions of a DM fermion χ can be written as follows Here, d M , d E , and A denote the magnetic, electric, and anapole moment, respectively. For Majorana DM, only the anapole operator is non-zero and can be written in effective field theory (EFT) as where Λ is the cutoff scale.
These effective descriptions have been studied in a variety of UV completions -for example, models of technicolor [8,12], composite DM [13], supersymmetry [6] and, recently, simplified models of leptophilic DM [17], [18]. In weakly coupled completions, DM is assumed to couple at tree level to heavy charged states and hence at one loop to the photon. Restricted to the case of supersymmetry, this could be a model of Bino DM coupling to sleptons [19]. In strongly coupled completions, the DM candidate can be a composite state of charged particles. The cutoff scale Λ corresponds to either the mass of heavy charged states running in the loop, or the scale of confinement in the strongly coupled hidden sector.
The purpose of this paper is to explore the detection prospects at the HL-LHC of DM with couplings to the SM shown in Eq. (1.1). We focus, in particular, on the anapole operator of Eq. (1.2) and calculate the reach of the high luminosity LHC (HL-LHC) in probing the cutoff scale Λ for different DM masses. A collider study of electric and magnetic dipole moments is left for future work 1 . Our collider study is conducted in the conservative and comparatively clean mono-Z channel, incorporating varying levels of systematic uncertainties. We utilize analysis methods developed recently by a subset of the authors. After carefully choosing kinematic variables that can discriminate between signal and SM background in Section 2, we select cuts using the Bayesian optimization method implemented in the Python algorithm Hyperopt [72]. A Boosted Decision Tree (BDT) is then used to classify events into signal and background classes, after a joint optimization of kinematic cuts and BDT hyperparameters.
We have two main motivations for this study. The first is to connect to the rather extensive body of literature on direct detection prospects for this class of DM models. Electromagnetic anapole and dipole DM has been investigated in the context of experiments like DAMA, CDMS, XENON and LUX by several groups [22] - [30], and projections for the future LZ experiment based on a simple scaling of the scattering cross section has been given in [17], [18]. At the level of the EFT, the reach in the cutoff scale Λ obtained from our collider study can be mapped on to a reach in the value of the anapole moment A through Eq. (1.2) (taking dimensionless couplings g ∼ O(1)). Since the anapole moment determines the scattering cross section of DM off nuclei, the collider reach can then be compared to the reach of DM direct detection experiments. We do this in Fig. 6, comparing the values of the anapole moment probed by LUX and LZ to the values probed by our collider study, assuming 5, 10, and 20 % systematic uncertainties at the HL-LHC with 3000 fb −1 of data.
A second motivation for our paper lies in applying this comparative study to a particular UV completion. Over the last few years, there has been steady interest in models of leptophilic DM interacting with the SM via heavy charged mediators [31] - [36]. A comprehensive one-loop analysis of the direct detection phenomenology of this class of models has been performed by [35]. The relic density and indirect detection rates have been calculated by [35] and their embedding within supersymmetry has been studied by [31]. In [37], constraints on this class of models were obtained under the assumption of a DM spike 1 We note that collider studies of magnetic dipole DM have been performed in [20] and [21], for both the LHC and the ILC, using techniques different from the ones employed in our paper.
near the supermassive black hole at the center of our Galaxy. For small mass separation between the charged mediators and the DM candidate, these models are difficult to probe at the LHC by direct production of the mediators themselves [38], [39]. On the other hand, an explicit calculation reveals that the anapole moment is enhanced precisely in these compressed regions of parameter space. Since the coupling to the photon increases, we obtain a corresponding enhancement in the performance of our collider study, as well as the scattering cross section with nuclei. The expectation, then, is that our collider study, in conjunction with projections from LZ, should be able to probe these compressed parameter regions. This aspect of our study is conducted in Section 6.
The paper is structured as follows. After performing our collider analysis in Section 2 -Section 4, the results of our EFT analysis are displayed in Fig. 5 in terms of the cutoff scale Λ, and in Fig. 6 in terms of the anapole moment A. In the latter figure, the limits from current and future direct detection experiments are also displayed, following a discussion of the methods used to calculate those limits in Section 5. The EFT results are then applied to the case of a simplified model with charged mediators in Section 6. We end with our Conclusions.

Collider Study
In this section, we first provide a brief overview of prior work on multipole DM. We then present the results of our collider study.
The relic density of anapole and dipole DM has been worked out by many authors [40] - [45]. In particular, [41] calculated the relic density in the DM mass range m χ ∼ O(100−500) GeV, incorporating annihilation channels like χχ → W + W − and χχ → tt. For DM with m χ ∼ 100 GeV, the correct thermal relic density is obtained for a value of the cutoff scale Λ ∼ 700 GeV. While noting, from our Fig. 5, that this critical value of Λ will be probed at the HL-LHC, we will in general remain agnostic about the relic density constraint. Depending on the cosmological history of the Universe prior to Big Bang Nucleosynthesis, a wide range of annihilation cross sections can in any case be allowed [46], [47].
Low mass multipole DM has been studied by several groups in the context of anomalies in direct detection experiments [22] - [30]. We note that low mass multipole DM has also found applications in addressing the longstanding Solar Abundance Problem (discrepanies between solar spectroscopy and helioseismology) [48]. The preferred anapole moment in such models turns out to be Λ ∼ O(1) GeV −2 . However, in our current paper, we restrict our attention to DM with mass ∼ O(100) GeV. Collider searches for low mass anapole or dipole DM will require different techniques that are left for future investigation.
We now turn to a discussion of the collider prospects of anapole DM with an emphasis on mono-X channels, where X = j, γ, Z, h [49] - [57]. These channels can serve as fertile places to search for anapole DM at colliders. The only requirement is that initial or final states contain charged particles which can emit a photon, which will then ultimately couple to the anapole DM. To our knowledge, the earliest appraisal of anapole DM in the context of the LHC appeared in [41]. The authors performed a monojet study with 19.5 (10.5) fb −1 of CMS (ATLAS) data at 8 TeV and obtained bounds on the cutoff scale Λ ∼ 350 GeV. We note that the monojet cross section is expected to dominate over the mono-Z signature for this class of models. Moreover, the mono-Z channel suffers from the usual branching ratio penalty of demanding Z decay to leptons. On the other hand, though, we also expect that systematic uncertainties on the background should scale more favorably for the mono-Z process at the HL-LHC. A comparative study of mono-Z and monojet signatures for this class of models at the HL-LHC is left for future work. For now, our collider results based on the mono-Z signature should be treated as a conservative estimate.
The mono-Z channel, moreover, offers a good compromise between the signal production cross section on the one hand, and the information available in the Z boson decay on the other. Our strategy will be to use the final state particle distributions to train a decision tree algorithm in order to efficiently separate signal and background events. While allowing Z boson decays to jets would increase the number of signal events, the SM background associated with jets plus missing E T is large. On the other hand, the leptonic decay mode is a viable alternative since it is a cleaner channel and the associated cross section is not too much smaller. Previous searches for dark matter in the mono-Z channel with and without machine learning tools showed good discovery prospects [4,55].
We therefore perform a mono-Z search in the leptonic channel at the 13 TeV LHC. Our signal is where = µ, e come from the Z boson and the dark matter pair from the virtual photon. The backgrounds considered in this work are the main irreducible ones • ZZ(γ * ) → + − + ν ν , and • W + W − → + − + ν ν , and the main reducible ones • ZW → ± ∓ ± + ν and The irreducible τ + τ − background is very small after τ decays to lighter leptons. The single top background W t has a final state similar to the tt background when the W boson and the top quark decay leptonicaly but with a somewhat smaller jet multiplicity. Yet, just like tt, as we are going to see in the next section, the larger jet multiplicity makes the W t identification by the BDT classifier very efficient. Because the W t rate is an order of magnitude smaller that tt, it can be safely neglected.
We require the following basic selection criteria for the mono-Z events: two opposite-charge leptons (electrons or muons) with transverse momentum greater than 20 GeV in the central region of the calorimeter and missing energy larger than 20 GeV for trigger purposes. These initial cuts are loose since ultimately we will tune the E T cut concurrently with tuning the hyperparameters of the machine learning (ML) algorithm. This approach proved to be very efficient in optimizing the performance of the decision trees algorithm in SM double Higgs production [58], a project undertaken recently by some of the authors. We also do not demand an explicit jet veto at this point. Instead, we pass the number of reconstructed jets and leptons to the decision trees algorithm in order to facilitate the identification of the reducible backgrounds, as discussed in the next section. The lack of a lepton invariant mass cut may be noticeable as well. The task of selecting the events by cutting on kinematic variables is the job of the BDT and may be delegated entirely to the ML training phase. We actually found that the BDT performs better when we keep the pre-selection of events at a minimum.
The DM effective operator in Eq. (1.2) was implemented in FeynRules [59]. Events were generated with MadGraph5 [60] at leading order with one extra jet. The hard and soft jet regimes were matched in the MLM scheme [61] at appropriate separation scales. Hadronization was performed with Pythia6 [62]. For the detector simulation and jet clustering we used Delphes3.3 [63] and FastJet [64] with the anti-kt algorithm. The luminosity was fixed at 3 ab −1 projecting the results to the end of the experiment.
The matched cross sections of signal and background processes after the basic selection criteria are displayed in table 2 for DM masses from 100 to 500 GeV and Λ = 1 TeV at the 13 TeV LHC.

Kinematic Variables for BDT Discrimination
In this Section, we describe the kinematic variables used to represent our simulated data. Each event is represented by a real-valued vector composed of the following ten kinematic variables, inspired by the mono-Z study performed in ref. [4]: • Missing energy E T . This variable is used both for cutting and BDT training. Events with heavier dark matter are characterized by harder E T spectrum.
where ∆φ is the angle between the two dimensional vector E miss T and the transverse momentum p Z T of the Z boson candidate. This variable is a measure of axial-E T , which is the projection of E miss T in the direction opposite to the Z candidate [65]. It is useful to differentiate among various DM operators in EFT [55].
• The variable | E T − p Z T |/p Z T , which is the fractional p T difference [65].
• The variable ∆φ( + , − ), which is the azimuthal separation of the two leptons.
• The variable α T = E T ( 2 )/M T , where E T 2 is the transverse energy of the softest lepton of the + − pair and • The variable cos(θ * ) [67] where θ * is defined as the boost invariant cos(θ * ) = tanh . This variable has been used in supersymmetric studies, and is correlated to the production angle of sparticles. It has been studied in decay processes such as slepton to leptons + E T [67], or sbottoms to bottom jets + E T [68].
• The variable M T c = 2 ( p T · p T + p T p T ), which is the contransverse mass [69].
• The variable M , which is the invariant mass of the leptons pairs. It is useful in distinguishing between leptons from Z and W decays.
• The variable n , which is the number of leptons identified in the event. The majority of W Z events present three charged leptons.
• The variable n j , which is the number of jets. Leptonic top quark pair production events present at least two hard jets.
Some of the variables defined above are shown in Fig. 1. In the left upper plot we show the missing E T distribution for the backgrounds (solid lines) and for the lighter (100 GeV) and the heavier (500 GeV) dark matter considered in this work. As expected, the heavier the DM particle is, the harder is its spectrum compared to the softer backgrounds. In the right upper plot, we display the difference of the azimuthal angle of the two hardest leptons in the event. As a consequence of the harder missing transverse momentum, DM events present more collimated leptons from the Z boson decay compared to the backgrounds. Two good variables to discern the reducible backgrounds ZW and tt are shown in the lower plots. On the left, we can see that the typical transverse momentum of the lepton pair is balanced by p T in events containing Z bosons, but less balanced for events containing W s. Finally, on the right is displayed the number of hard jets identified in the event. Events with top quarks contain at least two hard jets and a considerably large fraction of events with higher jet multiplicity compared to the other processes. As observed in ref. [4], the other variables not shown in Fig. 1 present good discrimination power as well, especially the lepton invariant mass. In contrast to ref. [4], however, we use all these variables to represent our simulated data in the learning process of the BDT classifier.

BDT Classification and Performance
In this Section, we discuss our decision tree analysis and give our results.
We utilized the XGBoost package [70] to train boosted decision trees. Approximately 1.5 million events were generated, with around 300,000 for each event class -i.e., the signal class and the four background classes. One hundred thousand events were singled out for optimization purposes and the remaining to train/validation and test the BDT in the proportion of 2/3 and 1/3 of the events, respectively. Since the DM mass is not known at the stage when the mono-Z signal is studied, it is not possible, in principle, to optimize the ML algorithm to discern the signal for a given mass. Yet, the DM mass in the signal class needs to be fixed to train the BDT. We chose to fix the DM mass at 100 GeV, the initial value of our mass scan. This choice was motivated by the fact that a 100 GeV particle signal is harder to discern from the background than a heavier one, as can be observed from the distributions shown in Fig. 1. Our expectation is that the algorithm also present a good performance for the heavier (and easier) dark matter signals. We checked that, in fact, heavier dark matter is more easily identified as a signal event in this approach. A more sophisticated approach, based on parameterized neural networks, is also possible [71].
The BDT hyperparameters, the E T threshold cut, and the number of leptons and jets in order a given event be vetoed were all adjusted jointly in a Bayesian optimization framework with HyperOpt [72]. The joint optimization of cuts and ML hyperparameters is advantageous once the kinematic cuts affect the performance of these algorithms in way that is hard to predict. Reducing the number of background events with hard cuts help  Figure 2. BDT performance plot: Normalized confusion matrix for the multi-channel analysis at m χ = 100 GeV, and systematic uncertainties of 5% for the test samples. The grid consists of the signal class (labeled "dm"), and the four major background classes (labeled ZZ, ZW , W W , and tt) corresponding to the backgrounds discussed below Eq. (2.1). It is evident that around 22% of the signal events are predicted to be ZZ events and 30% of ZZ events are predicted to be signal events.
to increase the signal significance, but has a deleterious effect on the BDT classification as the kinematic distributions of the various classes become more similar to each other. Delegating all the discrimination to the ML side, on the hand, might not suffice if the backgrounds are too much larger than the signal. The optimum point in this trade off is achieved by the joint optimization as described in details in ref. [58].
The signal significance was calculated taking a systematic uncertainty in the total backgrounds yield, ε B . We estimated four scenarios: one with almost no systematic uncertainties, taking ε B =1%, and several others with varying degrees of uncertainties, taking ε B = 5, 10 and 20%. The joint optimization described in the previous paragraph was performed taking these systematic uncertainties into account, in order that the optimization algorithm learns to increase the S/B ratio and tame the effects of these uncertainties as explained in detail in ref. [58].
We now show the performance plots for our BDT analysis. We will discuss, in turn, the confusion matrix, the Receiver Operating Characteristic (ROC) curves, the BDT output distributions for signal and backgrounds, and, finally, the Λ scale that can be probed with 3 ab −1 of data at the 13 TeV LHC, all obtained for the test samples.  . The tt has a larger background rejection rate due its larger jet multiplicity.
In Fig. 2, we show the normalized confusion matrix for the multi-channel analysis at m χ = 100 GeV, and systematic uncertainties of 5%. The grid consists of the signal class (labeled "dm"), and the four major background classes (labeled ZZ, ZW , W W , and tt) corresponding to the backgrounds discussed below Eq. (2.1). As expected, it is the ZZ background that is most difficult to discern from the signal. From the confusion matrix we see that around 22% of the signal events are predicted to be ZZ events and 30% of ZZ events are predicted to be signal events. The reducible ZW and tt are easily identified as background events due their larger lepton and jet multiplicities, respectively, as anticipated in the previous section.
In Fig. 3, we show the Receiver Operating Characteristic (ROC) curves for the multichannel analysis at m χ = 100 GeV, and systematic uncertainties of 5%. The red, green, blue, and yellow curves show the ZZ, W W , ZW , and tt backgrounds, respectively. The results of the confusion matrix are corroborated by the ROC curves of Fig. 3. Indeed, it is the ZZ background that shows the smallest area under curve (AUC), signifying the smallest background rejection for a given signal acceptance. We note, on the other hand, that tt has a larger background rejection rate, mainly because of its larger jet multiplicity. The lepton multiplicity helps to discern the ZW backgrounds which has the second largest AUC. We found that the optimum E T cut varies between 50 to 90 GeV depending on the systematics level and no jet or lepton vetoes. That is, the best performance was achieved by delegating the task of enhancing the classification performance more to the BDT, and less to the kinematical cuts. . BDT performance plot: The cutoff energy scale Λ for discovery (5σ) as a function of BDT score cut for a 100 GeV dark matter and assuming a 5% systematics. We chose a cut of 0.95 in the BDT output to separate signal and backgrounds.
In Fig. 4, we plot the cutoff energy scale Λ for discovery (5σ) as a function of BDT score cut. We see that the Λ scale which can be probed by the 13 TeV LHC increases very rapidly as the output cut approaches 1. We chose to keep the cut at 0.95 for the sake of the stability of the results. A harder cut probes regions with too few background events, which leads to larger fluctuations. In order to estimate the reach in Λ, we averaged the results of ten runs, randomly reshuffling the train and test samples at each run. The uncertainty of the LHC reach in Λ is around ±10-15% of the estimated Λ depending on the mass at the 0.95 output threshold. The results shown in Fig. 4 are from this averaging process, keeping a 5% systematics level. Results for other DM masses and systematic errors were obtained in the same way. After the BDT output cut, 31 (18) signal and 136(164) background events survive for a 100(500) GeV dark matter.
In Fig. 5, we display the 5σ reach in the cutoff energy scale Λ as a function of DM mass for several values of systematics uncertainties. The luminosity is fixed at 3000fb −1 . The performance of the BDT classifier improves as the dark matter gets heavier. This behavior can be understood when we look at the distributions of Fig. 1. For example, the 500 GeV DM presents kinematic distributions which are less similar to the backgrounds. Most importantly, the signal distribution in that case is distinct from the ZZ background, making the BDT classification more efficient. Of course, as the cross section decreases with the DM mass, the Λ scale that can be actually probed drops with the DM mass as shown in Fig. 5. In the case of the 300 GeV mass, the reach of the scale Λ only changed slightly compared to the 200 GeV case -a persistent effect up to the 10% systematics level. In this case, the gain in the BDT classification is competitive with the drop in cross section. For larger masses, however, the number of signal events produced was not enough to beat the better classification achieved with the ML algorithm and the estimated LHC sensitivity in Λ is degraded.
Assuming a very small systematic uncertainty of 1%, the scales for which the anapole DM can be discovered, at 5σ, are all above Λ ∼ 1.1 TeV, as can be seen in Fig. (5). Assuming systematic uncertainties at the level of 5%, a 100 GeV DM can still be discovered if Λ ≈ 1 TeV. From 5% to 10% and from 10% to 20% systematics, the discovery reach in Λ decreases approximately by 200 GeV for a given DM mass.

Comparison with Direct Detection Limits
In this Section, we compare the collider constraints on the anapole moment to constraints coming from direct detection experiments.
We are interested in DM masses O(100 GeV) and typical nuclear recoil energy ∼ 10-30 keV. This corresponds to a DM-nucleus momentum transfer where q q q is the three-vector of the four-momentum q. We assume that the DM has a coupling with the electromagnetic field given by Eq. (1.2). At small momentum transfer, the interaction of the DM and the nucleus can be described by the following effective Lagrangian where J µ is the nuclear current operator. We are neglecting the q 2 dependence in A; a detailed calculation performed in a previous paper by a subset of the authors revealed that this introduces at most a 0.6% error in the anapole moment for the energies considered here [18]. The differential cross section for the scattering of DM with nuclei is given by [17,22,30,73] Here, the mass of the target nucleus is denoted by m T , E R = q q q 2 /(2m T ) denotes the recoil energy of the nucleus, Z the nuclear charge, and v the velocity of the DM particle. F Z is the nuclear charge form factor, while the nuclear spin form factor is denoted by F s . The second term corresponds to scattering off the nuclear dipole moment d A , which is small for Xenon. The differential rate per unit target mass is where, ρ 0 is the local DM density. The minimal velocity of DM that is required for a recoil energy E R is given by v min = m T E R /2/M red , where M red is the reduced mass of the DM-nucleon system. The DM velocity distribution in the rest frame of the detector is given by f ⊕ ( v). Based on these calculations, several groups have calculated the scattering cross section of multipole DM. The constraints can be depicted as upper bounds on the anapole or dipole moments, and we will mainly use the results obtained in the m χ ∼ O(100) GeV range. The constraint on the scattering cross section for a given DM mass obtained from LUX 2016 [74] can be scaled to the corresponding projected constraint at LZ. We will take the most optimistic scenario, with one background event in 1000 days of exposure of 5.6 tonne fiducial mass [75]. The exclusion limit on the scattering cross section for this stringent projection of LZ is expected to be lowered by a factor of ∼ 7 × 10 −4 compared to the LUX 2013 results [76]. Clearly, this is a rough estimate and a careful analysis of future datasets will be needed.
In this context, we note that LZ projections for anapole DM have also been performed in [45] 2 . Our LZ limits are approximately a factor of 2-3 more stringent. For example, for DM with mass ∼ O(100) GeV, we obtain A/µ N ∼ 4 × 10 −7 fm, where µ N is the nuclear magneton. On the other hand, [45] adopt the projection limit of A/µ N ∼ 1 × 10 −6 fm. Clearly, using the LZ projection of [45] would only increase the relative importance of our HL-LHC study vis-vis direct detection.
The resulting comparative study between the collider reach and the reach from direct detection experiments is shown in Fig. 6. On the horizontal axis, we plot the DM mass in Figure 6. A comparison between the bounds coming from LUX 2016 [74] and projected LZ [75] limits on the anapole moment, versus those coming from our HL-LHC study with varying levels of systematic uncertainties. the range 100-500 GeV, while on the vertical axis we plot the anapole moment in units of the nuclear magneton. The LUX 2016 results are shown in the brown dash-dot-dot curve, while the black solid curve shows the projected LZ limit. The green dot-dashed, blue dotted, and red dashed curves show the results obtained from our collider study, with 20%, 10%, and 5% systematic uncertainties and 3000 fb −1 of data at the HL-LHC, respectively.
We can see that there is a vast improvement, almost by an order of magntiude, in the most optimistic projected LZ limit compared to the current LUX limit. This corresponds to the fact that the constraint on the scattering cross section is expected to become stronger by ∼ 100, although, as we have stressed, this is a rough estimate. On the other hand, the HL-LHC is expected to constrain the anapole moment by a factor of 2 − 6 compared to the current LUX results, depending on the level of systematics. It is possible that a study in other channels such as monojet, or a combination of channels, will strengthen the collider results.

Application to a Simplified Model
In this Section, we finally apply our EFT analysis to a specific simplified model. We choose a weakly coupled UV completion in which the DM is a Majorana fermion χ that couples to an uncolored fermion f (with mass m f ) and a pair of charged scalars f L,R . At one loop, the DM couples to the photon through a anapole moment interaction. The mass of the charged scalars will be taken to be ∼ 250 GeV, while the results taken from the EFT will correspond to a cutoff scale Λ ∼ 800 GeV. While the EFT results can be trusted to a first approximation, we note that a detailed collider study of the simplified model will yield more precise constraints.
The Lagrangian of the model is given by A nonzero mixing angle α is allowed between the scalar mass and chiral eigenstates The two scalar mass eigenvalues are denoted by m f 1 and m f 2 . The free parameters of the model are the four masses m χ , m f 1 , m f 2 and m f . A supersymmetric embedding of this model has been studied in [31,32]. Here, we briefly summarize the dependence of the anapole moment on the model parameters. For a full derivation, we refer to the Appendix of [18].
The relevant Feynman diagrams consist of a triangle loop with either two fermions f or two scalars f , and external legs given by two DM particles and a photon. Let us take the momentum of the incoming and outgoing DM particles to be given by p and p , respectively. The total off-shell scattering amplitude is given by where the momentum transfer is denoted by q = p − p and the anapole moment by A(q 2 ). The anapole moment A(q 2 ) can be expressed as A(q 2 ) = e |λ L | 2 cos 2 α − |λ R | 2 sin 2 α X 1 (q 2 ) + e |λ L | 2 sin 2 α − |λ R | 2 cos 2 α X 2 (q 2 ) , (6.4) where X 1,2 is the result of three-point loop integrals. The derivation of these equations, along with the full form of X i , is given in [18]. In the limit |q 2 | m 2 f and |q 2 | m 2 f i , the X i reduce to a simple expression, . This limit applies to DM direct detection for f = µ, τ . For very heavy mediators µ i 1, X i vanishes as µ −1 i log µ i . If the mass difference between f and the DM is small, the value of X i will increase; in the limit (µ i − 1) ∼ δ 1, . The bluish grey region shows the part of parameter space that will be constrained by mono-Z searches at the HL-LHC, assuming 3000 fb −1 of data and 5% systematic uncertainties. The red region is the region that will be constrained by LZ, assuming the most optimistic performance with one background event in 1000 days of exposure of 5.6 tonne fiducial mass. Current LUX limits and direct collider searches for the mediators f 1 are not able to constrain this part of parameter space.
For f = µ and τ , this model has a sizable anapole moment, which can be detected in direct detection experiments. We now present the limits from various experiments shown in Fig. 6 in the parameter space of this class of models.
In Fig. 7, we first plot the constraints on the (λ, m χ ) plane. On the vertical axis, we plot λ R = 2λ L , while on the horizontal axis we plot m χ in the range 100-200 GeV. We keep the mixing angle fixed at α = π/4. The mass of the lightest scalar mediator f 1 is kept at In the region of parameter space plotted, the only constraints come from our HL-LHC study and LZ projections. Current LUX constraints on the anapole moment are too weak to show up in this region, while direct search constraints for the uncolored mediator f 1 lie below 100 GeV. The bluish grey region shows the part of parameter space that will be constrained by mono-Z searches at the HL-LHC, assuming 3000 fb −1 of data and 5% systematic uncertainties. The red region is the region that will be constrained by LZ, assuming the most optimistic performance with one background event in 1000 days of exposure of 5.6 tonne fiducial mass. We note that lower values of µ 1 , corresponding to a greater degree of compression . The bluish grey region shows the part of parameter space that will be constrained by mono-Z searches at the HL-LHC, assuming 3000 fb −1 of data and 5% systematic uncertainties. The red region is the region that will be constrained by LZ, assuming the most optimistic performance with one background event in 1000 days of exposure of 5.6 tonne fiducial mass. Current LUX limits and direct collider searches for the mediators f 1 are not able to constrain this part of parameter space. between the lightest scalar mediator and the DM, have larger values of the anapole moment from Eq. (6.5). These regions are constrained both by our HL-LHC study and by LZ projections, indeed to a greater extent than is the case in Fig. 7. On the other hand, larger values of µ 1 are less suited to our collider search strategy as well as DM direct detection, due to a reduced value of the anapole moment. These regions start becoming constrained by collider searches for the mediators f themselves, as the mass gap between them and the DM increases leading to collider signals with large missing energy.
We show the range of DM masses between 100-200 GeV. Below m χ = 100 GeV, there are LEP constraints on the uncolored scalar mediators. Above 200 GeV, the anapole moment becomes smaller and the constraints become less stringent.
Finally, changing the value of α also changes the anapole moment. We turn to this dependence next.
In Fig. 8, we display the constraints on the (λ, α) plane. We keep the mass of DM and the lightest scalar mediator fixed at m χ = 200 GeV and µ 1 = m 2 f 1 /m 2 χ = 1.44, respectively. The color scheme is the same as in the previous figure. We see the presence of "blind regions" -regions near α = π/8, 7π/8 -where the anapole moment becomes highly suppressed, from Eq. (6.4). These regions are difficult to probe using any method that Figure 9. Constraints on the (λ R , λ L ) plane: We show constraints on the λ R versus λ L plane, keeping m χ = 200 GeV, µ 1 = m 2 f1 /m 2 χ = 1.44, and α = π/4. The bluish grey region shows the part of parameter space that will be constrained by mono-Z searches at the HL-LHC, assuming 3000 fb −1 of data and 5% systematic uncertainties. The red region is the region that will be constrained by LZ, assuming the most optimistic performance with one background event in 1000 days of exposure of 5.6 tonne fiducial mass. Current LUX limits and direct collider searches for the mediators f 1 are not able to constrain this part of parameter space. relies on the photon coupling; in [18], however, it was shown that indirect detection can constrain these regions. The effect of changing either µ 1 (and hence the light mediator mass) or the mass of the DM has been discussed before, and applies to this figure as well. Lowering µ 1 leads to larger values of the anapole moment and stronger constraints on the (λ, α) plane, and the "blind regions" get sharpened to values very close to α = π/8, 7π/8. Increasing m χ weakens the collider and direct detection constraints.
Finally, in Fig. 9, we plot our results in the (λ R , λ L ) plane, keeping m χ = 200 GeV, α = π/4, and µ 1 = 1.44. The color scheme is the same as in the previous figures. The corridor around the region where λ R ∼ λ L constitutes a "blind region" where the anapole moment is attenuated for α = π/4. These regions are difficult to probe for methods that rely on the photon coupling. Parts of this region can be explored by indirect detection.

Conclusions
In this paper, we have explored the HL-LHC detection prospects for DM that couples to the Standard Model through higher electromagnetic moments, particularly the anapole moment. The study is conducted at the level of EFT, with the aim of calculating the reach in the cutoff scale Λ. We have conducted our study in the mono-Z channel, taking into account varying levels of systematic uncertainties. We carefully choose kinematic variables that can discriminate between signal and SM background, and select cuts using the Bayesian optimization method implemented in the Python algorithm Hyperopt . A BDT is then used to classify events into signal and background classes.
The results of our collider study are shown in Fig. 5. We see that for a very small systematic uncertainty of 1%, the 5σ reach in Λ is above 1.1 TeV for DM masses in the range 100-500 GeV. Assuming a larger systematic uncertainty of 5% , a 100 GeV DM can still be discovered for Λ ≈ 1 TeV. From 5% to 10% and from 10% to 20% systematics, the discovery reach in Λ decreases by 200 GeV for a given DM mass.
The 5σ reach in the cutoff scale Λ obtained from our collider study can be mapped on to a reach in the value of the anapole moment A through Eq. (1.2). The resulting comparative study between the collider reach and the reach from direct detection experiments is shown in Fig. 6. The LUX 2016 results are shown in the brown dash-dot-dot curve, while the black solid curve shows the projected LZ limit. The green dot-dashed, blue dotted, and red dashed curves show the results obtained from our collider study, with 20%, 10%, and 5% systematic uncertainties and 3000 fb −1 of data at the HL-LHC, respectively.
Finally, the EFT analysis is applied to a specific simplified model. We choose a weakly coupled UV completion in which the DM is a Majorana fermion χ that couples to an uncolored fermion f (with mass m f ) and a pair of charged scalars f L,R . At one loop, the DM couples to the photon through a anapole moment interaction. The HL-LHC constraints on the parameter space of this class of models is presnted in Fig. 7 -Fig. 9. These constraints are juxtaposed with projected LZ constraints on the parameter space, assuming the most optimistic performance with one background event in 1000 days of exposure of 5.6 tonne fiducial mass.