Hunting for Hypercharge Anapole Dark Matter in All Spin Scenarios

We conduct a combined analysis to investigate dark matter (DM) with hypercharge anapole moments, focusing on scenarios where Majorana DM particles with spin 1/2, 1, 3/2, and 2 interact exclusively with Standard Model particles through U(1)$_{Y}$ hypercharge anapole terms for the first time. For completeness, we construct general effective U(1) gauge-invariant three-point vertices. These enable the generation of hypercharge gauge-invariant interaction vertices for both a virtual photon $\gamma$ and a virtual $Z$ boson with two identical massive Majorana particles of any non-zero spin $s$, after the spontaneous breaking of electroweak gauge symmetry. For complementarity, we adopt effective operators tailored to each dark matter spin allowing crossing symmetry. We calculate the relic abundance, analyze current constraints and future sensitivities from dark matter direct detection and collider experiments, and apply the conceptual naive perturbativity bound. Our estimations based on a generalized vertex calculation demonstrate that the scenario with a higher-spin DM is more stringently constrained than a lower-spin DM, primarily due to the reduced annihilation cross-section and/or the enhanced rate of LHC mono-jet events. As a remarkable outcome, the spin-2 anapole DM scenario is almost entirely excluded, while the high-luminosity LHC exhibits high sensitivities in probing spin-1 and 3/2 scenarios, except for a tiny parameter range of DM mass around 1 TeV. A significant portion of the remaining parameter space in the spin-1/2 DM scenario can be explored through upcoming Xenon experiments, with more than 20 ton-year exposure equivalent to approximately 5 years of running the XENONnT experiment.


I. INTRODUCTION
The nature of dark matter (DM) remains one of the greatest puzzles in particle physics and cosmology, accounting for approximately a quarter of the total energy content in the Universe.Exploring the non-gravitational interactions of DM with Standard Model (SM) particles offers a direct path to uncovering new structures and symmetries.This exploration can be pursued through three distinct experimental approaches: direct detection, indirect detection, and collider experiments.These three different types of experiments can be complementary to one another, and hence it is extremely helpful to combine the results and apply them to a single effective operator for DM-SM interactions.This is a reasonable and widely applicable approach in scenarios where the corresponding effective field theory (EFT) is valid [1][2][3][4].
Certainly, for a proper EFT description, it is required that the energy scale of all processes under consideration be well below the masses of the mediating particles, and that their interactions respect the established low-energy (global and/or gauge) symmetries of the SM [5][6][7][8].As the TeV energy scale, which is higher than the electroweak (EW) scale, is currently being probed and the SM with its SU(3) C × SU(2) L × U(1) Y gauge symmetry has been firmly established, particularly with the discovery of the Higgs boson [9,10], it is appropriate to maintain both the hypercharge U(1) Y symmetry and the other non-Abelian gauge symmetries, SU(3) C and SU(2) L .For instance, the importance of considering the hypercharge U(1) Y rather than the electromagnetic U(1) EM as a valid U(1) gauge symmetry has been clearly and persuasively demonstrated from various physics perspectives in recent work [11].
The non-gravitational interaction between DM and SM particles through an electromagnetic (EM) form factor, for DM with nonzero spin, is realized through a higher-dimensional operator.Consequently, this provides an interesting test bed for the EFT approach in DM studies.Additionally, the EM form factor induces unexpected (and suppressed) EM interactions of DM with SM particles, without the DM being directly charged.The phenomenological effects of such interactions were first discussed in Ref. [12].
Previously, several aspects of EM anapole DM interaction terms have been studied in the context of relic abundance measurements as well as in searches at direct detection, indirect detection, and collider experiments for the spin-1/2 case [62][63][64][65][66][67][68][69][70] and the spin-1 case [71,72].A comprehensive analysis of the U(1) Y anapole DM scenario with a spin-1/2 Majorana particle was recently conducted [11].On the contrary, the analysis for the spin-1 case was restricted so far to the EM anapole DM scenario, as the focus was only on the direct DM detection capability.In this context, it is worthwhile to perform a generic analysis accommodating and characterizing the hypercharge anapole DM particle of any nonzero spin with combined experimental probes all together.
In our general analysis, we construct effective U(1) Y gauge-invariant three-point vertices of a hypercharge gauge boson B and two identical on-shell Majorana particles with any non-zero spin. 1 These vertices accommodate an arbitrary spin s and non-zero mass m, generating interaction vertices with not only a virtual photon γ but also a virtual gauge boson Z after electroweak symmetry breaking (EWSB).A concise overview of the distinction between our analysis targets and those of other studies is shown in Table I of our summary and conclusion section.
Having outlined the general three-point vertices, our in-depth numerical analysis is specifically targeted towards four scenarios where the DM particle spin is set to be 1/2, 1, 3/2, and 2, while qualitatively exploring the implications for the scenarios with its spin larger than 2. We incorporate the relic abundance value determined by the Planck collaboration [74], the up-to-date result from the DM direct detection experiment XENONnT with approximately 1.1 ton-year exposure [75], and the LHC experiments with the integrated luminosity of 139 fb −1 [76][77][78], as well as the so-called naive perturbativity bound (NPB) making our EFT approach valid. 2 In addition, we estimate the projected sensitivities of the high luminosity LHC (HL-LHC) experiment with the full run of 3 ab −1 integrated luminosity [79] and those of the future XENONnT with the 20 ton-year exposure.This paper is organized as follows.In Sec.II, we derive the general effective hypercharge gaugeinvariant three-point vertices.These generate interaction vertices for a virtual photon γ and a massive 1 For clarity, the term ' Majorana particle' is conventionally used to represent a spin-1/2 self-conjugate particle; however, we will generalize it to include self-conjugate particles of any spin without loss of generality, as done in the work [73]. 2 Conceptually, there could also be unitarity bounds on the couplings for each annihilation mode, but these are quantitatively much weaker than the naive perturbativity limits.
The derivation is based on an efficient and systematic algorithm for constructing the covariant effective vertex for three particles of any spin and mass [80][81][82].In Sec.III, we calculate the annihilation cross sections of two Majorana particles into kinematically allowed pairs of SM particles.We then determine the constraints on the effective coupling strengths for the spin-1/2, 1, 3/2, and 2 cases from the observed DM relic abundance.Section IV is devoted to determining constraints on the coupling strengths from recent LHC and upcoming HL-LHC experiments [76,79].We particularly focus on how these constraints depend on the spin of the anapole DM Majorana particle.In Sec.V, we explore constraints from the up-to-date results of the DM direct detection experiment XENONnT [75] and the projected sensitivities from its highly enhanced exposure, and study the implications for the coupling strengths.Based on all the experimental constraints and an additional theoretical constraint from the naive perturbativity bound, we present an overall combined picture of the current constraints and sensitivities in Sec.VI.This comprehensive analysis places special emphasis on systematically delineating the distinct characteristics that vary depending on the spin values.We summarize the key points of our results and conclude in Sec.VII.Furthermore, Appendix A provides a concise and systematic algorithmic description for constructing all the general three-point anapole vertices efficiently and Appendix B provides a detailed explanation of the numerical calculation strategy employed to determine the DM relic abundance of dark matter.

II. ANAPOLE VERTICES FOR TWO MAJORANA PARTICLES OF ANY SPIN
In this section, we present an efficient and systematic algorithm for constructing covariant threepoint anapole vertices for two identical Majorana particles of any spin.This algorithm allows us to perform a complete and systematic characterization of all the spin values of the hypercharge anapole DM.Subsequently, we delve deeper into the specific case of covariant vertices, focusing on the spin-1/2, 1, 3/2, and 2 scenarios.This detailed analysis encompasses both analytic and numerical investigations to provide comprehensive insights and findings.

A. Aanapole three-point vertices of Majorana particles of any spin
A Majorana particle, by definition, is a CPT self-conjugate particle, which means it remains unchanged under the combined operations of charge conjugation (C), parity reversal (P), and time reversal (T).These particles do not possess any static charge or multipole moments, as all the terms in their interaction Hamiltonian are CPT-odd.
The only permissible U(1) gauge-invariant interaction vertices between a U(1) gauge boson and two identical massive Majorana particles of nonzero spin are known as anapole-type moments [73].
It is worth noting that no massless CPT self-conjugate Majorana particle can have any U(1) gaugeinvariant couplings unless its spin is 1/2.
The effective anapole three-point χχB vertex of two identical massive Majorana particles χ of any spin and a U(1) gauge boson B is in general given by the gauge invariant Lagrangian with the 4-vector current J µ comprised of two Majorana-particle fields and the field-strength tensor half-integer or integer spin-s Majorana particle.The curved arrow is for an arbitrary chosen fermion-number flow direction whose meaning is described in detail in Refs.[83,84].Equation (2.1) enables us to construct the effective conserved current for the annihilation of two Majorana particles into a virtual vector boson depicted in Fig. 1.Explicitly, the current can be cast into the form with p = k 1 + k 2 and q = k 1 − k 2 in terms of the incoming Majorana-particle momenta, k 1 and k 2 , where the vector current J µ is nothing but the momentum-representation version of the positionrepresentation current J µ in Eq. (2.1).Note that the vector current V µ automatically satisfies the U(1) gauge invariance condition p µ V µ = 0.In the covariant formulation, the vector current V µ can be written as the products between the wave tensors χ s (k 1 ) and χ s (k 2 ) of two Majorana particles and a covariant three-point vertex Γ with respect to the arbitrarily chosen fermion-number-flow arrow shown in Fig. 1, where the indices, α and β, stand collectively for the 4-vector indices, n = s for a half-integer or integer spin-s Majorana particle.All the details about the wave tensors [85][86][87][88][89][90] and the general covariant vertex are included collectively in Appendix A.
As two identical Majorana particles annihilate into a virtual gauge boson B * , the covariant threepoint vertex must satisfy the so-called identical-particle (IP) condition, (A23) or (A24), in the fermionic or bosonic case, respectively, as described in detail in Appendix A. By employing the general covariant three-point vertices and imposing the IP relation, the anapole three-point vertex can be cast into a compact square-bracket operator form: in terms of 2s independent f and b couplings for the spin-s Majorana fermion with s = n + 1/2 and boson with s = n, respectively.The cutoff scale Λ is set forth explicitly to indicate that the effective covariant vertex originates from a higher-dimensional operator term in the given effective Lagrangian or Hamiltonian.Here, the square-bracket notations are introduced for denoting the product of the basic helicity-related operators as well as two derived operators in a compact form with all the fourvector and spinor index symbols hidden.
Firstly, the two square-bracket operators, [A] and [V ], in Eq. (2.4) denote an orthogonal axial-vector and vector currents, of which the explicit forms are given by with an orthogonal gamma matrix γ ⊥µ = g ⊥µν γ ν and an orthogonal metric tensor g ⊥µν = g µν − pµ pν + qµ qν involving two normalized momenta, p = p/ p 2 and q = q/ −q 2 .Secondly, due to the totally-symmetric property of the wave tensors [85][86][87][88][89][90] over all the four-vector indices, the n-th power products of the metric tensor g and the basic scalar operator S 0 can be given in a compact square-bracket form: where the basic scalar operator S 0 is defined by of which the repeated appearance at the vertices increases the dimensions of the corresponding Lagrangians gradually.It is compensated by introducing the proper power of the cutoff scale Λ along with the operator as shown in Eqs.(2.4) and (2.5).Thirdly, the other two derived basic vector operators in terms of the normalized momentum p and the basic scalar operators S ± of which the explicit expression is with the angle-bracket notation ⟨αβ pq⟩ = ε αβρσ pρ qσ of a product between an anti-symmetric Levi-Civita tensor and two normalized momenta p and q.
B. Effective three-point anapole vertices for the spin of 1/2, 1, 3/2 and 2 Following the systematic derivation procedure for constructing the general anapole vertices and using the general properties of wave tensors described in Appendix A, we can recast the covariant three-point anapole vertices extracted from the general forms in Eqs.(2.4) and (2.5) effectively into the following form as in terms of a single coupling a 1/2 in the spin-1/2 case and two independent couplings, a 1 and b 1 , in the spin-1 case with an orthogonal antisymmetric tensor ⟨αβµq⟩ ⊥ = g ⊥µν ϵ αβ νσ q σ , and in terms of a single coupling a 3/2 in the spin-3/2 case and two independent couplings, a 2 and b 2 , in the spin-2 case up to the leading order in 1/Λ 2 .All the couplings are in general complex and the a i terms are parity-odd while the b j terms are parity-even where i = 1/2, 1, 3/2, 2 and j = 1, 2. Although our numerical analysis is confined up to spin 2, a simple extrapolation suggests that a single coupling exists in any half-integer spin scenario, while in the case of non-zero integer spin, there are two distinct and independent couplings at the leading order of 1/Λ 2 , as evident with Eqs.(2.4) and (2.5).
The effective U(1) gauge-invariant spin-1/2 and spin-1 anapole Lagrangians corresponding to the vertices in Eqs.(2.13) and (2.14) can be re-constructed by replacing each momentum with its corresponding derivative as in terms of the spin-1/2 and spin-1 Majorana fields, χ 1 2 and χ 1µ , respectively.Likewise, the U(1) gauge-invariant spin-3/2 and spin-2 anapole Lagrangians corresponding to the vertices in Eqs.(2.15) and (2.16) are given as: ) in terms of the spin-3/2 and spin-2 Majorana fields, χ 3 2 µ and χ 2µν , respectively.In the following, we specify the U(1) gauge boson to be the hypercharge U(1) Y gauge boson B in the SM.The hypercharge gauge field B µ is decomposed into a photon field A µ and a Z-boson field cos θ W and s W = sin θ W of the weak mixing angle θ W after the firmly-established EWSB.Furthermore, for simplicity and without loss of generality, we omit the spin index s of the Majorana particle χ s in the following discussion, as it applies universally across all spin cases.

III. DM RELIC ABUNDANCE
In this section, we calculate the relic abundance of our anapole DM of each spin from the thermal freeze-out mechanism.The overabundant region beyond the observed relic abundance [74] is simply considered to be excluded without introducing late time reduction possibilities.The corresponding areas are shown in the two-dimensional planes of mass of dark matter and the couplings a i , b j normalized by the cutoff scale squared Λ 2 .The self-conjugate hypercharge anapole DM particles can annihilate into the SM particles via the s-channel photon γ and Z boson exchanges.If kinematically allowed, the DM particles annihilate mainly via the processes χχ → f f , W − W + and/or ZH, where f is a SM quark q = u, d, s, c, b, t or lepton ℓ = e, µ, τ, ν e , ν µ , ν τ , as depicted in Fig. 2. The relic abundance of χ can be determined through the freeze-out of the annihilation processes in the figure.As noted previously in Ref. [73], every annihilation cross section is completely factored into a simple product of two independent parts, of which one corresponds to the DM annihilation into a virtual gauge boson and the other to the sequential decay of the virtual gauge boson into a pair of SM particles. 3Explicitly, the total cross 3 The angular distribution of each annihilation mode is uniquely determined independently of the DM particle spin.
This characteristic spin-independent angular distribution was demonstrated explicitly in the scattering process e − e + → γ * → χχ via a photon exchange in Ref. [73].
section of each annihilation mode can be written in the following compact form: ) for the spin-1/2 and spin-1 Majorana particles, respectively, and for the spin-3/2 and spin-2 Majorana particles, respectively, with the collision energy √ s and the speed of χ, β χ = 1 − 4m 2 χ /s in the center of mass (CM) frame of two χ particles.Here, the dimensionless SM pair-production term P SM = f P f f + P W W + P ZH is the sum of all the kinematically allowed production terms: with Here, for the sake of notation, the normalized propagator factor Π Z (s) and an effective SM vector coupling squared in terms of the SM vector and axial-vector couplings, , of the Z boson to a SM fermion pair f f with the isospin component I 3 f and electric charge Q f of the fermion f .We emphasize again that the SM production terms are independent of the DM particle spin and its couplings to the photon and Z boson.Consequently, the information on the characteristics of the anapole DM particle is encoded exclusively in the DM annihilation into a virtual gauge boson.
By calculating the DM relic abundance using the methodology described in detail in Appendix B, we can derive exclusion limits on the effective anapole couplings versus the DM mass.These limits

IV. LHC SEARCHES
In this section, we derive the exclusion limits on the couplings for the hypercharge anapole DM particle from its searches at the LHC experiment with the projected sensitivities on the couplings from the upcoming HL-LHC experiment [77,78,91,92].Although there are possibly various production channels at the LHC, we consider the most dominant production processes for the hypercharge anapole DM particle in our analytic analysis.
It was shown in a previous work [69] that the EM anapole DM particles in the U(1) EM gaugeinvariant framework can be produced dominantly through the di-jet processes via vector-boson fusion, especially when investigating them with strong experimental cuts.On the contrary, the di-jet processes cannot be dominant anymore in the hypercharge anapole DM case because not only the γ exchange diagram but also Z-boson exchange diagram contribute to the process, leading to a quite significant cancellation in the high-energy regime so that the unitarity problem is diminished extremely efficiently [11].As a result, the strongest LHC constraints on the hypercharge anapole DM couplings are expected to come from the so-called mono-jet processes pp → j + X with X standing for the collection of invisible particles including two DM particles, χχ.
Two parton-level processes, gq → qχχ and q q → gχχ, contribute dominantly to the mono-jet process pp → j + X at the LHC.Quantitatively, the former gq cross section is much larger than the latter q q cross section.In this light, we consider the process gq → qχχ in the present work as the most crucial mode for investigating the constraints from the LHC mono-jet events on the couplings versus the DM mass m χ in the spin-1/2, 1, 3/2, and 2 cases. 4There are two Feynman diagrams contributing to the sequential process gq → q γ * /Z * → qχχ, of which one is a s-channel quark-exchange mode and the other a t-channel quark-exchange mode as depicted in Fig. 4. As demonstrated in all the DM annihilation processes in Sect.III, the effective anapole three-point γχχ and Zχχ vertices allow each q q q q q q γ * ⊕ Z * γ * ⊕ Z * χ χ χ χ g g parton-level cross section to be completely factored into the SM 2-to-2 scattering process gq → qγ * /Z * and the DM pair production processes through the 2-body decay γ * /Z * → χχ.After performing the 2-body phase space integration over the invisible final two-body χχ system, we can obtain the partonlevel cross-section for the process gq → qχχ in a compact integral form as for the spin-1/2, 1, 3/2, and 2 cases, respectively, with the χχ invariant mass Q 2 corresponding to the virtual γ * or Z * invariant mass, the speed factor with the strong-interaction coupling g s , the normalized propagator Π Z (Q 2 ) and the modified vector coupling Vq of the quark q defined in Eqs.(3.8) and (3.9),where the parton-level q transverse momentum, which is invariant under the Lorentz boost along the beam direction, is pT = √ ŝ 2 (1 − Q 2 /ŝ) sin θ with the polar-angle θ between the momenta of the initial gluon g and the final quark q.The function F(ŝ, Q 2 ; p T ) with the transverse-momentum cut p T in Eq. (4.5) is given explicitly by with the maximal transverse momentum The transverse-momentum cut p T is introduced to regularize the forward singularity caused by neglecting the quark mass; this also allows us to ignore particles escaping detection along the proton-beam pipe directions.Compared to the spin-1/2 case, the parton-level production cross section in the spin-1 case has an additional kinematic enhancement factor Q 2 /4m 2 χ originating from the longitudinal mode of one of the two spin-1 DM particles as in Eq. (4.2).We note in passing that the power of the enhancement factor gets bigger for higher spin cases due to the larger number of longitudinal modes so that the mono-jet searches at the LHC impose much stronger constraints on the couplings gradually.Certainly, the parton-level cross sections should be folded with proper quark and gluon parton distribution functions for evaluating the mono-jet cross section at the LHC.
In this paper, we adopt a simple statistical analysis for deriving the LHC and HL-LHC limits on the hypercharge anapole couplings versus the DM mass, focusing on the effective characterization of the hypercharge anapole DM particle according to the spin.Let us calculate the simplest version of the signal significance z defined by with the numbers s and b of the signal and background events. 5We set the critical significance value z = 2 to determine the 95% confidence level (CL) exclusion limits on the couplings normalized to the cutoff scale squared Λ 2 .The events are selected according to the selection criteria p jet T > 250 GeV and |η| < 2.5, applied to the most recent mono-jet searches with an integrated luminosity of 139 fb −1 at the LHC energy of 13 TeV by the ATLAS collaboration in Ref. [77]. Figure 5 show the current exclusion limits from the LHC (shaded region bounded by the solid lines) and the expected sensitivities at the HL-LHC with the full running of the 3 ab −1 integrated luminosity (dashed lines).Comparing the panels from the top left one, it is clearly seen that the couplings in the higher-spin cases are constrained more strongly in the whole kinematically-available region due to the higher power of the enhancement factor Q 2 /4m 2 χ , especially in the low mass region as indicated in Eqs.(4.2), (4.3) and (4.4).
Close to the kinematical endpoint of the mass m χ ∼ 6 TeV, there is a slight reduction in the a 1 /Λ 2 (a 2 /Λ 2 ) constraint compared to the b 1 /Λ 2 (b 2 /Λ 2 ) constraint, due to the presence of the higher power of the kinematical suppression factor The dashed line in each panel shows the future sensitivity of the planned HL-LHC experiment with a slightly larger collision energy of 14 TeV and an integrated luminosity of 3 ab −1 , roughly 10 times larger than the present LHC luminosity [76,79].We require slightly stronger selection criteria of the mono-jet events for the HL-LHC: p jet T > 300 GeV and |η| < 2.5.As for the number of background events at √ s = 14 TeV, we consider the difference in the cross sections of the dominant background process pp → Zj → ν ν j, enhanced by 11.57/9.88≈ 1.17 compared to the √ s = 13 TeV case [91].
The projected sensitivities become stronger as m χ decreases.Due to the Q 2 dependence of the gq → qχχ production cross section, we expect that the 100 TeV future circular collider experiment under research and development (R&D) [93,94] enables us to cover a much larger region of the couplings versus the anapole DM mass.
Before closing this section, we emphasize once more that the power of the invariant mass square in the mono-jet cross section increases in proportion to the spin value of the anapole DM particle.The enhancement arises from the increased number of longitudinal modes of the higher-spin anapole DM particle.It strongly indicates that the constraints on the couplings versus the DM mass from the LHC and HL-LHC mono-jet searches become much stronger as the spin value of the anapole DM particle increases.

V. DM DIRECT DETECTION
In this section, we describe how the exclusion limits on the hypercharge anapole DM couplings versus the DM mass are extracted from the recent DM direct detection experiment XENONnT with the 1.1 ton-year exposure [75] which is at present the most powerful DM direct detection experiment.
Along with the exclusion limits, we consider the projected sensitivities of the future XENONnT with the 20 ton-year exposure.Non-relativistic dark matter is assumed to move with a typical velocity of order v ≃ 10 −3 c in the galactic halo. 6Thus, the recoil energy of a DM particle against a heavy target nucleus is expected to be in the keV energy scale, much smaller than the typical DM mass ≳ O(GeV) in consideration here as well as the Z-boson mass m Z = 91.2GeV.Hence, the Z-boson exchange contribution to the elastic DM scattering off the target nucleus can be safely ignored because the momentum transfer is significantly smaller than the Z boson mass m Z and the DM scattering off the nucleus is dominated by the photon exchange [11].
Taking the small recoil energy limit and ignoring the Z-boson exchange contribution safely we can cast the recoil-energy dependent differential cross section into the following factorized form: ) ) for the spin-1/2, 1, 3/2, and 2 cases, respectively, with the target nucleus mass m T and the recoil energy E R .Note that the momentum transfer q = √ 2m T E R is much smaller than our dark matter mass in consideration.One noteworthy feature is that the recoil-energy dependent function A(E R ) is factored out independently of the spin of the anapole DM particle: with the atomic number Z T of the target nucleus and the DM particle speed v χ relative to the nucleus.
The charge form factor F C is given by in terms of the first-kind spherical Bessel function j 1 of order 1 where r C = (c 2 + 7π 2 a 2 /3 − 5s 2 ) 1/2 with c ≃ (1.23A 1/3 − 0.60) fm, a ≃ 0.52 fm and s ≃ 0.9 fm, and the atomic mass A of the target nucleus, while the magnetic dipole moment form factor F M [99,100] is given by with the radius r M = 1.0A 1/3 fm.The recoil-energy dependent function A in Eq. (5.5) involves the nuclear magneton µ p = e/2m p with the proton mass m p and the weighted dipole moment μT for the target nuclei [101]: where f i , µ i , and s i are the abundance fraction, magnetic moment, and spin of the isotope i.
The recoil-energy dependent distribution of the DM direct detection process is given by integrating the differential cross section of each DM spin as in Eqs.(5.1) -(5.4) over the DM velocity with the distribution f LAB (⃗ v ) in the laboratory frame as where we use the local DM density ρ loc = 0.3 GeV cm −3 .In the present work, we adopt a simple Maxwell-Boltzmann distribution in the galactic frame truncated at the escape speed v esc of our Galaxy: with the velocity ⃗ v E of the Earth in the galactic frame and where the normalization constant N is By integrating out the recoil-energy distribution dR/dE R in Eq. (5.9) and taking into account the detection efficiencies of the XENONnT experiment [75], we can evaluate the number of recoil DM detection events.For the details of calculating the expected signal events we refer to Appendix B of Ref. [71].The top (bottom) left panel of Fig. 6 shows the 90% CL constraint on the normalized  from the different spin-averaged and polarization-weighted factors, 1/6 and 7/120, in the detection rates, respectively.For the expected sensitivities of the upcoming XENONnT with the 20 ton-year exposure shown with the dashed lines, we relied simply on scaled statistics without accounting for potential future improvements in background rejection and the control of systematic uncertainties.
As shown previously in Fig. 5, the LHC and HL-LHC mono-jet constraints on the couplings of a lower-spin particle are much weaker than those on the couplings of a higher-spin case, especially for the DM mass m χ ≲ 1 TeV.Hence, the LHC/HL-LHC and DM direct detection experiments can play quite complementary roles in imposing the exclusion limits on the couplings versus the DM mass.

VI. COMBINED CONSTRAINTS AND FUTURE SENSITIVITIES
In this section, we show the results combining all the aforementioned experimental constraints and future sensitivities on the effective couplings versus the DM mass coming from the Planck determination of the DM relic abundance, the LHC and HL-LHC mono-jet searches, and the present DM direct detection experiment XENONnT and its future data of 20 ton•yr exposure, which have been evaluated systematically in the previous three sections.On top of those experimental bounds and sensitivities, we include an additional theoretical constraint from the naive perturbativity bound (NPB) for guaranteeing the validity of the EFT formalism, which needs to be taken with a grain of salt. 7The energy-dependent NPBs on the couplings simply read for the spin-1/2, 1, 3/2, and 2 cases, respectively.As the CM energy √ s ≥ 2m χ , the NPB condition (6.1) applied to the asymptotically high-energy limit leads to the following inequality relations: on the effective couplings versus the DM mass m χ .The imposition of those constraints is a very loose statement on the tree-level perturbativity.If the limits are violated, we naively expect that higher-loop corrections must be included, which is beyond the scope of this paper.(2) case, while setting the other coupling to zero.
In the spin-1/2 case, two regions are still allowed.One tiny allowed region is near the Z pole with even though the permitted area is substantially reduced compared to the spin-1/2 scenario, due to the much more strengthed constraints from the relic abundance and LHC experiment.While the complete HL-LHC operational phase is not anticipated to entirely cover this small permissible region, the future 100 TeV circular pp collider [93] promises the capability to thoroughly investigate and potentially close off the remaining segments of this area, because of its vastly greater collision energy and significantly improved sensitivity.
Despite the slightly weaker XENONnT constraints on the coupling |a 3/2 |/Λ 2 in the spin-3/2 case than the smaller spin cases due to the spin-averaged and polarization-weighted factors, the currently allowed parameter region is smaller.This is because the LHC constraint gets much stronger stemming from the enhancement of the mono-jet production cross section by the larger number of longitudinal modes as explained in Sec.IV.
The full running of the HL-LHC is expected to probe nearly half or more of the very tiny allowed regions.This effect leads to a remarkable result for the higher spin case, i.

VII. SUMMARY AND CONCLUSION
Our investigation focused on a scenario where DM is characterized as a Majorana particle possessing a non-zero spin, interacting solely with SM particles via hypercharge anapole terms.This scenario renders us to pursue the combined analysis from various experimental results in a complementary way via the EFT approach.For the experimental/observational constraints and sensitivities, we applied the Planck measurement of the DM relic abundance, the direct detection experiment XENONnT, and the LHC mono-jet searches together with the expected sensitivities at the HL-LHC and the future XENONnT.Because of the straightforward calculations and strong theoretical foundations, we conducted a focused numerical analysis on DM spins s = 1/2, 1, 3/2, and 2, making comparative analyses among them for the first time within the realm of anapole DM studies.A succinct summary of the anapole Dark Matter scenarios and experimental searches analyzed in this paper is presented in Tab.I, juxtaposed with previous literature for comparison.The expectation for the higher spin DM scenarios will be briefly discussed at the end of this section.Considering the theoretically allowed range of the EFT approach together with a grain of salt, We demonstrate that the hypercharge anapole DM scenarios are currently on the verge of being discovered or ruled out.The main results of the present work can be summarized with the following key points: • As easily expected, the relic abundance imposes stronger constraints on d-wave terms than the • The LHC and HL-LHC mono-jet searches play the most crucial role in probing the higher-spin anapole DM lighter than about 1 TeV.This is due to the kinematic factor associated with the longitudinal modes of the DM particle.The mono-jet searches play an essential role, especially for probing the parity-even terms of the spin 1 or 2 DM, i.e., |b 1,2 |/Λ 2 , since their contributions to the DM direct detection cross sections are suppressed by the recoil energies.
• The momentum transfer in the DM direct detection process is sufficiently small and hence the dominant contribution is through the t-channel photon exchange.
It is noteworthy that the result from the near future XENONnT with the 20 ton-year exposure is expected to explore the allowed region of space that coincides with the HL-LHC expectations for the spin-1/2 case, as shown in the top left panel of Fig. 7. Remarkably, the upcoming XENONnT experiment has the potential to achieve the extended coverage within approximately 5 years, which is significantly sooner than the projected full operational timeline of the HL-LHC.Thus, the future DM direct detection experiments, such as the planned XLZD consortium involving multiple Xenon target experiments with more than 20 ton-year exposure, could eventually probe beyond the allowed region bounded by the HL-LHC expectation at a faster pace, in the spin-1/2 case.
• The NPB condition is equally and approximately applied to all the spin-1/2, 1, 3/2, and 2 couplings.The lower NPB bound is inversely proportional to the square of the DM mass, as presented in Eq. (6.2).It is important to note that a breach of this conceptual NPB suggests a breakdown of the EFT framework, potentially resolvable by a more fundamental UV theory [13][14][15].Hence, the NPB bounds need to be accepted cautiously.
Overall, the combined analysis shows that the hypercharge anapole coupling of a higher spin DM is more stringently constrained or expected to be probed sooner than that of a lower spin DM as evident in Fig. 7.We expect the 100 TeV proton-proton circular collider experiment (FCC), which is currently under R&D would have a potential to completely probe the whole remaining parameter space of the anapole DM scenarios considered here.Note that our analysis result combining the constraints from the observed relic abundance and the mono-jet searches at the LHC can be extended to a lighter DM mass down to O(10 MeV) as long as the DM relic is determined by the freeze-out mechanism, providing a powerful exclusion bound already.
The DM particle with its spin larger than 2 is regarded to be innately a composite particle in order to avoid various conceptual problems such as unitarity issues.Nevertheless, if the scale of compositeness is significantly high, we expect that scenarios with spins greater than 2 could potentially be completely ruled out because of the substantial kinematic factor linked to an increased number of longitudinal modes of the DM particle in the anapole vertices.although no conclusive comments can be made yet before dedicated studies. of phase transitions, the entropy of the Universe in a co-moving frame, S ≡ sa 3 = (2π 2 /45)g s T 3 a 3 , with the degrees of freedom g s contributing to the entropy density s is conserved throughout its evolution with the decreasing temperature T and, accordingly, the scale factor a. The freeze-out occurs during the radiation-dominated epoch, enabling us to write down the energy density ρ = (π 2 /30)g ρ T 4 , with the degrees of freedom g ρ , involved in the Hubble parameter H = (1/a) da/dt = 8πGρ/3 with the gravitational constant G.In this case, the yield Y (x) of DM particles, varying over x = m χ /T (that is, over time and/or temperature), satisfies the so-called evolution equation: where λ ann = (π/45) 1/2 M pl m χ ⟨ σv Møl ⟩ with the so-called Møller velocity v Møl of two annihilating DM particles [106] , and Y = n/s and Y eq = n eq /s including the DM number density n and the thermalequilibrium density n eq in the co-moving frame.Here, the yield Y eq of the DM particles in thermal equilibrium can be written as with the DM spin degrees of freedom g χ and the second-kind modified Bessel function K 2 of order 2.
The simple asymptotic expression at large values of x ≫ 1 is Y eq (x) ≈ 0.145 g χ g s x 3/2 e −x .(B3) In the present work, we calculate the thermally-averaged cross section ⟨σv Møl ⟩ by adopting its explicit form [106] as calculations.

1 pFigure 1 .
Figure 1.A diagram for the annihilation of two identical Majorana particles χ s into an off-shell hypercharge gauge boson B * .The combined momenta p = k 1 + k 2 and q = k 1 − k 2 are constructed by the combinations of incoming momenta k 1,2 of two Majorana particles.The k 1,2 -dependent χ α s and χ β s are the wave functions of two identical Majorana particles with the particle symbol denoted by χ in the main text.The indices, α and β, stand collectively for the 4-vector indices, α = α 1 • • • α n and β = β 1 • • • β n , with n = s − 1/2 or n = s for a

Figure 2 .
Figure 2. Feynman diagrams for the dominant annihilation processes of two identical Majorana particles into a pair of SM particles, χχ → f f (left), W − W + (middle) and ZH (right) where f is a SM quark q = u, d, s, c, b, t or lepton ℓ = e, µ, τ, ν e , ν µ , ν τ .Here, χ denotes the self-annihilating Majorana DM particle.The red open circle in each diagram indicates the effective three-point anapole vertex.The notation γ * ⊕Z * stands for the combined s-channel γ and Z exchanges.

Figure 4 .
Figure 4. Two parton-level Feynman diagrams contributing to the process gq → qχχ, of which the left one is for a s-channel quark exchange and the right one is for a t-channel quark exchange.If the quark mass is ignored, the t-channel diagram has a forward singularity to be regularized.The notation γ * ⊕ Z * stands for the combined s-channel γ and Z exchanges.
well as the gq collision CM energy √ ŝ and the transverse-momentum cut p T of the produced quark, invariant under any Lorentz boost along the proton beam direction.The parton-level 2-to-2 scattering processes gq → qγ * /Z * are encoded fully in the ŝ and Q 2 dependent cross section O. Explicitly the parton-level effective cross section O with the transverse-momentum cut p T is given by

Figure 5 .
Figure 5. Exclusion limits from the mono-jet events at the LHC (solid line) of 13 TeV energy and 139 fb −1 luminosity and projected sensitivities from those at the HL-LHC (dashed line) of 14 TeV energy and 3 ab −1 luminosity, respectively, derived mainly from the parton-level process gq → qχχ, on the effective couplings and the DM mass m χ .The top (bottom) left panel shows the limits on the normalized coupling a 1/2 /Λ 2 (a 3/2 /Λ 2 ) in the spin-1/2 (3/2) case.The top (bottom) middle panel shows the limits on the normalized coupling a 1 /Λ (a 2 /Λ 2 ) and the top (bottom) right panel shows the limits on the normalized coupling b 1 /Λ 2 (b 2 /Λ 2 ) in the spin-1 (2) case.

esc /v 2 0 ,
with the error function erf(z) = 2 z 0 e −t 2 dt/ √ π.The values of the three different speeds are set numerically to the escape speed v esc = 544 km s −1 , the speed of the Sun relative to the DM reference frame v 0 = 220 km s −1 and the Earth speed v E = 232 km s −1 in the galactic frame.
versus the DM mass m χ in the spin-1/2 (3/2) case.The top (bottom) middle panel shows the constraint on the normalized coupling |a 1 |/Λ 2 (|a 2 |/Λ 2 ) versus m χ , and the top (bottom) right panel shows the constraint on the normalized coupling |b 1 |/Λ 2 (|b 2 |/Λ 2 ) versus m χ in the spin-1 (2) case.In each plot, the excluded region of the corresponding coupling versus the DM mass is shown as the red-shaded area bounded by the exclusion limit with a red solid line.Note that each of the second terms proportional to the recoil energy E R on the couplings, |a 1 |/Λ 2 and |a 2 |/Λ 2 in Eqs.(5.2) and (5.4) is much smaller than unity, i.e., m T E R /(2m 2 χ ) ≪ 1, and hence the direct detection bounds on the |a i |/Λ 2 couplings with i = 1/2, 1, 3/2, 2 are dominantly determined by the spin-averaged factors.On the other hand, the direct detection cross sections are suppressed by m T E R /2m 2 χ once they are dominated by the |b i |/Λ 2 term with i = 1, 2, which results in the negligible sensitivities as shown in the right panels.The tiny difference between the limits on the couplings, |b 1 |/Λ 2 and |b 2 |/Λ 2 , arise

Figure 6 .
Figure 6.Exclusion limits from the DM direct detection experiment XENONnT on the effective couplings versus the DM mass m χ .The solid line is the current limit with the 1.1 ton-year exposure and the dashed line is the expected sensitivity with the 20 ton-year exposure.The top (bottom) left panel shows the limit on the normalized coupling |a 1/2 |/Λ 2 (|a 3/2 |/Λ 2 ) versus the DM mass in the spin-1/2 (3/2) case.The top (bottom) middle panel shows the limit on the normalized coupling |a 1 |/Λ 2 (|a 2 |/Λ 2 ) versus the DM mass in the spin-1 (2) case while setting the other coupling to zero.The top (bottom) right panel show the limit on the normalized coupling |b 1 |/Λ 2 (|b 2 |/Λ 2 ) versus the DM mass in the spin-1 (2) case while setting the other coupling to zero.The right panels clearly show that the the normalized couplings, |b 1 |/Λ 2 and |b 2 |/Λ 2 , for the spin-1 and spin-2 DM cases get much weaker constraints than the other cases.

Figure 7
Figure 7 shows the combined exclusion limits on the effective normalized couplings versus the anapole DM mass m χ from the Planck measurement of the DM relic abundance (black solid), the theoretical NPB condition (orange solid), the mono-jet searches at the LHC of 13 TeV with the integrated luminosity of 139 fb −1 (blue solid) and the full running of the HL-LHC of 14 TeV with the 3 ab −1 integrated luminosity (blue dashed), and the present DM direct detection experiment XENONnT (red solid) with the 1.1 ton-year exposure along with the future XENONnT with the 20 ton-year exposure (red dashed).The top (bottom) left panel is the combined constraint on the normalized coupling |a 1/2 |/Λ 2 (|a 3/2 |/Λ 2 ) versus the mass m χ in the spin-1/2 (3/2) case.The top (bottom) middle panel shows the combined limit on the normalized coupling |a 1 |/Λ 2 (|a 2 |/Λ 2 ) and the top (bottom) right panel shows the combined limit on the normalized coupling |b 1 |/Λ 2 (|b 2 |/Λ 2 )versus m χ in the spin-1 (2) case.Note that the future sensitivities of XENONnT can be reached within about 5 years of running or the XLZD consortium of many Xenon target experiment plans[102].

Figure 7 . 2 (
Figure 7. Combined exclusion limits on the effective couplings versus the DM mass m χ from the measured DM relic abundance of the Planck satellite (black solid), the theoretical NPB (orange solid), the LHC (blue solid) and HL-LHC (blue dashed) mono-jet search experiments and the DM direct detection experiment XENONnT current limit with the 1.1 ton-year exposure (red solid) and the future prospect with the 20 ton-year exposure (red dashed).The top (bottom) left panel shows the combined constraint on the normalized coupling |a 1/2 |/Λ 2 (|a 3/2 |/Λ 2 ) versus the DM mass in the spin-1/2 (3/2) case.The top (bottom) middle panel shows the combined constraint on the normalized coupling |a 1 |/Λ 2 (|a 2 |/Λ 2 ) versus the DM mass, and the top (bottom) right panel shows the combined constraint on the normalized coupling |b 1 |/Λ 2 (|b 2 |/Λ 2 ) versus the DM mass in the spin-1 where the relic abundance constraint becomes much weaker due to the sharp Z-resonance effect.This tiny region is expected to be completely excluded by the upcoming DM direct detection experiments and the HL-LHC experiment, as indicated by the red and blue dashed lines.The other allowed region is a triangle-shape area centered near |a 1/2 |/Λ 2 = 10 −6 GeV −2 and m χ = 1 TeV, nearly the half of which can be probed in the near future by XENONnT after the 20 ton•yr exposure.This shows the effectiveness of the complementary DM EFT approach.In the spin-1 case, the Z-resonance region with m χ ≃ m Z /2 is entirely ruled out by the combined constraints arising from both relic abundance measurements and LHC search experiments, illustrated in the upper middle and right frames of Fig.7.This complete exclusion arises from two synergistic effects, independent of the constraints posed by both the DM direct detection experiment and the NPB bound.Compared to the spin-1/2 case, the relic abundance constraint is bolstered by 1.5 times, attributed to a smaller spin-averaged factor than that in the spin-1/2 case.Furthermore, the constraints imposed by the LHC experiments are significantly heightened, primarily due to the longitudinal mode of the spin-1 DM particle, particularly noticeable for m χ less than 1 TeV.Moreover, the spin-1 case with a non-zero |a 1 | but |b 1 | = 0 is expected to face total exclusion, as depicted in the top middle panel of Fig.7.Here, the combined constraints from the relic abundance measurements, the (HL-)LHC and (upgraded) DM direct detection experiments and the theoretical NPB synergistically contribute to this complete exclusion.Conversely, the spin-1 case with a non-zero |b 1 | but |a 1 | = 0 features a triangular-shaped allowed region centered around |b 1 |/Λ 2 = 10 −6 GeV −2 and m χ = 1, TeV, e., spin-2 DM.As shown in the bottom-middle and bottom-right panels, the LHC mono-jet searches so far have provided extremely strong constraints: full exclusion for the |a 2 |/Λ 2 dominated case and nearly exclusion except for the region m χ ∼ 1.2 TeV for the |b 2 |/Λ 2 dominated case.The remaining region is expected to be fully probed at the HL-LHC in particular with more effective search strategies.Due to the gradual tightening of constraints from relic abundance and LHC data for higher-spin cases, it becomes evident that all hypercharge anapole dark matter particles with masses below the GeV scale are excluded, irrespective of their spin values, provided they adhere to the thermal freeze-out scenario.In order to cover a general model set-up for spin-1 and 2 DM including both non-negligible parity odd terms (|a 1,2 |/Λ 2 ) and even terms (|b 1,2 |/Λ 2 ) in the cross sections, we display the combined constraints and sensitivities in the 2-dimensional (|a i | and |b i |) plane after fixing the DM mass m χ = 1.25TeV and the cutoff scale Λ = 2 TeV in Fig.8.Each of the constraints is given by an ellipse due to the cross sections are proportional to the combination of the absolute squares of two couplings in the spin-1 and spin-2 cases.We ignore the constraints on |b i | from XENONnT which are too weak.The left panel is for the spin-1 and the right panel is for the spin-2 DM.Note that the future XENONnT sensitivity with the 20 ton•yr exposure (red dashed vertical line) is comparable to those of the HL-LHC full running for the spin-1 DM, while the latter becomes more powerful for the spin-2 DM.We expect that a higher spin s > 2 DM scenario would suffer from even stronger bounds although further dedicated studies are needed.The right panel of Fig.8shows that the full running of the HL-LHC can probe fully the spin-2 scenario with m χ = 1.25 TeV.

Figure 8 .
Figure 8.The parameter spaces of two couplings, |a 1 | and |b 1 | (left), and two couplings, |a 2 | and |b 2 | (right) for the specific values of the DM particle mass, m χ = 1.25 TeV, and the cutoff scale, Λ = 2 TeV, in the spin-1 and 2 cases, respectively.The red and blue dashed lines indicate the projected sensitivities from the XENONnT with the 20 ton-year exposure and the HL-LHC experiment after the full running with the integrated luminosity of 3 ab −1 , respectively.

p
-wave terms since larger couplings are required to obtain the right relic abundance.The coefficients of the p-wave terms decrease as the spin of DM due to the spin-averaged and polarizationweighted factors in the annihilation cross sections: 1/4 (|a 1/2 |), 1/9 (|b 1 |), 5/72 (|a 3/2 |), and 1/20 (|b 2 |).This induces stronger constraints as the DM spin increases.In each case, the relic abundance constraint in the Z-pole resonance region is about a factor of 10 times weaker than the others.

3 2 χd
B4)    involving the annihilation cross sections in Eqs.(3.1)  or (3.2) with the second-kind modified Bessel function K 1 of order 1.This is because the conventional series expansions of the cross sections over the relative velocity between two annihilating DM particles do not converge when the annihilation energy √ s is close to the Z boson mass due to the sharp Z-boson pole contribution[107].The exponentiallystiff suppression of the modified Bessel functions for large values of x ≫ 1 renders the numerical calculation unreliable.Instead, we use its analytic asymptotic expression directly as⟨σv Møl ⟩ ∼ x ds s 1/4 (s − 4m 2 χ ) exp x 2 −for x ≫ 1.The exclusion limits on the couplings satisfying the observed DM relic density Ω χ h 2 ≈ 0.12[74] can be derived by integrating out the equation (B1) from x = 1 to x = 10 3 with the initial condition Y = Y eq at x = 1.Numerically, the function λ ann dependent on the annihilation cross section varies stiffly with the DM mass m χ in the GeV region, reducing the accuracy of the numerical calculation[108].The numerical accuracy can be enhanced greatly by taking the integration after recasting the evolution equation as ln g s d ln T e W − e (2Weq−W ) , (B6) with a more slowly varying logarithmic function W = ln Y .The relic abundance calculations described so far are obtained semi-analytically, and the results are consistent with the MicrOMEGAs[109]

Table I .
Comparison between our current research and various previous studies on EM anapole DM and hypercharge anapole DM conducted through three main experiments: relic abundance analysis, searches at the LHC, and/or direct detection experiments.The works are referenced with corresponding citation numbers provided in the references.In the case of spin-1 EM anapole DM, each cross mark (✘) denotes that the relic abundance and LHC search aspects have not been explored to date.On the other hand, the red check mark (✔) signifies the comprehensive quantitative investigation of all the clarified experimental aspects in the spin-1.3/2and 2 cases in addition to the spin-1/2 case, undertaken in our current work.