Dark matter from strong dynamics: the minimal theory of dark baryons

As a simple model for dark matter, we propose a QCD-like theory based on SU(2) gauge theory with one flavor of dark quark. The model is confining at low energy and we use lattice simulations to investigate the properties of the lowest-lying hadrons. Compared to QCD, the theory has several peculiar differences: there are no Goldstone bosons or chiral symmetry restoration when the dark quark becomes massless; the usual global baryon number symmetry is enlarged to SU(2)B, resembling isospin; and baryons and mesons are unified together in SU(2)B iso-multiplets. We argue that the lightest baryon, a vector boson, is a stable dark matter candidate and is a composite realization of the hidden vector dark matter scenario. The model naturally includes a lighter state, the analog of the η′ in QCD, for dark matter to annihilate into to set the relic density via thermal freeze-out. Dark matter baryons may also be asymmetric, strongly self-interacting, or have their relic density set via 3 → 2 cannibalizing transitions. We discuss some experimental implications of coupling dark baryons to the Higgs portal.


Introduction
The mass and stability of luminous matter in the Universe are largely a byproduct of QCD. Around 99% of the mass of baryonic matter arises from the strong interaction and its stability is a consequence of an accidental U(1) B baryon number symmetry in the Standard Model (SM). Yet baryons represent only about one sixth of the total mass budget of matter in the Universe. The remainder is dark matter (DM). Although unknown in its properties, many different models and detection strategies have been proposed for DM, often motivated by other issues associated with the SM (e.g., the stability of the weak scale or the strong CP problem) [1,2].

JHEP12(2018)118
On the other hand, dark and luminous matter may come from two separate particle physics sectors, orthogonal to any problems of the SM [3]. DM may be the lightest stable state within a dark sector that has its own gauge group and matter representations. If these fields are singlets under the SM gauge group, the interactions between the two sectors only arise through higher dimensional operators and may be very feeble. However, interactions within the dark sector are generically much larger, especially for a non-abelian gauge theory that is strongly coupled. Moreover, long-standing puzzles on the galactic-scale structure of DM provide an astrophysical motivation for strong self-interactions between DM particles [4] (see ref. [5] for a review). Composite dark sectors are a natural framework for self-interacting DM [6,7].
It is an appealing hypothesis that the mass and stability of DM arise through strong dynamics, similar to luminous baryons. Early realizations along these lines include technicolor baryons [8][9][10] and mirror baryons [11,12]. For strongly coupled gauge theories, lattice field theory is the main calculational tool in the nonperturbative regime. While mainly used for QCD [13], recent studies of non-abelian dark sectors have turned to the lattice to investigate the basic properties of these theories, such as the spectrum of states and form factors for interactions with the SM [14][15][16][17][18][19][20][21][22]. We refer the reader to ref. [23] for a recent survey of different models along these lines.
In this work, we propose a minimal model realizing these ideas and compute its basic properties on the lattice. The DM candidate in our model will be the lightest baryon in a strongly coupled Yang-Mills theory, which is stable due to an accidental symmetry. By minimal, we mean the fewest number of colors N c and flavors N f , and the smallest nontrivial representation for matter fields. Hence, we consider SU(2) gauge theory with one Dirac fermion q (dark quark) in the fundamental representation. We do not consider the case of a single Weyl fermion due to Witten's anomaly [24]. Hambye and Tytgat proposed a similar DM model based on SU(2) gauge theory with scalar quarks [25].
In the space of gauge theories, N c = 2 theories have long been useful as a simplified version of QCD [26][27][28][29]. However, an important distinction is the fact that the fundamental representation of SU(2) is pseudo-real, unlike in SU(N c ) with N c > 2. As a consequence, two-color theories have an enlarged global symmetry that reflects transformations between quarks and antiquarks. At the hadron level, there is a unification of baryons, antibaryons, and mesons.
In our one-flavor theory, the quark q and antiquarkq fields form a doublet, written schematically as Q ∼ q q . As we show in section 2, the theory has an unbroken global SU(2) B symmetry acting on Q. 1 This symmetry is a non-abelian generalization of a U(1) B baryon number symmetry for q; it is clear that U(1) B is the diagonal subgroup of SU(2) B . Since the whole setup is analogous to isospin, we refer to this symmetry as baryonic isospin.
We argue below that SU(2) B is not violated by chiral symmetry breaking or a finite mass for q. Hence, the hadronic spectrum of the theory will fall nicely into SU(2) B iso-multiplets.

JHEP12(2018)118
We envision that the lightest baryon in our theory will be a suitable DM candidate. The lightest qq state is part of a spin-1 iso-triplet Borrowing an analogy from QCD, this state is akin to the ρ meson. However, the superscripts in eq. (1.1) refer not to electric charge, but to U(1) B charges: baryon (+), antibaryon (−), and meson (0). All three components are stable DM candidates provided SU(2) B remains unbroken. Another peculiar feature of our N f = 1 model is the absence of Goldstone bosons. Once the axial U(1) A anomaly is considered, no chiral symmetries are present in the "chiral" limit, where q becomes massless (see discussion in ref. [30]). The would-be Goldstone boson from the U(1) A symmetry, which we denote η (analogous to the η in QCD), acquires a mass through the anomaly. This stands in contrast to SU(2) gauge theory with N f = 2, which has an enlarged pion sector compared to QCD and the lightest baryons are themselves Goldstone bosons [14].
In the early Universe, strong interactions in the dark sector populate a thermal plasma of dark quarks and gluons, which later are confined into hadrons after a cosmological phase transition, similar to QCD. The DM relic density may be frozen-out before or after the transition, depending on the dark quark mass m q and the confinement scale Λ MS . In the latter case (m q Λ MS ), an appealing feature of our model is that there is a builtin annihilation channel ρρ → ηη for setting the relic density provided m ρ > m η (with the η subsequently decaying into SM particles). It is one of our key lattice results that this inequality holds for any value of m q , unlike QCD where m ρ < m η . Annihilation is important for standard freeze-out [31] or asymmetric freeze-out where the dark sector has a dark baryon asymmetry [32,33]. The precise details depend on the relative temperature of the dark sector [34], its coupling with the visible sector (see, e.g., [35]), and the possible role of cannibalizing transitions [36,38], such as ρρρ → ρρ. We defer an analysis of the cosmology of our model to future study.
The remainder of this work is organized as follows. In section 2, we present our dark sector model, including the leading non-renormalizable operators with SM fields. We discuss the SU(2) B symmetry properties of the Lagrangian and other bilinear operators that will be relevant for the lattice computations. We also discuss recent arguments that SU(2) gauge theory does not provide a suitably stable DM candidate [20] and argue that the ρ meson in our model avoids these pitfalls. Section 3 describes the lattice ensembles that we use, how quark propagators are constructed, and provides a first look at the hadron spectrum. We devote particular attention to defining the "chiral" m q = 0 limit in this Goldstone boson-less theory. We also compute Λ MS as a convenient scale to normalize dimensionless lattice quantities into physical units. Section 4 presents our main results: the calculation of the dark hadron mass spectrum and decay constants. In section 5, we discuss couplings between the dark sector and the SM and implications for DM detection. In particular, we use the Feynman-Hellman theorem to provide a determination of the

JHEP12(2018)118
Higgs coupling to our dark matter candidate. Conclusions are provided in section 6. The appendices describe two complementary methods to determine the "chiral" point, provide an alternative and more precise approach to defining physical scales, and summarize our lattice ensembles.

Renormalizable Lagrangian and bilinear operators
The Lagrangian for SU(2) gauge theory with one Dirac fermion q, with mass m, is We assume that q is in the fundamental representation of SU(2) and is a singlet under the SM gauge symmetries. The covariant derivative is D µ = ∂ µ + i 2 gA a µ σ a , where g is the gauge coupling and σ a represents the Pauli matrices acting on SU(2) color indices. Although there are no renormalizable interactions between the dark sector and the SM, the two sectors may couple through higher dimensional operators, which we discuss below.
By analogy with QCD, eq. (2.1) has a U(1) L × U(1) R chiral symmetry for m = 0. However, the two-color theory is different from QCD since the fundamental representation of SU(2) is pseudo-real. Our theory possesses an enlarged U(2) global symmetry. 2 To see this, we can write the fermion part of eq. (2.1) in the following form where C is the charge conjugation matrix acting on Dirac spinors and 3) The kinetic term in eq. (2.2) is manifestly invariant under U(2) transformations acting on Q. For the mass term, let us decompose the global symmetry as U(2) = U(1) A × SU(2) B , since rotating Q by an overall phase is equivalent to an axial U(1) A transformation on q.
As mentioned in the introduction, SU(2) B is a baryonic isospin symmetry, with U(1) B as a subgroup, that plays a similar role as isospin in QCD. While U(1) A is broken for m = 0, SU(2) B remains intact since E is an invariant tensor. In lattice calculations, local operators constructed from q,q create and annihilate states in the hadronic spectrum with the same quantum numbers. In this work, we consider states with J P = 0 ± and 1 ± . The relevant mesonic operators are

JHEP12(2018)118
On the right-hand side, we have expressed these operators in terms of Q to make clear the SU(2) B isospin properties of these states. The scalar, pseudoscalar, and axial vector operators are iso-singlets. To write the vector operator, we introduce Pauli matrices τ a acting on isospin indices. The vector operator is part of an iso-triplet that includes both meson and diquark operators. From eq. (2.5), we see that the lightest baryon in the theory has J P = 1 − and forms a triplet under SU(2) B , described in eq. (1.1). We also write the tensor bilinear as By analogy with QCD, we expect the chiral condensate qq to receive a nonzero value through spontaneous symmetry breaking. However, sinceqq is an iso-singlet operator, its vacuum expectation value does not violate baryonic isospin. On the other hand, the chiral condensate breaks the global U(1) A symmetry, potentially leading to a pseudo-Goldstone boson η that becomes massless for m = 0. However, just as in QCD, U(1) A is anomalous, which gives an additional contribution to the η mass.

Nonrenormalizable interactions and CP violation
The dark sector and SM may be coupled through higher-dimensional operators. The leading operators, arising at dimension five, are O S,P |H| 2 , where H is the SM Higgs field. When the Higgs field gets its vev H = v/ √ 2, there is an additional contribution to the dark quark mass. In general, this term need not be aligned with the Dirac mass m and there may be a relative CP-violating phase between them. 3 If we start in a basis where only O S |H| 2 appears, we must allow m to be complex: Here, M is the mass scale parametrizing the coupling between the two sectors. Performing a chiral rotation to make the total quark mass real and positive, we have where m q = |m + 1 2M v 2 | and φ = arg(m + 1 2M v 2 ). CP violation manifests as a coupling between both operators O S,P and the Higgs boson h. The pseudoscalar coupling is particularly important since it causes the η meson to be cosmologically unstable. Phenomenological consequences of eq. (2.8) are explored in section 5.

Dark matter stability
Let us now discuss the question of whether the lightest baryon in our theory, the ρ meson, provides a suitable DM candidate. Ref. [20] argued that if DM is stabilized by an accidental symmetry, the symmetry must be preserved including operators of dimension five, not just at the renormalizable level. Dimension-five operators, even if suppressed by the Planck scale, may induce DM to decay much more rapidly than the age of the Universe. On the other hand, dimension-six operators lead to a cosmologically acceptable DM lifetime if the suppression scale M is large enough (but below the Planck scale). This is the same situation as the proton in the SM: since the leading operators contributing to proton decay arise at dimension six, protons are cosmologically stable for M 10 13 GeV. According to ref. [20], this argument disfavors SU(2) dark sectors since the global U(1) B symmetry may be violated by dimension-five operators of the form ∼ qq|H| 2 .
However, these arguments do not apply to our N f = 1 model. The only dimension-five operators are O S,P |H| 2 and neither allow for ρ decay since they do not violate SU(2) B . The leading operators that violate SU(2) B must involve O aµ V and arise at dimension six or higher by Lorentz symmetry. Therefore, we conclude that the ρ meson is a viable DM candidate in terms of its stability, while iso-singlet states, such as the η meson, are not.
In fact, the scale of physics connecting the dark and visible sectors need not be extremely high (M TeV) to preserve ρ stability. For example, if the two sectors are coupled through a singlet scalar field, SU(2) B is still preserved since the scalar may only couple to O S,P . Alternatively, if a Z gauge boson mediates the coupling, it may couple to O aµ V . This will break the SU(2) B down to its U(1) B subgroup and, while the ρ 0 will be destablized, the ρ ± remains a stable DM candidate. Hence, DM stability is robust in the face of these simplest mediators between sectors.

Lattice ensembles and propagators
For our lattice study, we discretize the one-flavour SU(2) theory of section 2 to arrive at the familiar Wilson action, where U µ (x) is the SU(2) gauge field and ψ(x) is the 4-component Dirac spinor for the dark quark. The sum over x covers the entire lattice and in this work we primarily use V = L 3 × T = 12 3 × 32, but for certain topics we will use 12 3 × 48 lattices also. We choose to use this arguably small volume for our exploratory study compared to those of typical lattice QCD simulations as we want to cover a large range of bare input masses m 0 on reasonable resources. For the light quark spectrum, we have to compute costly disconnected contributions requiring large numbers of propagator inversions to extract a JHEP12(2018)118 signal. As the computational cost of these inversions scales like some power (V n , n > 1) of the volume and require more iterations as the quark mass decreases, a small volume was deemed a necessity to broadly and accurately map the spectrum for a large range of quark masses. The investigation of the finite volume effects from our volume and the finite lattice spacing effects will be left for a future study. The bare gauge coupling β = 4/g 2 is a function of the lattice spacing, which serves as the ultraviolet cutoff. For this study we choose to work at fixed bare gauge coupling of β = 2.2. The physical scale of our theory can be defined by matching to a known phenomenological scale. The bare quark mass m 0 (which is typically a negative number) gets shifted by additive renormalization, so the massless limit can only be found from the results of numerical simulations. We calculate with several different values of m 0 as listed in table 3. Also shown in table 3 are the number of configurations in each ensemble. Ensembles were generated using the RHMC algorithm [39].
After these ensembles have been generated, the largest remaining expense is the calculation of quark propagators, which requires inversion of a large-but-sparse matrix, For many applications only one row of the inverse is required and then it is sufficient to solve the eigenvalue problem and we choose the source to be a time-diluted [40] Z 2 -stochastic wall (Z2SEMWall) [41] source. Unfortunately, one row is not sufficient whenever two quarks within a single operator can annihilate. For these disconnected diagrams, we use an unbiased stochastic estimator, i.e. time and spin dilution [42,43]. We find that 64 stochastic "hits" per configuration was beneficial for reducing noise at reasonable cost, which is a finding similar to [44]. From the configurations listed in table 3 for our lightest quark masses we do approximately O(64, 000 → 330, 000) inversions for each ensemble's spectrum measurement.

The lightest hadrons
To create a hadron on a lattice, we select the appropriate operator from eqs. (2.4). We create the state at some initial Euclidean time and destroy it at some different time, so the resulting correlation function is given by Hadron masses m n can be obtained by fitting the lattice data to this functional form. The calculation for our dark matter candidate, the ρ ± of eq. (1.1), is straightforward as we can choose an operator where each quark propagator runs from source to sink, but other hadrons are much more costly due to the ability ofqq to annihilate within a single operator.

JHEP12(2018)118
We also compute the decay rates of dark sector states. Notice that eq. (3.4) contains 0|O 1 |n , which is proportional to the hadron's decay constant and can be extracted from the lattice data up to a multiplicative renormalization factor. For this project, we compute where i (t) is a spatial (temporal) Lorentz component andê is a unit polarization vector. We find it beneficial to simultaneously fit for the pseudoscalar to determine the ground state mass and amplitude, and for the vector. We find that over the whole temporal range (excluding t/a = 0) a multicosh/sinh fit with three states (N = 3) does a good job of describing our data. n = 1 is our lightest state, and so we can use the fit parameters in the following way to define the decay constants of eq. (3.5), For the renormalisation factors we take the results from [45] who determined them using 1-loop perturbation theory (see also [46]), with coefficients C A = 15.7, C V = 20.62, and C P = −6.71.
In a theory with more than one quark flavour, there would be a pseudo-Goldstone boson like the pion of QCD for which lattice calculations do not require disconnected contributions and so can give a very precise determination of the mass. In our single flavour theory, this state is absent. Nevertheless, we calculate this fictitious mass, which we denote as m π , by neglecting disconnected contributions to eq. (3.6). Even though m π is not a state in our theory, it provides a convenient alternative to the renormalized quark mass and allows us to determine the "chiral" point at which the quark mass vanishes. The extrapolation to m 2 π = 0 defines a critical value of bare quark mass m 0 that we will call m c . Numerically, we find m c = −0.9029(4). (3.10)

JHEP12(2018)118
Another way to determine m c uses a calculation of the topological susceptibility and leads to a compatible result, as shown in appendix B. We define the quark mass to be m q = m 0 −m c .

String tension and confinement scale
In lattice calculations, dimensionful quantities are given in units of the lattice spacing a and must be determined by fixing a physical scale. Since there are no known fixed scales to normalize our dark sector model, we will use the dark confinement scale Λ MS to define the overall scale of our theory. We will then report masses and decay constants in units of Λ MS . As in QCD, Λ MS is the characteristic energy scale of strong interactions. While Λ MS is purely defined in perturbation theory, we can perturbatively match it to the string tension σ, which can be directly measured in simulations. Following eqs. (4.60) and (4.61) of [47] and perturbative factors from [48], the result for our SU(2) theory with one fundamental quark flavour is The string tension is the slope of the linear potential between two color charges at large separation. The potential between a static quark and static anti-quark can be measured by tracing over Wilson loops connecting points x and x + r which have temporal length τ (see [47] and references therein), (3.12) Since generating these Wilson loops for each possible separation r/a is somewhat expensive, following [49], we fix the fields to Coulomb gauge and then need only compute open-ended Polyakov line correlators because the gauge condition will connect the ends of the line spatially. We start by measuring all matrix-valued Polyakov lines P , of length τ , from timeslice T , over L 3 , as Ut(x, t). (3.13) We can then directly compute the quantity for all separations r (and all of their translations over the L 3 volume) cheaply by performing the convolution with fast Fourier transforms, one direction is therefore L/2. We will average over equivalent r 2 values to further boost statistical precision. We can investigate where the static potential has saturated its ground state by looking at an "effective mass", The signal degrades for large values of τ and suffers from excited state contamination at very small τ , so an appropriate middle-ground must be taken. The static potential is often fit to the Cornell-type model [50], The dimensionless quantity that is extracted from the lattice simulation is a 2 σ. Based on our calculation of the massless limit in eq. (3.10), we fit simultaneously the mass dependence of eq. (3.17) with V (r) = B(1 + c 1 m q ) + a 2 σ(1 + c 2 m q )r. (3.18) For the string tension, we only need to fit the constant and linear terms of the Cornell potential. The fit parameters B and a 2 σ then give the potential in the massless-quark limit. We note that at small r there are significant discretisation effects, and at large r we expect significant signal deterioration and finite volume effects. Hence, we have performed our fits between these extremes, performing a fit-window analysis where we varied the upper and lower ends of the window looking for both a minimum in χ 2 and stability in the fit parameter a 2 σ. We then used a representative fit window to obtain our quoted results. Fits for our 12 3 × 32 lattices are displayed in the left panel of figure 1 and numerical results for both volumes are listed in table 1. Our final result for the string tension is 19) and that combines with eq. (3.11) to give aΛ MS = 0.249 (8) .  Table 1. Global fit results for the static potential with T = 32. We used τ = 2 in eq. (3.16) to determine aV (r).
This auxilliary scale allows us to quote all of our results in terms of the physical confinement scale that phenomenologists can choose, Λ MS , instead of the dimensionless quantities we directly compute. When high precision is required, lattice QCD studies typically set the scale with quantities called t 0 and w 0 rather than using the string tension. Appendix C presents our calculation of those quantities. However, for our exploratory study, Λ MS is convenient, perhaps more phenomenologically relevant, and entirely sufficient. Figure 2 shows the masses of some of our lightest hadrons from simulations with our lightest bare quark masses: the pseudoscalar η, vector ρ and the axial vector hadron a 1 . The η appears about a factor of 2 lighter than the ρ and a factor of 3 lighter than the a 1 for our lightest simulated quark mass. The a 1 is noticeably heavier than the ρ but approaches it in our massless-quark limit. Disconnected diagrams have been omitted from the axial vector calculation because they were found to be too noisy to make any quantifiable contribution.

Hadron masses and decay constants
All three hadrons' mass dependence for small quark masses can be described quite well by a linear fit in m q , as is illustrated in figure 2. However, the mass-dependence of the η appears to slightly prefer the form m 2 η = a + bx, giving a χ 2 /d.o.f ≈ 1.5. This is in line with the naïve expectation that the η receives a constant shift due to the anomaly even at vanishing quark mass [30].
It is worth noting that the determination of m η for m q /Λ MS = 0.233 is particularly low, and this is probably the reason for poor fits at larger quark mass. We have also noticed that the disconnected contribution to the η becomes more difficult to measure at larger quark masses.
The decay constants f η , f ρ , and f P of eq. (3.5) are displayed in figure 3. f ρ is approximately 3 4 Λ MS and f η is much smaller, ranging from 1 3 to 1 5 the size of the ρ decay constant over the handful of masses in table 2. The decay constant f P is of importance to DM phenomenology (c.f 5.2) and we find that its value is roughly twice the size of the η decay constant over the quark mass range considered here, and shows little sign of curvature.
We found that a simple linear fit describes the data of f P /Λ MS very well over our range of lightest quark masses and this is evident in the plot. However, some level of curvature appears present in both f ρ and f η . We believe this to be the onset of higherorder corrections of m q or some other functional dependence affecting our extrapolation and so find the data to be reasonably well described by the quadratic form f ρ/η = a+bm q +cm 2 q .  An overview of the broad mass spectrum in this minimal dark theory is given by figure 4. It appears as though this theory is somewhat reminiscent of QCD in the hierarchy of its spectrum; it has a light pseudoscalar meson, a heavier vector meson, and heavier still axial and scalar mesons. The decay constant for the pseudoscalar is smaller than that of the vector by about a factor 4, which in QCD is about a factor 2 different. Our dark matter candidate sits at roughly 2Λ MS , which is in the same ballpark as the ρ-meson in QCD. With dark matter phenomenology in mind, we note that m η < m ρ appears true for any value of the quark mass.

JHEP12(2018)118
The extension of figure 4 to larger m q , however, should be viewed with caution. Lattice artifacts can become large where am q > 1. Nevertheless, the right side of figure 4 shows a phenomenon familiar from heavy quark physics: hyperfine splittings shrink to produce a degeneracy of pseudoscalar with vector and also scalar with axial vector.

Direct detection and Higgs decay
The lightest vector meson is our DM candidate and is represented by an iso-triplet vector field ρ a µ , where a labels baryonic isospin. At lowest order, there are two operators, O S,P , that may couple to the Higgs field, as given in eq. (2.8). Since the CP phase φ is arbitrary, we may treat the corresponding mass scales M S = M/ cos φ and M P = M/ sin φ as separate parameters.

JHEP12(2018)118
Including the scalar operator, the low-energy effective Lagrangian for DM is where ρ a µν is the field strength tensor and λ S = ρ|qq|ρ /M S is a coupling determined below. 4 We have omitted purely dark sector interactions, e.g., with the η meson, that are beyond the scope of this work.
The low-energy theory of dark baryons in eq. (5.1) is reminiscent of models of hidden vector DM coupled via the Higgs portal [51][52][53][54][55]. In these models, one typically assumes that the Higgs interaction governs the DM relic abundance, implying a lower bound on λ S . This parameter space is strongly constrained by a combination of direct detection and Higgs decay limits [53]. In our framework, however, this assumption is not necessary since strong dynamics within the dark sector determine the relic density.
The spin-independent DM-nucleon cross section is [52] where m h is the Higgs boson mass and f N ≈ 0.3 is the Higgs-nucleon coupling [56]. The coupling λ S depends on a matrix element determined by our lattice results. The Feynman-Hellman theorem allows us to write where f S = ∂m ρ /∂m q . We determine f S from our lattice results using two methods. We found that an empirical fit of the form F (x) = (a 1 +a 2 x+a 3 x 2 )e −a 4 x +(a 5 +a 6 x)(1−e −a 7 x ) describes our data for the whole range of m ρ , and we use this fit to take the derivative. Second, we compute the derivative using the finite differences of the points. Both methods, shown in figure 5, are in good agreement and yield values in the range f S ≈ 1−3. However, lattice artifacts are present for m q /Λ MS 1 (corresponding to m ρ a −1 ), likely leading to an underestimate of f S . We expect f S ≈ 2 at large quark mass since m ρ ≈ 2m q . We additionally caution using the empirical fit results to extrapolate to m q = 0 as this fit has little predictive power; likewise extrapolating the low-mass finite differences lacks significance. If this point is of interest then we suggest using the slope of the light mass extrapolation from table 2, which is 2.1.
Direct detection limits are most constraining for weak-scale DM mass. Recently, XENON1T obtained the most stringent upper bound on the spin-independent cross section, 4.1 × 10 −47 cm 2 for 30 GeV DM mass [57], implying M S > 28 TeV. For larger DM mass, the XENON1T bound weakens while σ ρN is nearly constant.
Higgs studies at the LHC provide the most stringent constraints for low mass DM. In our model, the Higgs boson may decay into dark sector states that are long-lived and escape 4 We expect the pseudoscalar term in eq. (2.8) to induce a Higgs-DM interaction of the form ε αβµν ρ a αβ ρ a µν (|H| 2 − 1 2 v 2 )/MP . However, this leads to a velocity-suppressed direct detection cross section that is much less constrained compared to the scalar interaction. the detector. If we assume m q , Λ MS m h /2 and that all dark states escape invisibly, it is straightforward to compute the Higgs invisible width from a quark-level calculation. We have which is independent of m ρ , the CP phase φ, or any other dark sector parameters. Present limits constrain the Higgs invisible branching fraction to be below 23% [58,59]. For our model, this implies M > 40 TeV. We note that the invisible Higgs constraints are very different compared to hidden vector DM models where DM is a gauge boson, not a composite state. In that case, the Higgs invisible width scales as Γ(h → inv) ∝ m −4 DM [54] and becomes very constraining for light DM (see, e.g., figure 9 of [58]).

Fate of the lightest dark hadron
The lightest state in the dark spectrum is the η meson. If it were stable, it would constitute an O(1) fraction of the DM density. However, the η meson is not a worthy DM candidate since it can mix with the Higgs boson through a dimension-five operator O P |H| 2 , inducing it to decay to the SM. Even if this operator is suppressed by the Planck scale, the η lifetime would be much shorter than the age of the Universe. Its decay products, moreover, are fixed by the SM Higgs couplings.
The lifetimes of meta-stable dark states are strongly constrained if they decay into visible SM particles. Cosmic microwave background measurements exclude an O(1) fraction of meta-stable DM unless it decays prior to recombination, before ∼ 10 13 s [60]. Decays occurring between ∼ 0.1 − 10 12 s affect primodial abundances of light nuclei [61]. In particular, decays via Higgs mixing are largely constrained to occur before ∼ 0.1 s, otherwise the injection of hadrons into the plasma alters the neutron/proton ratio after weak interactions have frozen out [62]. However, the limits depend on the cosmological abundance of η mesons before they decay, which we defer to future work. Here, to be conservative, we require the lifetime to be τ η < 1 s. The total η width can be written as The first term represents η decays through Higgs mixing, where the mixing angle θ hη is defined by and Γ h is the total SM Higgs width (evaluated at m η , not m h ). We have adapted results from ref. [63] to get Γ h as a function of mass below bottom threshold, while for larger mass we take results from ref. [64]. The second term in eq. (5.5) is an additional decay channel that opens for m η > 2m h . The combination of τ η < 1 s and invisible Higgs decay yields a lower limit m η > 228 MeV and m ρ > 320 MeV for the range of m q in table 2. This conclusion is further bolstered by astrophysical constraints on self-interactions, discussed below. Figure 6 illustrates the complementarity between different constraints. For definiteness, we have taken m q /Λ MS = 0.1 and CP phase φ = π/4. The remaining parameters of the model are the DM mass m ρ and the interaction scale M . Other parameters of the model are determined according to our lattice results: m η ≈ 0.57 m ρ , f P ≈ 0.39, and f S ≈ 1. We have truncated the invisible Higgs limits at 10 GeV since the assumptions leading to eq. (5.4) eventually breakdown. With the exception of tuning φ = 0, taking other parameter choices does not greatly shift the shaded regions.

Self-interactions
In our model, DM particles are not collisionless and elastically scatter with one another through strong interactions. If the scattering rate is large enough, self-interactions can leave an observable imprint on DM halos of galaxies and clusters. The relevant figure of merit is σ el /m, the cross section for DM elastic scattering per unit DM mass, which is typically expressed in units of cm 2 /g ≈ 2 barn/GeV. While self-interacting DM is often motivated in terms of explaining various small scale structure issues [5], here we simply make a conservative constraint on the parameter space of our model. Actually calculating σ el /m is a challenging prospect for the lattice that we defer to future work.
By dimensional analysis, we expect σ el ∼ 4πΛ −2 MS since Λ MS sets the typical size of ρ. Since m ρ > 2Λ MS for any dark quark mass, we can therefore set a lower bound Observations of relaxed massive clusters [65,66] provide the strongest constraint on selfinteractions, favoring σ el /m ≈ 0.1 cm 2 /g or less [67]. If we take σ el /m < 0.5 cm 2 /g as a conservative upper limit [68], we have Merging cluster constraints, such as the Bullet Cluster [69], are comparatively weaker. In particular, recent simulations have found offsets for self-interacting DM halos to be much smaller than previously thought [70]. Our dimensional analysis estimate breaks down if DM scattering has an s-wave resonance, corresponding to a di-baryon (ρρ) that is a nearly zero energy bound state. In this case, σ el /m can be far larger than the lower bound implied by eq. (5.7), approaching the s-wave unitarity limit when the mass gap and scattering energy go to zero [71]. This is analogous to proton-neutron scattering, which is enhanced owing to the smallness of the deuteron binding energy. Eq. (5.8) is still satisfied in this case. On the other hand, antiresonances (the Ramsauer-Townsend effect) may act to suppress DM scattering for certain choices of parameters [72,73], evading our limit, but without a detailed calculation it is not possible to say anything further.

Conclusions
Since strong dynamics explains the mass and stability of visible baryons, it is possible that similar physics is realized for DM as well. In this work, we have studied the simplest model of dark baryons: SU(2) gauge theory with one flavor of dark quark. Unlike QCD, the theory has no spontaneously broken chiral symmetries and no pseudo-Goldstone bosons. Instead, there is an unbroken global SU(2) B baryon symmetry resembling isospin, which unifies baryons and mesons into degenerate iso-multiplets. The lightest baryon is one component of a iso-triplet vector ρ, which is our DM candidate. Dark hadrons may couple to the SM through non-renormalizable interactions and we have considered the leading dimension-five operators involving the Higgs field.

JHEP12(2018)118
In this initial and exploratory study, we have used lattice simulations to compute the spectrum of the lightest dark hadrons. The overall mass scale of the theory is unknown a priori. Hence, with an eye towards phenomenology, we have presented all dimensionful parameters normalized with respect to the confinement scale Λ MS (computed from the string tension σ). The dark quark mass m q is a free parameter and our simulations focus on the quark mass regime with m q /Λ MS ≈ 0.1 → 1. In this range, the lightest hadron is the iso-singlet pseudoscalar meson η. We have included the effect of disconnected diagrams, which causes the η to remain massive according to our extrapolation to m q = 0, as expected from the U(1) A anomaly. The iso-triplet vector ρ is the next-to-lightest state. We have also presented results for the lightest axial vector and scalar, which remain heavier still.
We note that several sources of systematic error have not been accounted for. As our volume is quite small, we expect significant finite volume effects, particularly for light quark masses. We also expect finite lattice-spacing artifacts to be present since our lattice spacing is somewhat coarse and our action is correct only up to O(a) discretisation effects. However, for a first study, the broad brush strokes of this theory are what is important and we anticipate our results to be accurate at around the 10% level with these systematics in mind. Now that we better understand the parameter space and the model's feasibility as a DM candidate, dedicated finite volume and continuum limit studies beyond fixed L and β will be necessary to refine our numerical predictions.
In our opinion, there are three nice features of our model worth re-emphasizing, apart from its minimality.
• DM stability: the accidental SU(2) B baryon number symmetry is preserved up through operators of dimension-five. From an effective theory point of view, our DM candidate is as stable as the proton (and a counterexample to arguments in ref. [20]).
• CP violation and η decay: including dimension-five operators, the dark quark receives a mass contribution from the Higgs field in addition to its bare mass. Since both terms need not be aligned in general, there appears a CP phase that mixes the η with the Higgs boson, allowing the η to decay rapidly in the early Universe before nucleosynthesis.
• Annihilation channel: our model has a built-in mechanism for efficient annihilation to set the DM relic density, ρρ → ηη, with the η mesons later decaying to the SM. Since our lattice results show that m ρ > m η for any quark mass, this process is always kinematically allowed.
On the phenomenology side, we have arrived at the following conclusions. There is a lower limit on m ρ , m η of a few hundred MeV from combining Higgs invisible decay constraints with bounds on the η lifetime from nucleosynthesis. A similar limit, m ρ > 280 MeV, is required from constraints on DM self-interactions in clusters. For larger DM masses, the parameter space is constrained by Higgs invisible decays and direct detection, implying that the scale M connecting the dark sector with the Higgs field must be larger than 1 − 40 TeV depending on m ρ . We have used our lattice results to extract the η decay constant and DM scalar form factor needed for these calculations. At the same time, other possibilities remain for coupling our SU(2) theory to the SM (e.g. through a Z ), which will change many of these conclusions.

JHEP12(2018)118
. We select where m 2 π → 0 to be the point of our vanishing quark mass. This is consistent with another method to define this point through the topological susceptibility, described below. Figure 7 illustrates that the physical η is approximately a constant shift above the π in mass. Although this constant shift appears to be fairly small, indicate that the η remains massive in the "chiral" limit.

B Topological susceptibility
In lattice QFT, topological charge Q can be defined from the gauge fields and the topological susceptibility χ can then be obtained from Our calculation of F µν (x) is the average of all four plaquettes in the µ − ν plane that touch the point x, the standard clover definition. An important issue to note with this discretisation is that the lattice values for Q do not tend to be integers, due to shortdistance effects which must be reduced by some smoothing procedure. For this smoothing we will use HYP smearing [75], monitoring the stability of Q 2 as the number of smearing iterations is increased. We can expect for N f light quark flavors that [76] χ = Σ where Σ is the chiral condensate. This implies that the limit χ → 0 occurs when the quark mass vanishes. We will measure the topological susceptibility by the "slab method" [77,78], computed on sub-volumes V = L 3 ∆, For 0 < ∆ < T , the translationally-invariant sum is best performed using convolutions over the slab. The left panel of figure 8 shows our numerical determination of the topological susceptibility. The right panel illustrates the improvement of the slab method relative to the standard method, which is simply the slab method with ∆ = T . If we fit the slab method determinations only up to L/2 we find good ( 1.5× reduction in error) statistical improvement over using the full volume determination. This indicates that the full-volume sum is  noisy, and a truncated sum over a sub-volume contains less noise but still captures the relevant physics. We observe stable results after approximately 21 HYP smearing iterations. Figure 9 confirms that our simulations are not getting stuck in a particular topological sector. We find that the integrated autocorrelation time for the topological charge is less than our chosen spacing for measurements in Monte-Carlo time.
The calculation of the topological susceptibility for a range of bare quark masses permits an extrapolation to zero as shown in the left panel of figure 10, representing the limit of a massless quark for that lattice volume. We have a few lighter quark masses here compared to those listed in table 3. These however were very difficult to invert for the meson spectrum and we suspect that they contribute large finite volume systematics to hadronic measurements. However, for this noisy gauge-field quantity they seem to be acceptable to use. Repeating this procedure on a second lattice size allows an extrapolation to the limit T → ∞ and our result is plotted in the right panel of figure 10. which is in good agreement with eq. (3.10). This consistency from two different methods is reassuring. We will use eq. (3.10) to define the massless limit since it has a slightly smaller error bar.
C The lattice scales t 0 and w 0 Our results have primarily used √ σ to set the physical scale of this dark matter theory, due to its direct phenomenological interpretation. In lattice QCD calculations, however, it has become common to invoke standardized parameters named t 0 and w 0 because they can be determined much more precisely than the string tension. We report our calculations of these quantities here to facilitate comparison with future lattice studies of this theory.
To begin, we generalize the gauge link U µ (x) → U µ (x, t) where t represents the flow time. The original, un-flowed link value is obtained at t = 0. The flow time does not have units of physical time, and the dimensionless quantity that emerges from a lattice simulation is a 2 t. Gradient flow is defined by dU dt = Z(U )U .

(C.1)
U is shorthand for the gauge field at a particular flow time and Z(U ) is chosen to be the "force term" which is essentially the factor within the lattice action that multiplies this particular link. The equation is solved by performing an iterated flow with (small) step size , We can use this technique to very accurately define a scale through [79,80] G(t) = t 2 F µν F µν , G(t 0 ) = N/10. The lattice spacing derived from these two definitions should be consistent up to discretisation effects. The factor of N originates from the correct identification of the t'Hooft limit in comparison to the commonly used value of 0.3 for SU(3) [82,83].
The left panel of figure 11 shows our numerical results for one ensemble, and the right panel shows the fit to quark mass and the massless limit, leading to √ t 0 a = 1.357 (7), w 0 a = 1.416 (10).  Table 3. The bare mass m 0 and the number of configurations generated N conf for the ensembles used in this work.