Progress applying density of states for gravitational waves

Many models of composite dark matter feature a first-order confinement transition in the early Universe, which would produce a stochastic background of gravitational waves that will be searched for by future gravitational-wave observatories. We present work in progress using lattice field theory to predict the properties of such first-order transitions. Targeting SU(N) Yang--Mills theories, this work employs the Logarithmic Linear Relaxation (LLR) density of states algorithm to avoid super-critical slowing down at the transition.


Introduction
Standard Markov-chain Monte Carlo importance sampling techniques have been the workhorse for many applications in modern theoretical physics, yet they have some drastic intrinsic shortcomings for certain situations.These include lattice field theory studies of first-order phase transitions, where super-critical slowing down characterized by uncontrollable autocorrelations can result from the difficulty of tunnelling between the two coexisting phases on large lattice volumes [1,2].In such situations alternatives, such as density of states approaches, may perform better.
First-order phase transitions are currently attracting a great deal of interest due to the possibility that such phenomena in the early Universe could produce an observable background of gravitational waves -see Ref. [3] and references therein.In this work we focus on gravitational waves from first-order confinement transitions in composite dark matter models, which has also been considered by Refs.[4][5][6].In particular, we consider the Stealth Dark Matter model proposed by the Lattice Strong Dynamics Collaboration [4,7,8].This is an SU(4) gauge theory coupled to four fermions that transform in the fundamental representation of the gauge group.While these fundamental fermions are electrically charged, they confine to produce a spin-zero, electroweak-singlet 'dark baryon'.
In addition to guaranteeing the stability of a massive dark matter candidate through an analogue of baryon number conservation, the symmetries of Stealth Dark Matter strongly suppress its scattering cross section in direct-detection experiments [7,8], especially for heavy dark matter masses M DM ≳ 1 TeV.Collider searches for Stealth Dark Matter also become challenging for such heavy masses [9,10], which helps to motivate the ongoing work [4][5][6] investigating how models of this sort could be constrained or discovered by future gravitational-wave observatories such as LISA [3].
In this proceedings we summarize the progress of our ongoing work applying the Logarithmic Linear Relaxation (LLR) density of states algorithm [11,12] to analyze SU(4) Yang-Mills theory.This pure-gauge theory can be considered the 'quenched' limit of Stealth Dark Matter corresponding to infinitely heavy fermions.It is a convenient context in which to consider the LLR algorithm, which is challenging to apply to systems with dynamical fermions [13].The SU(N) Yang-Mills confinement transition is first order for N ≥ 3, with significantly stronger transitions for N ≥ 4 [14].This makes the SU(4) theory a promising first target for application of the LLR algorithm, which could form the basis for future studies of either the SU(3) case relevant for QCD, or of N ≥ 5 to explore the large-N limit.
In the next section we begin by summarizing the LLR algorithm.We then discuss in Section 3 how we tested our LLR code by reproducing some results from Ref. [12] for compact U(1) lattice gauge theory.In Section 4 we present current results from our ongoing LLR analyses of the SU(4) Yang-Mills confinement transition, updating Ref. [15].We wrap up in Section 5 with a brief overview of our planned next steps.

Linear Logarithmic Relaxation
Observables of SU(N) Yang-Mills theories on the lattice are given by where S [ϕ] is the lattice action and ϕ is the set of field variables attached to each link in the lattice.Analytically evaluating these extremely high-dimensional path integrals is practically impossible for the systems we're interested in.Standard Monte Carlo techniques instead aim to give an approximation, by only using a modest number of configurations ϕ, sampled with probability ∝ e −S [ϕ] .An alternative approach is to calculate the density of states and use it to reconstruct the observables O in Eq. 1 by performing a simple one-dimensional integration over the energy, This quantity ρ(E) is not easily accessible.To calculate it we use the LLR algorithm [11,12], which begins by defining the reweighted expectation value Here E i is a fixed energy value, θ E i ,δ is the modified Heaviside function (1 in the interval E i ± δ 2 and 0 everywhere else), and for now 'a' is just a free parameter not to be confused with the lattice spacing.The next step is to set the reweighted expectation value ⟨⟨E − E i ⟩⟩ δ (a) to zero and approximate the integral with the trapezium rule: After expanding the exponential e ±a δ 2 and the density ρ(E i ± δ/2) in a Taylor series and neglecting O(δ 2 ) terms in the limit δ → 0, we arrive at the following expression for the parameter a: We can see from this, that a(E i ) is a linear approximation of the derivative of the logarithm of the density of states ρ(E i ).The full density of states ρ(E) can now be reconstructed, with an exponentially suppressed error [11,12,16], by performing a numerical integration of the values of a(E i ) across many intervals E i ± δ 2 and exponentiating the result.To compute a for a given E i , we use the Robbins-Monro algorithm to solve the equation This series has a fixed point at the correct value of the LLR parameter a = a (n+1) = a (n) .At each iteration j we evaluate the reweighted expectation value ⟨⟨E − E i ⟩⟩ δ (a ( j) ) using standard importance-sampling Monte Carlo techniques, but with the probability weight e −a ( j) S rather than the usual e −S .The modified Heaviside function θ E i ,δ in Eq. 4 means we reject all Monte Carlo updates that produce a configuration with an energy outside of E i ± δ 2 , potentially causing low acceptance rates for small energy intervals δ.Alternatively we can substitute this hard energy cut-off with a smooth Gaussian window function [13,16]: This effectively changes the probability weight in the Monte Carlo simulation to exp e −aE .Unlike the modified Heaviside function, the Gaussian window function is differentiable, enabling the use of the hybrid Monte Carlo (HMC) algorithm in the evaluation of the reweighted expectation value.In our work we test and compare both ways of constraining the energy interval.

U(1) bulk phase transition
As an initial test of our implementation of the LLR algorithm, we reproduced some results from Ref. [12] for four-dimensional pure-gauge U(1) theory, defined by the action ,   with is attached to the lattice site x in direction μ and the sum runs over each link in the lattice.β = 1 g 2 0 with g 2 0 being the bare gauge coupling.For this simple system, we use the Metropolis-Rosenbluth-Teller algorithm with randomly updated angles θ µ (x), imposing a hard energy cut-off (Eq.4) to evaluate the reweighted expectation value ⟨⟨E − E i ⟩⟩ δ (a (n) ).As shown in Ref. [15], our results for the values of a are consistent with the results from Ref. [12], providing confidence in our underlying implementation of the LLR algorithm.Figure 1 shows our results for a, using two different values of the energy interval size, δ = 0.01/V and 0.001/V.While we are interested in the limit δ → 0 from Eq. 8, we can see that smaller δ results in larger statistical uncertainties.In both plots the error bars are determined through a jackknife analysis of N j = 10 completely independent calculations for each energy interval.From these results for a we can reconstruct the probability density P β (E) = ρ(E)e βE via a simple trapezium-rule numerical integration.Figure 2 illustrates the sensitivity of this probability density to the value of β.While P β (E) has only a single peak for β = 1.02, it has a clear double-peak structure for β c = 1.00735, signalling the coexistence of the two phases at the first-order bulk phase transition.The latent heat ∆E/(6V) ≈ 0.045 of this first-order phase transition can be directly read off as the separation between the two peaks.

SU(4) phase transition
We now turn to our ongoing work applying the LLR algorithm to SU(4) Yang-Mills theory.The action is with the plaquette , with g 2 0 the bare gauge coupling, the sum runs over all lattice sites and U µ (x) is the SU(4)-valued link variable attached to lattice site x in direction μ.
We implemented the LLR algorithm in a fork of Stefano Piemonte's LeonardYM software. 1For the calculation of the reweighted expectation value ⟨⟨E − E i ⟩⟩ δ (a (n) ), we have compared three different updating schemes: overrelaxation updates in the full SU(N) group [17], the Metropolis-Rosenbluth-Teller algorithm with SU(N) updates generalized from Ref. [18], and HMC updates.While our use of the HMC algorithm requires the smooth Gaussian window function Eq. 10, in the other two cases we have further compared both this Gaussian window approach and the hard energy cut-off of Eq. 4. We obtain consistent results from all of these updating schemes.
For SU(N) Yang-Mills theories on an N 3 s × N t lattice, there are two distinct transitions we could consider.The physically relevant confinement transition is first order for N ≥ 3 (but only weakly so for N = 3) [14].This transition occurs at a fixed temperature T c = 1/(aN t ) (with this a the lattice spacing), implying β c → ∞ as N t → ∞.In addition, there is a bulk transition at an N t -independent coupling β bulk , which is first order for N ≥ 5 (but only weakly so for N = 5) [14].If N t is too small, β c ≈ β bulk and the confinement transition is distorted by the nearby bulk transition -even for N ≤ 4 where the latter is a continuous crossover.Based on prior work including Refs.[4,19], we consider N t ≥ 6 in order to avoid this problem.
As a brief aside on the bulk transition, in Fig. 3 we show results for the LLR parameter a across a wide range of energies for 6 4 lattices.Although we can see that the decrease of a decelerates around the energy of the bulk crossover, the function a(E) remains monotonic in contrast to Fig. 1 for the U(1) case.As a result the reconstructed probability density P β (E) does not feature a double-peak structure for any value of β, confirming a continuous crossover with no phase coexistence.In ongoing work soon to be reported, we will directly contrast this with the first-order bulk transition of SU(6) Yang-Mills theory.
To investigate the physically interesting first-order confinement transition of pure-gauge SU(4), we have to use a lattice with an aspect ratio r ≡ N s /N t > 1.While we have carried out LLR calculations using a range of lattice volumes, here we focus on V = 12 3 × 6, also fixing the energy interval size δ = 0.001/V and N j = 5 independent runs of the Robbins-Monro algorithm for each energy interval.Based on N t = 6 importance sampling results in Ref. [19], in Fig. 4 we focus on the energy range 13.2 < E/V < 13.9.There is no sign of a non-monotonic a(E), which continues to be the case throughout larger ranges of E/V, implying that we are not yet able to resolve the first-order confinement transition.We have confirmed this by reconstructing the probability density P β (E) using both naive trapeziumrule integration and a more robust polynomial fit technique [20,21].As expected, neither method produces a double-peak structure for any β.
We can gain some appreciation for why we have not yet observed a first-order SU(4) confinement transition by noting that the N t = 6 importance-sampling results in Ref. [19] imply a latent heat of ∆E/V ≈ 0.004 at β c ≈ 10.8.This latent heat is roughly two orders of magnitude smaller than that in the U(1) case, so it is not necessarily surprising that our LLR calculations have not yet resolved it.To test this hypothesis, we generated synthetic values of a that would correspond to a first-order phase transition with a similar pseudo-critical β c ≈ 10.8 and latent heat ∆E/V ≈ 0.004.In Fig. 5 we overlay our numerical data on top of this synthetic a(E) (blue line), which allows us to conclude that our statistical precision is currently insufficient to resolve the non-monotonocity corresponding to such a first-order phase transition.These considerations provide some clear next steps for our work, which we discuss in the next section.

Outlook and next steps
In this proceedings we have presented our progress applying the LLR density of states algorithm to investigate first-order transitions in lattice gauge theories.By analyzing the density of states we aim to avoid the super-critical slowing down that plague Markov-chain importancesampling techniques at such first-order transitions.Motivated by the Stealth Dark Matter model and the stochastic gravitational-wave background that would be produced by a firstorder dark-sector confinement transition in the early Universe, our ongoing work focuses on SU(4) Yang-Mills theory.
Our current results from 12 3 × 6 lattices are insufficient to resolve the SU(4) confinement transition.While a smaller energy interval size δ could help, given the small latent heat of this transition, this will also make it more challenging to control statistical uncertainties as shown by Fig. 1.We are currently improving our analyses by analyzing larger N 3 s × N t lattice volumes, both increasing N s with fixed N t to study the thermodynamic limit and increasing N t with fixed aspect ratio N s /N t to extrapolate to the continuum limit.Once we have resolved the first-order confinement transition, we will be able to predict observables such as the latent heat and surface tension that affect the resulting spectrum of gravitational waves.
From our earlier results for compact U(1) lattice gauge theory, we can appreciate that, in contrast to standard importance sampling approaches, the LLR method may be easiest to apply to study very strong first-order phase transitions characterized by a large latent heat.Because the first-order SU(N) bulk transition for N > 4 is notably stronger than the finitetemperature transition, we are separately using the LLR algorithm to study such large-N bulk transitions and explore the expected scaling of the latent heat ∝ N 2 [14].We will soon present preliminary results for the SU(5) and SU(6) bulk phase transitions.
Based on these studies of both the large-N bulk transition and SU(4) confinement transition, we will be in a good position to consider the confinement transition for SU(N) Yang-Mills theories with larger N. We will also be able to estimate the computing resources that would be required to analyze the weaker first-order confinement transition of SU(3) Yang-Mills theory corresponding to quenched QCD.Finally, we are exploring broader applications of the LLR algorithm, including to phase transitions in bosonic matrix models of interest in the context of holographic gauge/gravity duality.

Figure 1 .
Figure1.Results for a obtained for compact U(1) pure-gauge theory with lattice volume V = 84  and energy intervals of size δ = 0.01/V (left) or δ = 0.001/V (right).We can see that the statistical uncertainties increase in the limit δ → 0 considered in Eq. 8.

,Figure 3 .
Figure 3. SU(4) results for a from 6 4 lattices with an energy interval size of δ = 0.01/V.The statistical uncertaintiess are obtained by performing N j = 5 independent runs per interval.

Figure 4 .
Figure 4. SU(4) results for a from V = 12 3 × 6 lattices with an energy interval size of δ = 0.001/V.The error bars are estimated by performing N j = 5 independent runs per interval.There is no sign of a first-order phase transition, which would correspond to a non-monotonic a(E).

,Figure 5 .
Figure5.The left plot compares numerical data for a (black) with the synthetically generated values (blue line) that would be needed to reconstruct a similar first-order confinement transition to the one in Ref.[19].Both the synthetic and the measured a use lattice volume V = 12 3 × 6.The right plot shows the probability density P β (E) = ρ(E)e βE reconstructed from the synthetic a(E) with β ≈ 10.8.A small double-peak structure is visible, confirming a first-order transition.