Maximally self-interacting dark matter: models and predictions

We study self-interacting dark matter (SIDM) scenarios, where the $s$-wave self-scattering cross section almost saturates the Unitarity bound. Such self-scattering cross sections are featured by strong velocity dependence in a wide range of velocities. They may be indicated by observations of dark matter halos in a wide range of masses, from Milky Way's dwarf spheroidal galaxies to galaxy clusters. We pin down the model parameters that saturates the Unitarity bound in well-motivated SIDM models: the gauged $L_{\mu} - L_{\tau}$ model and composite asymmetric dark matter model. We discuss implications and predictions of such model parameters for cosmology like the $H_{0}$ tension and dark-matter direct-detection experiments, and particle phenomenology like the beam-dump experiments.


Introduction
parameters. In Section 3, we consider DM with a light mediator based on the gauged L µ − L τ model. We identify the mediator and DM masses that explain the strong velocity dependence of DM self-interaction. Interestingly, we find that the discrepancy in muon anomalous magnetic moment g − 2 and tension in H 0 are mitigated in the very parameter points. The light L µ − L τ gauge boson are also subject to various experimental searches such as non-standard interactions of neutrinos and missing-energy events in colliders. In Section 4, we consider composite asymmetric dark matter (ADM) based on the dark quantum chromodynamics (QCD) and electrodynamics (QED). The dark QCD scale and dark pion masses are constrained so that dark nucleon DM can have the indicated velocitydependent cross section. The binding energies of two nucleons and vector resonances are also predicted, which have important implications for the dark nucleosynthesis in the early Universe and dark spectroscopy measurements in lepton colliders, respectively. The dark photon, which is required to make cosmology viable, is subject to intensive experimental efforts such as beam-dump experiments. The dark proton may be found in DM directdetection experiments. We give concluding remarks in Section 5.

Effective-range theory
The effective-range theory is first developed in the attempt to understand low-energy scatterings of nucleons [45,46]. See Ref. [47] for a review of the effective-range theory in the context of effective-field theory. Here we refer to Ref. [48] that revisits the effective-range theory in the context of SIDM. The 2-body scattering cross section can be decomposed into the partial-wave (multipole ) contributions: Here we ignore a possible quantum interference for the identical particles. 2 . The momentum k = µv rel is given by the reduced mass µ (µ = m/2 for the identical particles with the mass of m) and relative velocity v rel . The analytic properties of the wavefunction determine the low-energy behavior of δ as Here a and r e are called the scattering length and effective range, respectively. 2 One may also choose to use the transfer cross section σ T , which may be a more suitable quantity to parameterize the momentum-transfer effect among DM particles through elastic scattering. For the scattering of identical particles, σ T is written as [49] σ T = dΩ(1 − | cos θ|) dσ dΩ , where σ is the standard cross section, σ = dΩ(dσ/dΩ). Since we focus on the s-wave scattering, this amounts to an additional factor of 1/2, σ T = σ/2.
In the following, we focus on the s-wave ( = 0), which is expected to be dominant at the low-energy. 3 Then, the low-energy cross section is given by σ = 4πa 2 1 + k 2 (a 2 − ar e ) + a 2 r 2 e k 4 /4 .
When |a| |r e |, in the range of 1/|r e | > k > 1/|a|, the cross section saturates the Unitarity bound, A large |a| indicates the existence of the shallow bound state (a > 0) or virtual state (a < 0; bound state with a non-normalizable wavefunction), depending on its sign. The bound (virtual) state energy is given by the positive (negative) imaginary pole k pole : In Fig. 1a, the data points show the preferred values of σ/m as a function of the relative velocity; they are inferred from the observations on central DM densities of MW's dwarf spheroidal galaxies (cyan), field dwarf spheroidal galaxies (red)/low surface brightness spiral galaxies (blue), and galaxy clusters (green). In the left panel, we show σ/m (orange) in the effective-range theory with (a, r e /a) = (−7.5 × 10 3 fm, −38) (solid) and (a, r e /a) = (−7.5 × 10 3 fm, −3.5) (dashed); the former data set has larger |a/r e |, which means that it is closer to the quantum resonance than the latter (and hence more saturate the Unitarity bound). The solid curve exhibits strong velocity dependence and simultaneously explains the preferred σ/m ∼ 0.1 cm 2 /g at v ∼ 1000 km/s and σ/m ∼ 100 cm 2 /g at v ∼ 30 km/s. As one can see, σ/m nearly saturates the Unitarity bound (brown), so that σ/m ∝ 1/v 2 rel at the presented velocity scales. Interestingly, such velocity dependence provides a good fit to the inferred σ/m from MW's dwarf spheroidal galaxies (cyan) at low-velocities. As we take the parameters away from the quantum resonance, σ/m desaturate like the dashed curve, but may still explain the cyan data points. Meanwhile, there can be seen that the cyan data points disagree with the inferred values from the field dwarf galaxies (red). Hereafter we take the former seriously, while being ignorant to the latter. We give further discussion in Appendix A.
It is illustrating to consider the scattering under the Hulthén potential: The Hulthén potential approximates the Yukawa potential,  Figure 1: (a) Velocity dependence of σ/m (orange) in the effective-range theory (left) and the Hulthén potential (right) for given DM mass and the effective-range theory parameters (a, r e ). The Hulthén-potential parameters took in the right panel exhibits the same effective-range theory parameters as the left panel; (α, m φ ) = (7.7 × 10 −3 , 0.1 GeV) (solid) and (α, m φ ) = (7.1 × 10 −4 , 0.01 GeV) (dashed). The brown curve is the Unitarity bound for the DM self-scattering cross section, σ max /m. The data points are the inferred values of σ/m from the observations on field dwarf (red)/LSB (blue) galaxies, and galaxy clusters (green) [7]. The data points in cyan are the values that may explain the anticorrelation between the central DM densities and their orbital pericenter distances [44]. The effective range theory deviates from the Hulthén potential results at high velocity due to the effective-range expansion. (b) s-wave scattering length (left) and effective range (right) in units of m −1 φ as a function of φ for the Hulthén potential. The stars indicate the benchmark parameters took in Fig. 1a, which are near the quantum resonance.
with a proper choice of δ. Here m φ is the mass of a mediator and we take δ = 2ζ(3)m φ with the ζ(z) being the Riemann zeta function. 4 The Hulthén potential enjoys an analytic expression of the phase shift: Here Γ(z) is the gamma function and This phase shift results in the scattering length and effective range of respectively. Here ψ (n) (z) is the polygamma function of order n and γ E = 0.577 . . . is the EulerMascheroni constant. Above we consider an attractive potential, by taking positive α. One can use the same expressions with an analytic continuation to negative We note that while effective-range theory reproduces well the analytic results of the Hulthén potential in the k → 0 limit, it starts to deviate for k −1 m −1 φ . In other words, when the de Brogile wavelength of the incoming particles is shorter than the range of the Hulthén potential. This is expected since the low-momentum effective-range expansion fails at high k. This feature can also be seen by comparing the left and right panel of Fig. 1a; the right panel is the analytic result of Hulthén potential with the parameters that give the same effective-range theory parameters as in the left panel. Fig. 1b shows a (left) and r e (right) as a function of φ in the units of m −1 φ . The resonances appear at φ = n 2 (n = 1, 2, . . . ). In the following, we focus on the first resonance, n = 1, while discussing higher resonances in Appendix B. The points depicted by stars corresponds to the parameters took in Fig. 1a; they are both near the first resonance, nearly saturating the Unitarity bound.
Above we consider the elastic scattering and see that we need a light mediator and induced relatively long-range force to reproduce the strong velocity dependence. The Sommerfeld enhancement for DM annihilation [50][51][52][53][54][55][56][57][58] is controlled by the same potential. If the Sommerfeld enhancement is also almost on the quantum resonance, this would have multiple implications for the cosmology: strong constraints from the indirect-detection experiments if DM annihilation ends up with the electromagnetic energy injection [59]; and the second stage of the DM freeze-out (i.e., re-annihilation) [58,[60][61][62][63]. The enhancement factor in the Hulthén potential is given by [24] for s-wave and p-wave annihilation, respectively. For annihilation, we take δ = ζ(2)m φ . The left panels of Fig. 2 shows the Sommerfeld-enhancement factor as a function of the relative velocity for s-wave (top) and p-wave (bottom). The Hulthén-potential parameters took corresponds to the parameters depicted as stars in Fig. 1b, which are close to the quantum resonance for elastic scattering. For these parameters, the s-wave Sommerfeldenhancement factor is also resonantly enhanced towards low-velocity as ∝ 1/v 2 rel and saturates; as shown in the left panel of Fig. 2a, the parameters are close to the s-wave quantum resonance for annihilation. This may be dangerous in terms of the constraints from observations on the CMB and indirect-detection experiments. The constraints may be evaded if we consider neutrinos as the final DM annihilation product or consider the p-wave annihilation [25,59]. Indeed, the p-wave Sommerfeld-enhancement factors for the parameters are not huge, since they are sufficiently far from the p-wave quantum resonance for DM annihilation (see the right panel of Fig. 2b).
In Section 3, we consider a model with a light mediator based on the gauged L µ − L τ model, where DM annihilation is p-wave and ends up with neutrinos. ADM is another good candidate, since it involves a light mediator to deplete the thermal relic abundance of the symmetric component. ADM may also be expected to be safe from the indirect-detection bounds at the first sight. In Section 4, nevertheless, we see that even a (inevitable induced) tiny DM-anti DM oscillation leads to a significant annihilation in the late Universe.
3 Dark matter with a light mediator: gauged L µ − L τ model In this section, we consider DM with a light mediator based on the gauged L µ − L τ model [25]. In this model, a new gauge boson Z couples to the muon and tau lepton (and their neutrinos) [64,65]. Interestingly, this new force contributes to the muon (g − 2) µ , alleviating a possible tension between the SM prediction and measurement [66][67][68] (See also the latest review on (g − 2) µ in the SM, Ref. [69]). 6 After other constraints are 5 Note that our φ is different from Ref. [58]. This expression coincides with [50] 6 Two new experiments would shed light on the longstanding tension in (g − 2) µ in the near future: E989 experiment at Fermilab [70] and the E34 experiment at J-PARC [71]. The lattice studies of hadronic  Fig. 1a. The s-wave annihilation is near the quantum resonance for the benchmark parameters, as shown in the right panel. The Sommerfeld-enhancement factor grows as ∝ 1/v 2 rel (brown) towards low velocity and saturates. (Right): φ -dependence of the swave Sommerfeld-enhancement factor. We take v = 1, · · · , 10 −3 from bottom to top. The dashed lines show the Sommerfeld-enhancement factors in the Coulomb potential that correspond to φ → ∞. The points depicted by stars correspond to the benchmark parameters took in Fig. 1a for fixed v = 10 −3 . (b) Same as Fig. 2a, but for the p-wave Sommerfeld-enhancement factor. taken into account, the (g − 2) µ discrepancy is ameliorated for the gauge coupling of g 4 × 10 −4 -10 −3 and the mass of the new gauge boson of m Z 8-200 MeV. Ref. [25] extends this model with the a vector-like pair of fermions, N andN . We assume that the L µ −L τ breaking Higgs Φ carries a unit charge and N (N ) carries a (minus) half charge in units of the muon and tau-lepton charges. 7 The resultant Z 2 symmetry, which is unbroken by the Higgs vacuum expectation value (VEV), guarantees the stability of the lightest mass eigenstate N 1 among N andN , namely, N 1 being DM. We consider the pseudo-Dirac dark matter, whose Dirac mass m N is larger than the Majorana mass induced by the VEV of Φ. The masses are given by the Yukawa coupling y > 0: 8 The relevant interactions are Here N 1 and N 2 denote Majorana fermions and ϕ is a real scalar resulting from Φ. The thermal relic abundance of N 1 is predominantly determined by co-annihilation of N 1 N 1 , N 2 N 2 → Z Z , ϕϕ and N 1 N 2 → Z ϕ. The thermally averaged annihilation cross sections are given by the dimensionless temperature x = m N /T : Here we drop terms proportional to the gauge coupling g , because it is subdominant compared to those with Yukawa coupling y. Note that self-annihilation channels are pwave, while the co-annihilation channel is s-wave. The effective annihilation cross section is given by the relative yield r [83]: 9 σv eff = (1 − r) 2 σv 11 + r 2 σv 22 + 2r(1 − r) σv 12 , r = e −∆m/T 1 + e −∆m/T .
We compare the effective cross section at x = 20 with the canonical cross section, (σv) can = 3 × 10 −26 cm 3 /s. We fix the Yukawa coupling in this way. The above expression of the effective cross section is derived under the assumption that N 1 and N 2 are in chemical equilibrium. We check that in parameters of interest, the chemical equilibrium is achieved by contributions to (g − 2) µ may also resolve the tension in the future. The latest study on the hadronic vacuum polarization contribution [72] claims that the discrepancy between the SM prediction and the experimental result disappears. Their result is also in tension with the other result (such as the R-ratio determination), and thus it is expected to scrutinize their result in detail by the other lattice groups. 7 The charge of Φ determines the neutrino mass matrix by the see-saw mechanism [73][74][75][76]. Our choice is a minimal viable one [77] (see also Appendix C) and may explain the baryon asymmetry of the Universe [78] through the leptogenesis [79][80][81][82]. 8 Hereafter we assume the CP symmetry to restrict the Yukawa coupling. For more general discussion, please see Ref. [25]. 9 As we see below, we consider an O(10) MeV ϕ, by taking a O(10 −6 ) quartic coupling. This may delay the phase transition of the L µ − L τ breaking, possibly inducing a small thermal inflation. Please see Ref. [26] for a further discussion.  Parameter regions for σ/m at the velocity scales of MW's dwarf spheroidal galaxies (cyan), field dwarf galaxies (blue), MW-size galaxies (red), and galaxy clusters (green) for given (m Z , g ). The region where m ϕ < m Z (gray) may be disfavored by the BBN observations [25]; lifetime of ϕ is longer than 1 s. The yellow region depicts the parameters close to the first quantum resonance for elastic scattering, i.e., 0.85 < φ < 1.15 where φ = αm N 1 /( 2ζ(3)m ϕ ); this range includes the benchmark parameters shown as stars in Fig. 1b. (Right): The velocity dependence of σ/m (orange) for the benchmark parameter points (depicted in the left panel) near the quantum resonance. The brown curves are the Unitarity bound, σ max /m. The data points are the same as in Fig. 1a.
Again since the gauge coupling g is tiny, the 2-body potential of DM is also dominated by the Yukawa coupling: where α = y 2 /(8π). 10 We approximate the Yukawa potential by the Hulthén potential. We assume the Higgs mass m ϕ to be larger than the gauge boson mass m Z , namely, m ϕ > m Z ; otherwise the Higgs lifetime exceeds 1 s and cause cosmological problems. The left panels of Fig. 3 shows the parameter regions for σ/m = 100-200 cm 2 /g in MW's dwarf spheroidal galaxies (cyan), 0.1-1 cm 2 /g (light blue) and 1-10 cm 2 /g (blue) in field dwarf galaxies, 0.1-1 cm 2 /g in MW-size galaxies (red), and > 0.1 cm 2 /g (green) in galaxy clusters, where we take v rel = 20, 30, 200, and 1000 km/s, respectively. There, we focus on g and m Z that ameliorates the (g − 2) µ discrepancy. We used the analytic results of the Hulthén potential [19] to calculate the elastic scattering cross section of N 1 ; we multiply additional 1/2 factor to Eq. (1) since we consider the elastic scattering among identical N 1 particles. It should be noted that the analytic results of the Hulthén potential is valid only in the resonant regime, i.e., αm N 1 /m ϕ 1 and m N 1 v rel /m ϕ 1; the analytic formula for the cross section in the classical regime, i.e., αm N 1 /m ϕ 1 and m N 1 v rel /m ϕ 1, is given in Ref. [58], and the Born approximation can be used in the Born regime, i.e., αm N 1 /m ϕ 1. However, we remark that using the corresponding formulas in each regimes do not change our conclusion. We indicate the classical and Born regimes (brown) in the left panels of Fig. 3; for the classical regime, we require the condition for the velocity range of interest, v rel 1000 km/s. From Fig. 3a and Fig. 3b, we see that with a light mediator, i.e., m ϕ 20 MeV, σ/m can be as large as ∼ 100 cm 2 /g in MW's satellites, while diminishing towards galaxy clusters. In order to exhibit such large cross section at low-velocities, the parameters should lie very close to the quantum resonance (yellow); the velocity dependence near the quantum resonance, σ/m ∝ 1/v 2 rel , fits well the cyan data points. As an example, we show the velocity dependence of σ/m for the benchmark parameters depicted (as star and diamond) in the right panels of Fig. 3.
With Z heavier than m Z 20 MeV, it is impossible to achieve such large cross sections (σ/m ∼ 100 cm 2 /g) at MW's satellites even at the quantum resonance (if we took the error bars at face values). This is demonstrated in Fig. 3c. There, we take m Z = 50 MeV as an example. The cosmological constraints restricts m ϕ m Z . Thus, with heavier Z , the constraint also restricts N 1 to be heavier near the quantum resonance (yellow). Eventually, it becomes impossible to achieve ∼ 100 cm 2 /g in MW's satellites even if σ/m saturates the Unitarity bound; see the right panel of Fig. 3c.
DM annihilates into Z and ϕ predominantly in the p-wave, which subsequently decay predominantly into neutrinos. It seems to follow that this DM is safe from stringent constraints from indirect-detection experiments. On the other hand, considering the large Sommerfeld-enhancement factor for the s-wave, we need to be careful about the subdominant-modes of annihilation. A leading s-wave annihilation channel is Considering the current (future) bound of DM annihilation into neutrinos mainly from Super-Kamionkande (Hyper-Kamionkande) [84], (σv) 10 −24 (10 −25 ) cm 3 /s for 20 GeV DM, the Sommerfeld-enhancement factor of S 10 11 needs attention. The electromagnetic decay of Z → e + e − via the kinetic mixing between Z and the SM hypercharge gauge boson needs another attention. Its natural value is = g /70 from the muon-tau lepton loop. The electromagentic branching ratio is Br ∼ ( e/g ) 2 2 × 10 −5 . We may compare (σv)Br with the constraints on DM electromagentic annihilation, e.g., from the CMB measurements [85], (σv) 10 −26 cm 3 /s for 20 GeV DM.
Before concluding this section, we describe the implications and prospects of our results. We again assume the natural value of = g /70 for the kinetic mixing, unless otherwise noted. Since the light L µ − L τ gauge boson predominantly decays into neutrinos, it heats only neutrinos after the neutrino decoupling, increasing the effective number of neutrino degrees of freedom N eff [25,86]. Interestingly, its slight deviation from the SM value, ∆N eff 0.2-0.5 mitigates the tension in the Hubble expansion rate H 0 between the local measurements (i.e., local ladder) [87][88][89] and CMB-based measurements (i.e., inverse ladder) [85,90]. The corresponding Z mass is m Z = 10-20 MeV [91], which coincides with our "prediction" surprisingly. Such a sizable ∆N eff can be examined by CMB-S4 experiments [92,93]. The non-standard interactions of solar neutrinos with electrons and nuclei, in neutrino experiments (e.g., COHERENT) and DM experiments (e.g., LZ and Darwin), examine this light Z region [94,95]. The Z -resonant non-standard interactions of highenergy neutrinos with cosmic neutrino background lead to a "dip" in the high-energy neutrino spectrum [86,96,97], which may explain (though not statistically significant) the null detection of astrophysical neutrinos with 200-400 TeV in IceCube [98][99][100]. The light Z can be probed as missing-energy events in colliders [101,102].

Composite dark matter: asymmetric dark matter
In this section, we consider composite ADM based on the dark QCD×QED dynamics with light dark quarks, up U (+2/3) and down D (−1/3) (dark QED charge) [36]. DM consists of dark nucleons, whose asymmetry has the same origin as the SM baryon asymmetry (i.e., co-genesis). As a consequence, the mass of the dark nucleons is similar to but slightly heavier than that of SM nucleons as indicated from the coincidence of the mass densities Ω dm ∼ 5Ω b . In other words, the dynamical scales of dark QCD and SM QCD are similar to each other: Λ QCD ∼ (10/N g )Λ QCD , where N g is the number of the generations of U and D .
The simple model of Ref. [36] is featured by the (intermediate-scale) neutrino portal and (low-scale) dark photon portal. Non-equilibrium decay of the right-handed neutrinos with the masses of M R > 10 9 GeV generates the whole baryon asymmetry (i.e., thermal leptogenesis [79][80][81][82]). The asymmetry is shared among the dark and SM sectors via the higher-dimensional portal operator. The portal operator originates from the righthanded neutrinos and dark colored Higgs, H C , with the mass of M C . The dark photon is assumed to be the lightest particle in the dark sector. It carries all the dark sector entropy in the late Universe, and transfers it through decay into SM particles via a kinetic mixing with the SM photon. To this end, the dark photon mass is to be m A 10 MeV and the kinetic mixing parameter to be 10 −10 [36,103]. This phenomenological model enjoys compelling ultraviolet physics such as SU(5) SM ×SU(4) dark [104] and mirror SU(5) SM ×SU(5) dark [105]. They involve an intermediate-scale supersymmetry and are discussed in Appendix D.
To discuss the phenomenology, we follow the simplified version of the model given by Ref. [106]. The intermediate-scale portal operators are 11 Here Y C is the Yukawa coupling of the dark colored Higgs to dark U D andŪ D . y N and Y N are the Yukawa couplings of the right-handed neutrinos to SM LH and dark H CD , respectively. The former is related with the observed tiny neutrino masses m ν through the see-saw mechanism [73][74][75][76]: These portal operators are to be relevant at least around the time of the thermal leptogenesis: The first inequality is to prohibit direct decay of the right-handed neutrinos into the dark colored Higgs, which may affect the thermal leptogenesis [107].
Since we consider the QCD-like dark sector, we may employ the QCD-calculation results in the DM phenomenology. The 2-body dark nucleon scattering is described by the effective range theory as σ = 1 2 for the spin singlet and triplet (deuteron) channels, respectively (in SM QCD). 12 We only consider the dark neutron-dark proton (n -p ) scattering even though we also have the other scattering processes: n -n and p -p scattering. 13 We assume that the half of the total DM consists of dark protons. The DM can interact with both of dark proton and dark neutron, and thus we divide each scattering cross sections by two. Here we follow Ref. [32] and fit the QCD-calculation results given in Ref. [108] with Λ = 250 MeV.This is valid for m π /Λ 0-1.4. We approximate the nucleon mass as m N 1.25N c Λ, where the color number is N c = 3 in the present model [114].
As shown in Ref. [36,115,116], if the asymmetry is fully shared among the dark and SM sectors, i.e., if the portal interactions are relevant in the course of the thermal leptogenesis, the DM mass is to be m N 8.5/N g GeV. This is because the B − L asymmetries of the SM and dark sectors follow where N g is the generation of SM fermions and m is the number of light Higgs doublets. The relation between the present B and B − L asymmetries within the SM sector depends 13 The quantum resonance of these scattering processes may require a different dark pion mass from that for n -p scattering since the dark isospin symmetry, among n and p , is broken by U(1) dark and the current dark quark masses. Therefore, we expect that we can ignore the other dark nucleon scattering on the quantum resonance of n -p scattering. When the quantum resonance of p -p and n -n scattering coincides with that of n -p , we have to take into account each scattering cross section, Here the additional factor of 2 originates from the difference in the decomposition of the scattering amplitude into the isospin irreducible representations and from the symmetric factor of the final states. It is beyond the scope of this paper to make a further study of the dark nucleon scattering in detail.
on the nature of the electroweak symmetry breaking (EWSB) [117,118]: A SM , sphaleron decoupling after the EWSB , A SM , rapid sphaleron after the EWSB , A SM , rapid sphaleron + top decoupling after the EWSB .
In the first relation, we assume that the B − L number and hypercharge are conserved, while in the second and third, we assume that the B − L number and electromagnetic charge are conserved. Hereafter, we use the third relation with N g = 3, m = 1, and As seen in Fig. 4, m N 8.5 GeV has a constant cross section below v rel 100km/s, barely reproducing the velocity dependence of σ/m inferred by the MW's dwarf spheroidal galaxies. The binding energies are E bs 0 , E bt 60 MeV for resonant singlet scattering , E bt 0 for resonant triplet scattering , where E b 0 means a shallow bound/virtual state. Meanwhile, m N 16 GeV has a σ/m ∝ 1/v 2 rel in the low-velocity region, providing a better fit. The binding energies are E bs 0 , E bt 110 MeV for resonant singlet scattering , E bt 0 for resonant triplet scattering .
These 2-body binding energies are important inputs for the bound-state formation in the Universe (i.e., the dark nucleosynthesis) [119][120][121][122][123][124][125][126]. In the present model, the dark photon is typically heavier than the binding energy, and thus the 2-body bound-state formation proceeds only through the electron/positron emission and is suppressed by the kinetic mixing parameter. Though it is beyond the scope of this paper, we may also be able to predict the strength of the dark QCD phase transition. If it is the strong first order, we have a chance to examine this model in near-future gravitational-wave detectors [127].
To make m N 16 GeV consistent with Ω dm ∼ 5Ω b , we need to make portal operators partially irrelevant. We leave more general analysis, including the direct decay of righthanded neutrinos into the dark colored Higgs, for a future work (see Ref. [107] for a similar analysis). Hereafter, instead, we assume that the both cases (m N 8.5 GeV, 16 GeV) are realized with M C 10/x In the present model, the dark photon plays a key role in releasing the dark sector entropy to the visible sector via the kinetic mixing with the SM photon. Fig. 5 summarizes the current constraints on the dark-photon parameters (m A , ) and the future prospects in beam-dump and collider experiments. In collider experiments, dark photons can be produced through the kinetic mixing and subsequently decay into SM particles. For m A < 5 GeV, BaBar [132,133], LHCb [134], and KLOE [135][136][137][138] Figure 5: Constraints on the dark-photon parameters: dark photon mass m A and kinetic mixing . The gray meshed region is the current bounds from supernova 1987A constraints (Refs. [128,129]), beam-dump experiments, and collider experiments (taken from Ref. [130]). The pink dashed line shows the future prospects for beam-dump and collider experiments. The orange (red) lines show the upper limits on from the PandaX-II experiment [131] for m N = 16 GeV (m N = 8.5 GeV) and fixed U(1) dark coupling, α D = g 2 D /(4π).
limit on . Belle-II [139] will place the upper limit on mainly for m A = 10-5000 MeV. The displaced vertex searches at the LHCb experiment will also explore the region between the collider and beam-dump experiments [140,141]. The beam-dump experiments exclude a broad parameter space for m A < 1 GeV: current constraints are mainly from CHARM [142], LSND [143], and U70/νCal [144,145]. The dark-photon parameters in this region will be further explored by projected facilities like SHiP [146,147], FASER [148], and SeaQuest [149,150]. The dark photon decay rate is given by where R( √ s) = σ(e + e − → hadrons) σ(e + e − → µ + µ − ) takes into account the hadronic resonances (i.e., vector meson dominance). In the following, we consider the latter for m A > 100 MeV taken from Ref. [151]. Since the late-time decay of the dark photon heats only the electronphoton plasma after the neutrino decoupling, N eff is lowered and the upper bound is put on the dark photon lifetime. The dark photon lifetime τ A = 1/Γ A exceeds 1 s in the blue hatched region of Fig. 5. In the right shaded region of Fig. 5, the dark photon is no longer the lightest particle in the dark sector (the orange region for m N = 16 GeV, while the red region for m N = 8.5 GeV). In particular, the dark photon is heavier than the dark pion, whose mass is indicated by the strong velocity dependence of the 2-body dark nucleon scattering. The constraints on the dark photon decay would change in this region. In the left shaded region, the dark photon mass is less than the binding energy E bt (the color similarly indicates the DM mass as that in the right shaded region). In this region, the dark nucleosynthesis would proceed by emitting the dark photon, and then the DM would consists of the dark nuclei.
When the DM consists of the dark proton charged under U(1) dark , DM interacts with SM nuclei through the kinetic mixing between the SM photon and dark photon. The differential cross section between the dark proton and target nucleus is Here the momentum transfer is q 2 = 2µ 2 T v 2 rel (1 − cos θ cm ), where µ T is the reduced mass of DM and the target nucleus and θ cm is the center-of-mass scattering angle. Following Ref. [152], we set q 2 = 2µ 2 T v 2 , where v = 232 km/s is the typical DM speed at the rest frame of the Earth [153]. F T is the nuclear form factor. We discuss the dark neutrontarget nucleus scattering through the magnetic moment in Appendix E, which is negligible. We place the upper limit on the kinetic mixing parameter from the direct-detection constraints by the 54 ton-day exposure of the PandaX-II experiment [131] (orange lines and red lines in Fig. 5). We assume that the half of the total DM consists of dark protons. γ -r a y c o n s t r a i n t s S u p e r -K a m i o k a n d e The DM mass is set to be 16 GeV (orange lines) and 8.5 GeV (red lines), and we take U(1) dark coupling to be α D = α EM = 1/137 (dot-dashed lines) and α D = α EM /100 (solid lines).
The intermediate-scale portal operators Eq. (20) lead to the DM decay into the dark meson and anti-neutrino [116]. No excess of a neutrino signal has been found at the Super-Kamiokande experiment [154], and this sets the lower limit on the DM lifetime [155]. The decay rate is given by Here, v H is a VEV of the SM Higgs and W is a matrix element for DM to a dark meson. We assume the matrix element to be W 0.1 GeV 2 (m N /m n ) 2 , to which a lattice calculation result [156] is naïvely scaled. For the decaying DM with a mass of 16 GeV, the lower limit is roughly τ DM = 1/Γ DM 10 23 s. We naïvely extrapolate the results of Ref. [155] below the DM mass of 10 GeV. The black-shaded region in Fig. 6 is excluded by no excess of a neutrino signal in the Super-Kamiokande experiment. The same ultraviolet origin of the operator in Eq. (20) induces the Majorana mass term for the dark neutrons [106]: Some portion of the dark neutron DM is converted into its anti-particle at the late time t through the dark neutron-anti neutron oscillation: Pn ∼ (m M t) 2 . Then dark nucleons n and p annihilate withn into dark pions with the effective cross section, Here S eff is the Sommerfeld-enhancement factor S s-wave times a fudge factor. When the dark photon is the lightest particle in the dark sector, dark pions decay into dark photons, and then dark photons eventually decay into e + e − pair. The final state radiation of this process can emit visible photons. The late-time annihilation is constrained by the observations of the γ-ray from the dwarf spheroidal galaxies by the Fermi-LAT and from the interstellar e + e − flux by the Voyager-1. We use the lower limits on the oscillation time scale given by Ref. [106] with S eff = 1. The green-hatched region in Fig. 6 shows the constraints from the observations of the γ-ray flux. The constraints from the γ-ray observations depend on the annihilation cross section of dark nucleons. The strong velocity dependence of the self-scattering cross section implies a sizable Sommerfeld-enhancement factor as seen in Fig. 2a, and thus we also show the lower limit on M R when we take the factor to be S s-wave = 10 3 (long-dashed line) and S s-wave = 10 5 (short-dashed line).

Conclusion
We have studied the possibility that self-interacting dark matter has a maximal selfscattering cross section from the Unitarity point of view. The maximally self-interacting dark matter can explain the observed structure of the Universe ranging from σ/m ∼ 0.1 cm 2 /g at galaxy clusters to ∼ 100 cm 2 /g at MW's dwarf spheroidal galaxies. We have demonstrated that a general requirement of the quantum (zero-energy) resonant scattering, by employing the effective-range theory and analytic results with the Hulthén potential. It requires a light force mediator and model parameters to be the special (fine-tuned) values. We have also demonstrated that DM annihilation tends to be largely boosted by the Sommerfeld enhancement with the same parameters.
It requires DM annihilation to end up with neutrinos or to be the p-wave. As such a model, we consider the gauged L µ − L τ model extended with Dirac DM. The Dirac DM is slightly split into two Majorana DM by the VEV of the gauged L µ −L τ Higgs. We take the Higgs mode to be light so that it mediates the (relatively) long-range force between DM particles. We have pinned down the parameter points that lead to the DM self-scattering saturating the Unitarity bound. With them, we can explain the discrepancy in the muon g − 2, predict the sizable ∆N eff mitigating the tension in H 0 , and have rich implications for non-standard neutrino interactions and collider searches.
ADM is another option to evade the indirect-detection bounds from the boosted DM annihilation (at the first sight). Since ADM requires a lighter state to deplete its symmetric components through efficient annihilation, such a light state is a natural candidate for a force mediator. We have considered dark nucleon ADM with dark pions being the force mediator, based on the dark QCD×QED dynamics. The dark pion mass is very predictive to explain the quantum-resonant self-scattering between dark nucleon DM. For viable cosmology, we assume that dark photon is the lightest state in the dark sector, which provide various ways to probe the dark sector: collider and beam-dump experiments and DM direct-detection experiments. Furthermore, the binding energies of two nucloeon states are also predicted, which are important for the dark nucleosynthesis that proceeds with a light dark photon. The predicted dark vector resonances can be measured by lepton colliders.
We have shown that the recent astronomical data for structure formation of the Universe are already interesting and have rich implications for DM physics. Though it is to be clarified which data are robust or suffer from astronomical uncertainties such as supernova feedback, we are optimistic that the situation will get quickly improved with the fast development of hydrodynamic simulations and more precise observations. Gravitational probes of DM and related beyond-WIMP DM scenarios enjoy an exciting era.

A MW's dwarf spheroidal galaxies vs field dwarf galaxies
As seen in Fig. 1a, the inferred values of σ/m, in the MW's dwarf spheroidal galaxies (cyan) [44] and field dwarf galaxies (red) [3], seem not compatible with each other, though they have similar collision velocities, v rel ∼ 20-60 km/s. The former prefer a larger value, σ/m ∼ 30-100 cm 2 /g, while the latter prefer a smaller value, σ/m ∼ 0.1-20 cm 2 /g. It might be because Ref. [3] considers the long mean free path regime. If we repeated their analyses in the short mean free path regime, we might find a larger value of σ/m. In the main text, we discuss implications for DM phenomenology, taking the former seriously. This would be reasonable, if we took the error bars at face values. It is beyond the scope of this paper to explore a possible origin of this tension between the MW's dwarf spheroidal galaxies and field dwarf galaxies. In this section, instead, we discuss how our conclusions changed, if σ/m in the MW's dwarf spheroidal galaxies was lowered.
In the left panel of Fig. 7a, we take the effective-range theory parameters so that σ/m crosses the inferred values from the field dwarf galaxies (red), rather than the ones from MW's dwarf spheroidal galaxies (cyan) as in Fig. 1a. In this case, we do not need a strong velocity dependence of σ/m as in Fig. 1a. Notice that the |a/r e |'s took are relatively smaller than the ones in Fig. 1a; this indicates that the parameters are relatively away from the quantum resonance. This is explicitly shown in Fig. 7b, where the points depicted by stars represent the Hulthén-potential parameters that corresponds to the effective-range theory parameters in left panel of Fig. 7a. Notice that the φ 's do not need to be so close to the first quantum resonance, i.e., φ 1, to cross the red data points in Fig. 7a. This is in contrast to the parameters took in Fig. 1a, where the effective-range theory parameters need to be fine-tuned to be near the quantum resonance to cross the cyan data points, as shown in Fig. 1b. At the same time, the parameters of Fig. 7a are also relatively away from the quantum resonance for annihilation; compare the points depicted by stars in Fig. 8 with the ones in Fig. 2. Consequently, they exhibit relatively smaller Sommerfeld-enhancement factors than the ones of Fig. 1a, making them less constrained from the indirect-detection experiments.
The lack of need for the strong velocity dependence amounts to relaxation of the "prediction" of narrow model parameter regions of the gauged L µ − L τ model in the main text. In the left panels of Fig. 9, we depict the example parameters of the model that explains the inferred values of σ/m from the field dwarf galaxies, i.e., σ/m ∼ 1 cm 2 /g at v rel ∼ 20 km/s. We show the corresponding velocity dependence of σ/m in the right panels. In Fig. 9a, we present the case of m Z = 10 MeV, which lies in the Z mass range preferred to explain the inferred σ/m from MW's dwarf spheroidal galaxies, m Z 20 MeV, as discussed in the main text. Contrary to the parameter points depicted in Fig. 3a, we see that the parameters here do not need to be so fine-tuned to be close to the quantum resonance (yellow) to explain the red data points. Furthermore, the viable Z mass range extends towards larger values, compared to the case discussed in the main text. For heavier Z , e.g., m Z = 50 MeV, it is possible to find parameters to explain the inferred σ/m from field dwarf galaxies, while they need to lie near the quantum resonance.
For the ADM based on the dark QCD×QED dynamics, the notable change would be that the dark nucleon mass of m N 8.5 GeV also becomes preferable. In Fig. 4, we showed that m N 8.5 GeV poorly fits the inferred velocity dependence of σ/m inferred by the MW's dwarf spheroidal galaxies. But if we focus on σ/m inferred by the field dwarf galaxies, m N 8.5 GeV also provides a good fit to the inferred velocity dependence (see Fig. 10), where we do not need to consider additional entropy production to dilute the DM abundance or consider asymmetric generation of B − L number. Meanwhile, the preferred dark pion masses in both Fig. 4 and Fig. 10 do not differ much, i.e., the pion mass (or m π /Λ) still has to be reasonably tuned near the quantum resonance for elastic scattering. Similarly, the binding energy of nucleons, E bt , decrease by a factor of ∼ 0.9 as we take pion mass slightly away from the quantum resonance to explain the velocity dependence of σ/m inferred by the field dwarf galaxies. This would slightly retreat the left shaded region in Fig. 5 to the further left, in which the dark nucleosynthesis may occur by emitting dark photons.

B Higher resonances
As seen in Fig. 1b, the quantum resonances for elastic scattering appear at φ = n 2 (n = 1, 2, . . . ). In the main text, we focused on the parameters near the first one, n = 1. While it is certainly encouraging to discuss the higher resonances as well, there are a couple of reasons we try not to. In this section, we discuss these reasons and remark on the possible difficulties and issues in investigating higher resonances.
The first aspect of the difficulty to explore higher resonances is related to the method of comparing astronomical data with velocity-dependent σ/m. For an elaborate discussion on this aspect, we comment on how astronomical data are interpreted by SIDM to infer the data points presented in Fig. 1a. The horizontal axis of Fig. 1a, v rel , is the average relative velocity between two DM particles. Consider the collision between particle 1 and 2. Assuming the Maxwellian velocity distribution for both particles, but with different 1d velocity dispersion: σ 2 1 and σ 2 2 , respectively. After integrating out the center-of-mass velocity, whose 1d velocity dispersion is σ 2 cm = , we obtain the distribution function for the relative velocity: where the 1d velocity dispersion is related with the average relative velocity as In the gravitationally interacting system, we expect σ 1 = σ 2 = σ and thus v rel = 4σ/ √ π as found in the literature [3]. Note that this differs from the kinetically equilibrium system, where m 1 σ 2 1 = m 2 σ 2 2 . Meanwhile, the inferred values of σ/m at given v rel are achieved by assuming a constant self-scattering cross section inside a galaxy of interest. In the case of velocity-dependent σ/m, a fairer comparison may be done by taking the distribution average for the cross section.  Figure 7: (a) Same as Fig. 1a, but for the effective-range theory parameters, (a, r e ), that explains the inferred values of σ/m from the observations on field dwarf galaxies (red). The Hulthén-potential parameters took in the right panel exhibits the same effectiverange theory parameters as the left panel; (α, m φ ) = (7.0 × 10 −3 , 0.1 GeV) (solid) and (α, m φ ) = (3.5 × 10 −4 , 0.01 GeV) (dashed). (b) Same as Fig. 1b, but the stars indicate the benchmark parameters took in Fig. 7a. Compared to the parameters took in Fig. 1a, the parameters here are relatively away from the quantum resonance at φ = 1.  Figure 8: Same as Fig. 2, but with the Hulthén-potential parameters took in Fig. 7a. The s-wave annihilation is relatively away from the quantum resonance compared the ones in 1a; the enhancement factors are smaller, and the s-wave enhancement factors grow towards low velocity slower than the resonant behavior, i.e., ∝ 1/v 2 rel (brown).  Figure 9: Same as the left panel of Fig. 3, but the with different benchmark parameter points; here, the depicted points explain the red data points, rather than the cyan ones as in Fig. 3.  Figure 10: Velocity dependence of the 2-body dark nucleon scattering for given dark nucleon mass m N and dark pion mass m π . We take the dark pion mass to make the velocity dependence of scattering consistent with the field dwarf galaxies. We use the same coloring of lines as in Fig. 4, and the data points are the same as in Fig. 1a.
In the main text, we do not take the distribution average for the cross section. This is partially just for the simplicity. Another reason is that when astronomical data are interpreted by SIDM, a self-scattering cross section is assumed to be constant. Thus, to be consistent, we need to reanalyze the data by taking into account the velocity (dispersion) dependence of the (distribution averaged) cross section. This is beyond the scope of this paper, and we do not expect a significant difference, as long as we consider the effectiverange theory parameters near the first resonance.
At higher resonances for elastic scattering, σ/m could exhibit diminishing values at specific momentum. This is demonstrated in the right panel of Fig. 11a, where we show the analytic result of σ/m for the Hulthén potential. There, we see that for the effectiverange parameter set near the second resonance, i.e., the blue star of Fig. 11b, σ/m is diminishing around v rel ∼ 700 km/s. This is related to the behavior of the phase shift δ 0 at the second resonance. At the vicinity of the second resonance, δ 0 starts out from 0 at high-k limit and approaches 3π/2 in the k → 0 limit; between the two limits, δ 0 = π for some specific value of k and σ/m diminishes at such k. One may expect that the existence of this "dip" generally makes the parameters hard to be compatible with the data points in Fig. 11a. However, this incompatibility may be mitigated once we take the distribution average for the cross section, i.e., the "dip" becomes relaxed. Again, we do not attempt to explore this issue for the simplicity of the discussion. Furthermore, in order to consistently study the compatibility of the parameters with observations, the existence of the "dip" in the velocity dependence may require a more delicate reanalysis on astronomical data as discussed above.
The second aspect is related to the Sommerfeld enhancement of DM annihilation. As we consider the resonances beyond the first one for elastic scattering, higher partial wave resonances for annihilation may become important. As an example, in Fig. 12, we show the velocity dependence of the Sommerfeld-enhancement factors for the parameter point near the second resonance (blue star), which corresponds to the parameter took in the   Figure 11: Same as Fig. 1a, but comparing the parameter points close to the first (cyan star) and the second (blue star) quantum resonance for elastic scattering, i.e., n = 1 and n = 2 for φ = n 2 , respectively. In the right panel of Fig. 11a, for the parameter point near the second resonance (dashed), the analytic results for the for the Hulthén potential exhibits diminishing σ/m at specific v rel ∼ 700 km/s.  Figure 12: Same as Fig. 2, but for the parameter points depicted in Fig. 11b, which are close to the first (cyan star) and the second quantum resonance (blue star) for elastic scattering. In the left panel, the brown curve corresponds to ∝ 1/v 2 rel . Contrary to the cyan data point, which was discussed in Fig. 2, the blue data point is close to both s-wave and p-wave resonances for annihilation. right panel of Fig. 11a. Coincidentally, the parameter lies near the s-wave and p-wave resonance for annihilation. For this parameter point, taking p-wave annihilating DM to evade the constraints from indirect-detection experiments (as for the cyan parameter point) may not help much. Likewise, the implications on the constraints may highly depend on the value of the took φ for higher resonances of elastic scattering. A dedicated investigation on this aspect for higher resonances may be done elsewhere, while we try to focus on the simplest case in this work.
C Neutrino masses in the gauged L µ − L τ model The renormalizable Lagrangian density can be written as Here, L i (i = 2, 3) denotes the left-handed leptons in the flavor basis, whileμ andτ denote the right-handed charged leptons. After the SM Higgs H and L µ − L τ breaking Higgs Φ develop the VEVs, v H / √ 2 and v Φ / √ 2, the resultant neutrino mass sector takes a form of Here (y eµ , y eτ ). We fix the phases of N l so that M eµ , M eτ , M µτ > 0, while leaving the phases of ν l for the PMNS parametrization [159,160]. The see-saw mechanism [73][74][75][76] provides the neutrino mass term at low energy as where the PMNS matrix is parametrized by with θ ij ∈ [0, π/2] and δ, α i ∈ [0, 2π). and the mass eigenvalues are written as for the normal ordering, while m 3 < m 1 < m 2 for the inverted ordering. The neutrino oscillation parameters are summarized in Table 1 Table 1: Neutrino oscillation parameters from Ref. [77,161].
It is remarkable that the U(1) Lµ−Lτ symmetry restricts M D to be a diagonal matrix. Thus if M D has a inverse matrix, the following relation holds: Since M −1 D is also diagonal, the flavor structure of M −1 ν should follow the structure of M N . Two zero entries of M N (see Eq. (36)) provides four constraints for the PMNS parameters [77]: δ π/2 or 3π/2, m 0.05-0.1 eV, α 2 /π 0.6, and α 3 /π 1.4 from δm 2 , ∆m 2 , and θ ij given in Table 1. The other entries provide 8 constraints, while there are 11 model parameters (6 from Y l and 5 from M ll ). This is because Eq. (40) is invariant under the following transformation: , c being real. We obtain Table 2 by fitting the other eight model parameters in Eq. (40). The mass matrices are rewritten as Parameter best fit in Table 1 Y

D Supersymmetric realization of composite ADM
A supersymmetric realization of composite ADM scenario [104,105] may cause the fast DM decay as in nucleon decay in supersymmetric grand unified theories [162,163]. Heavy particles in ultraviolet physics induce the intermediate-scale portal operators with the mass dimension-six rather than the mass dimension-seven, The equilibrium condition now reads The decay rate is given by M S and M 1/2 denote typical mass scales of supersymmetric scalar particles and supersymmetric fermionic particles, respectively. The intermediate-scale portal operator Eq. (42) vanishes when the generation N g = 1 due to anti-symmetrization over color indices. Hereafter, we assume N g = 2.
The supersymmetric particles change the relation between the B − L asymmetries in the SM and dark sectors: γr a y c o n s tr a in ts Figure 13: Constraints on M R -M S plane in a supersymmetric realization of composite ADM. The green-hatched region and the green-dashed lines are the same as in Fig. 6. The LSP abundance from the gravitino decay exceeds the observed DM abundance in the meshed-purple region in thermal leptogenesis scenario.
As with non-supersymmetric realization, we also have the Majorana mass term for the dark neutrons in the supersymmetric realization. The corresponding superpotential has the mass dimension seven, and is given by This operator induces the dimension nine operator at the mass scale of supersymmetric particles with two-loop diagrams [164], and gives the Majorana mass term for the dark neutrons. A typical size of the Majorana mass is Fig. 13 shows the indirect-detection constraints on the M R -M S plane from a neutrino signal and γ-ray observations. In this figure, we assume that the loop factor is L = 0.01, and M 1/2 = 10 −4 M S that is a spectrum motivated by split supersymmetry scenario [165][166][167]. The black shaded region in Fig. 13 is excluded by no excess of a neutrino signal at the Super-Kamiokande experiment. We also show the future prospect of the neutrino detector experiments, Hyper-Kamiokande, with scaling by a factor of ten of the current bound, which corresponds to the orange-dotted line. The green-hatched region is excluded by the γ-ray observations. In supersymmetric realization, the lightest supersymmetric particle (LSP) can be also stable due to the R-parity. We have the LSPs in each of visible and dark sectors, and then the heavier state can decay into the lighter one through portal interactions at the late time. The supersymmetric counterpart of the kinetic mixing between the dark photon and SM hypercharge gauge boson leads the prompt decay of the heavier state to the lighter one when the higgsinos are the lightest particles in each of sectors [104]. Here, the higgsino in the visible sector, called visible higgsino, is the supersymmetric partner of the SM Higgs, while the one in the dark sector, called dark higgsino, is that of the U(1) dark breaking Higgs. Even though the dark higgsino is less constrained than the visible higgsino, the LSP abundance from the gravitino decay would cause the overclosure of the Universe [168,169]. The purple-meshed region in Fig. 13 shows the LSP abundance exceeds the observed DM abundance when T R M R and the LSP with a mass of M 1/2 .
The constraints from the LSP abundance may be naturally relaxed in the models with an intermediate-scale dark grand unification: SU(5) SM ×SU(4) dark [104] and mirror SU(5) SM ×SU(5) dark [105]. They associate an intermediate-scale dark monopole in the SU(4) dark →SU(3) dark × U(1) dark phase transition. Their abundance is determined by the pair-annihilation as [170,171] Here we use g * = g * MSSM + g * dark + g * N and g * Note that this expression is valid only for S after /S before > 1; otherwise, S after /S before = 1 (i.e., no entropy production). Here T c ∼ m A /g D is the critical temperature of the U(1) dark phase transition. T ann is the temperature at which monopole annihilates. In the second equality, we take g * (T c ) = g * (T ann ). T eq is the temperature at which monopole domination begins: In the above discussion, we assume that the SM and dark sectors are in equilibrium in the course of the monopole domination, while it may not be valid for a tiny kinetic mixing. Monopole annihilation may also produce the LSP [172], while we leave further analysis for a future work.

E Direct detection through the magnetic moment of the dark neutron
Similarly to the dark proton scattering with SM nuclei, the dark neutron carries the magnetic moment under U(1) dark : Here we take the magnetic moment to be µ n = (g n /2)(e D /2m N ) with the g-factors being g n = −3.83 (g p = 5.59 for the dark proton) in analogy to the SM nucleons. The matrix element for the dark neutron-SM proton scattering is 15 where M p is the matrix element for the dark proton-SM proton scattering, and Here, s p(n ) is the spin vector of the SM proton (dark neutron, i.e., DM) and v ⊥ = v − q/(2µ N ). We can further evaluate the matrix element for the dark neutron-target nucleus scattering with the form factors of F N 1 N 2 i,j [173,174], while we leave it for a future work. In orders of magnitude, the direct-detection bounds on from dark-neutron scattering is weakened by a factor of 1/v 2 ∼ 10 6 compared to the dark proton-target nuclei scattering.