Stealth dark matter and gravitational waves

In recent years there has been increasing interest in the possibility of searching for a stochastic background of gravitational waves produced by a first-order phase transition in the early universe. This involves future space-based gravitational-wave observatories as discussed by Refs. [1, 2, 3] and references therein. Here we report initial investigations into the extent to which these gravitational-wave searches may constrain or discover the stealth dark matter (DM) model introduced and investigated through non-perturbative lattice calculations in Refs. [4, 5].

Stealth DM involves a composite bosonic ‘dark baryon’ DM candidate (analogous to the neutron) that arises from an SU(4) gauge theory coupled to four flavors of moderately heavy fermions in the fundamental representation (i.e., with masses roughly comparable to the confinement scale). In order to generate the correct cosmological abundance of DM, these fermions transform in non-trivial (vector-like) representations of the electroweak group, in such a way that the lightest (spin-0) dark baryon is a singlet under the entire standard model (SM) gauge group. This DM candidate is automatically stable on cosmological time scales due to the conservation of dark baryon number, in parallel to proton stability in the SM. Its mass arises from confinement and from the technically natural masses of its ‘dark fermion’ constituents, with natural mass scales of order $1\lesssim M_{DM}\lesssim 100$ TeV.

Even though those ‘dark fermions’ are electrically charged and couple to the Higgs boson, the symmetries of the stealth DM model dramatically suppress direct detection cross sections. Direct detection searches only mildly constrain the dark baryon’s effective Higgs interaction [4, 6], while the leading (dimension-5) magnetic moment and (dimension-6) charge radius effective operators governing photon exchange are both forbidden. The cross section due to the dimension-7 electromagnetic polarizability interaction then provides an unavoidable lower bound on the entire class of DM models involving dark baryons with charged constituents (reviewed in Ref. [7]). The resulting direct detection constraint $M_{DM}\gtrsim 0.2$ TeV [5] is comparable to the collider constraint $M_{DM}\gtrsim 0.3$ TeV [4, 8] set by searches for charged ‘dark pions’. These bounds are roughly two orders of magnitude weaker than those for direct detection of an SU(3) model with unsuppressed magnetic moment and charge radius interactions [9].

The connection to gravitational waves comes from the possibility that the confinement transition in the SU(4) sector could be first order. For SU($N$) gauge theories with $N\geq 3$ and a small number of fundamental fermions, such first-order transitions occur in two regimes: where the fermions are all sufficiently heavy or all sufficiently light. This is illustrated in the ‘Columbia plot’ sketched in Fig. 1. In the context of stealth DM we are interested in the region of heavier fermion masses, connected to the pure-gauge theory recovered in the ‘quenched’ limit where the fermions become infinitely massive. (In the chiral limit the baryon-to-pion mass ratio diverges, allowing collider searches to rule out the entire parameter space.)

Our investigations proceed in two stages. First we need to determine roughly how heavy the dark fermions must be in order for the stealth DM confinement transition to be first order and produce a stochastic background of gravitational waves. If the dark fermion mass must be much larger than the confinement scale, then dark glueballs may be stable and contribute to the relic density, requiring reconsideration of the phenomenology and constraints reported by Refs. [4, 5].

The second stage, after locating a first-order transition, is to study it in more detail in order to predict the spectrum of gravitational waves it would produce. The key parameters that need to be computed or estimated for this purpose are the latent heat (or vacuum energy fraction), the phase transition duration, and the bubble wall velocity [1, 2, 3]. Only the first of these is straightforward to determine through lattice calculations.

Completing this work is necessary to translate future searches for stochastic gravitational waves (resulting in either detections or exclusions) into novel constraints on stealth DM and related models. For example, the gravitational-wave spectrum also depends on the transition temperature $T_{*}$ (which may differ from the equilibrium critical temperature $T_{c}$ used to set the scale of our lattice calculations, due to possible supercooling). If we can assume $T_{*}\simeq T_{c}$ or estimate how they differ, observations of gravitational waves will indicate both the approximate mass scale of the dark baryons as well as the minimum mass of the dark pions being searched for at colliders.

This proceedings focuses on the first step described above, investigating the nature of the SU(4) confinement transition as we vary the mass of the four dynamical fermions. We first consider the limiting case of SU(4) pure-gauge theory in the next section, confirming that we can observe its known first-order thermal transition. In Section 2 we then couple the theory to four degenerate fundamental fermions, and observe a clear change in the qualitative features of the transition upon decreasing the fermion mass $\href{http://www-spires.slac.stanford.edu/spires/find/hep/www?j=ANMAA\%2C$\%2CW}{{\it Ann.\ Math.\ }{\bf$} (.) W}econcludeinSection~{}\ref{sec:conc}bydiscussingtheremainingworkneededtopredictthegravitational-wavespectrumofstealthDM.\par\par\par\vspace{-8 pt}$

There have been many lattice investigations of the confinement transition of SU($N$) Yang–Mills theory over the years. See Refs. [10, 11] for work with a focus on $N=4$ and Ref. [12] for a broader review. Our main goals in revisiting this calculation are first to confirm that our code and algorithms are working correctly, and then to estimate appropriate lattice volumes $L^{3}\!\times\!N_{t}$ to use for the more expensive calculations with dynamical fermions. For both pure-gauge and dynamical configuration generation with the HMC algorithm we use QHMC/FUEL [13], which is built on the USQCD SciDAC software stack and provides efficient performance for arbitrary SU($N$).

Even though each individual lattice ensemble generated for this project is modest in size and computational expense, many ensembles are needed, so it is worthwhile to consider the smallest viable $L$ and $N_{t}$. For each fermion mass and the $\href{http://www-spires.slac.stanford.edu/spires/find/hep/www?j=ANMAA\%2C\to\%2C$}{{\it Ann.\ Math.\ }{\bf\to} (\infty) $}pure-gaugelimit,wewantatleastthree$N_t$inordertoenablecontinuumextrapolations,ineachcasewithatleastthreeaspectratios$L / N_t$inordertoenableinfinite-volumeextrapolations.Fixing${Ann. Math. , (N) _t, L / N_t}$,wescanbothfromhightolowtemperaturesandfromlowtohightemperatures,tocheckforpossiblehysteresis.Intotal,withfourfinite$Ann. Math. = () {0.05, 0.1, 0.2, 0.4}$, four $N_t = {4, 6, 8, 12}$ and five $L / N_t = {2, 3, 4, 6, 8}$ we have generated 1,310 ensembles each with 5k–50k molecular dynamics time units (MDTU). We monitor the autocorrelation times $τ$ for (non-topological) observables of interest, finding that the longest of these can reach $τ≈2$k~MDTU in the vicinity of the transition. We accumulate enough data to ensure at least 14 statistically independent measurements after thermalization.

It is well known that severe discretization artifacts are likely to be present for small $N_t ≤4$, related to the thermal confinement transition merging with a bulk transition into a lattice phase~[14, 15]. Following Ref. [16] we have attempted to avoid this bulk transition by using a gauge action that includes both fundamental and adjoint plaquette terms, with negative adjoint coupling $$\beta$ _A = -$\beta$ _F / 4$. This differs from the $$\beta$ _A = 0$ case considered in Refs.~[10, 11], requiring continuum extrapolations in order to carry out quantitative comparisons with that past work. From the left panel of Fig. 2 we conclude that it is still unsafe to rely on $N_t = 4$ results in continuum extrapolations, since $N_t ≥6$ is needed to clearly separate the $N_t$-dependent thermal confinement transition (signalled by the Wilson-flowed Polyakov loop susceptibility $χ_PL_W$) from the $N_t$-independent bulk transition (signalled by the plaquette susceptibility $χ_□$). See Ref. [17] for background on our use of the Wilson flow to improve signals for the Polyakov loop and its susceptibility; we carry out measurements with $8t / N_t = 0.5$, corresponding to flow time $t = 2$ for $N_t = 8$.

In the right plot of Fig. 2 we confirm that the confinement transition we see is first order, by measuring a double-peaked histogram for $|PL_W|$ on $24^3$\!\times\!$ 8$ lattices with $$\beta$ _F = 15$. The lattice volumes we have generated do not show any hysteresis, and the statistics we have accumulated do not produce clear volume dependence either in the height of the Polyakov loop susceptibility peak or in the Polyakov loop kurtosis~[18]. The presence of two peaks in the $|PL_W|$ histogram is therefore the main evidence we use to distinguish between first-order transitions and crossovers.

Having completed the exercise of mapping out the finite-temperature phase diagram for SU(4) Yang–Mills theory, we now repeat this work including dynamical fermions with a range of bare masses $Ann. Math. $ (.) Although stealth DM requires at least a small splitting between two pairs of degenerate fermions (to guarantee that charged `dark mesons' decay before Big Bang nucleosynthesis~[4]), for simplicity we carry out these finite-temperature studies with four degenerate flavors, corresponding to the ``$N_f = 4$'' diagonal line in Fig. 1. This allows us to use unrooted staggered fermions, and we continue to follow Ref. [16] by working with an improved lattice action featuring nHYP smearing with parameters $(0.5, 0.5, 0.4)$. Since continuum extrapolations remain work in progress, in this proceedings we focus on the largest $N_t = 8$ for which a large amount of data is available.

We parameterize the fermion masses by computing the corresponding ratio of dark pion and dark vector meson masses, $M_P / M_V$, which involves additional calculations of the meson spectrum on zero-temperature $24^3$\!\times\!$ 48$ lattices at the critical $$\beta$ _F^(c)$ for the given $Ann. Math. $ (.) This simplifies comparisons with previous quenched lattice studies of stealth DM~[4, 5], which used valence Wilson fermions with masses corresponding to $0.55 $\lesssim$ M_P / M_V $\lesssim$ 0.77$. The four staggered $Ann. Math. = () {0.05, 0.1, 0.2, 0.4}$ we consider correspond to $M_P / M_V = {0.65(3), 0.80(3), 0.91(1), 0.96(1)}$, respectively, with uncertainties set by varying $$\Delta$ $\beta$ _F = ±0.2$ around $$\beta$ _F^(c)$. We chose the smallest $Ann. Math. = (0) .05$ in order to overlap the mass range considered in Refs.~[4, 5], while the larger masses are needed to obtain a clearly first-order transition.

The left plot of Fig. 3 shows the real part of the Wilson-flowed Polyakov loop $PL_W$ vs. $$\beta$ _F$ to illustrate how the dynamical fermions shift the confinement transition to stronger critical $$\beta$ _F^(c) ≈{12.7, 13.1, 13.6, 14.2}$ as Ann. Math. d (e) creases. Even though the heaviest fermion mass $Ann. Math. = (0) .4$ produces a rather large mass ratio $M_P / M_V ≈0.96$, the corresponding $$\beta$ _F^(c) ≈14.2$ is still significantly different than the pure-gauge $$\beta$ _F^(c) ≈15.0$ for our fundamental–adjoint lattice action. Even at this heaviest mass, the dynamical fermions have significant effects. So far it is only for this $Ann. Math. = (0) .4$ that we observe two-peak $|PL_W|$ histograms like the one shown in the right plot of Fig. 3, which establish that this mass is sufficiently heavy to remain in the upper-right region of Fig. 1 where the transition is first order and can produce gravitational waves.

As mentioned at the end of the previous section, the statistics that we are able to accumulate have not yet allowed us to clearly distinguish between first-order transitions and crossovers using the most traditional observables such as the volume dependence in the height of the Polyakov loop susceptibility peak or in the Polyakov loop kurtosis. A different observable that currently appears more promising is the `deconfinement fraction' discussed in Refs.~[10, 19], which measures the proportion of $arg(PL_W)$ measurements that fall within a certain (tunable) angle $θ$ around any of the $Z_4$ symmetry axes. Although this quantity was originally developed in the context of pure-gauge theories, it remains well-defined in the presence of fundamental fermions that favor the positive real axis. With $N_in$ of $N_tot$ measurements suitably aligned along the $Z_4$ axes, defining $f ≡π/ 4π/ 4 - θ[N_inN_tot - θπ/ 4]$ shifts and normalizes the deconfinement fraction so that $f →1$ in the deconfined phase and $f →0$ in the confined phase. This is shown in Fig. 4, where the key feature is the contrast between the results for $Ann. Math. = (0) .4$ (which behave the same as the first-order transition of the pure-gauge theory) and those for $Ann. Math. = (0) .2$, which show qualitatively different dependence on the spatial lattice volume. This volume dependence becomes even more pronounced as Ann. Math. d (e) creases further, consistent with the transition changing to a crossover between $0.2 < Ann. Math. < (0) .4$.

Stealth dark matter features an early-universe confinement transition that—if it is first order—can produce a stochastic background of gravitational waves that will be searched for using future space-based observatories. To investigate this possibility, the Lattice Strong Dynamics Collaboration is analyzing the finite-temperature dynamics of stealth DM, first determining the range of dark fermion masses for which the transition is first order and then studying this first-order transition in more detail to predict the characteristics of the gravitational waves it would produce. This proceedings focused on the first step, finding that a crossover persists up to large masses corresponding to $M_P / M_V ≃0.9$, with evidence for a first-order transition when $M_P / M_V ≈0.96$.

Work on the second step is underway, focused on determining the latent heat of the transition, and its extrapolation to the infinite-volume continuum limit. While it will be more challenging to establish robust non-perturbative constraints on the phase transition duration and bubble wall velocity for the first-order stealth DM transition, these lattice calculations should be able to provide new insight into those quantities. It will also be worthwhile to explore alternative means of analyzing such first-order phase transitions, such as the density of states approach~[20].