The Benefits of Diligence: How Precise are Predicted Gravitational Wave Spectra in Models with Phase Transitions?

Models of particle physics that feature phase transitions typically provide predictions for stochastic gravitational wave signals at future detectors and such predictions are used to delineate portions of the model parameter space that can be constrained. The question is: how precise are such predictions? Uncertainties enter in the calculation of the macroscopic thermal parameters and the dynamics of the phase transition itself. We calculate such uncertainties with increasing levels of sophistication in treating the phase transition dynamics. Currently, the highest level of diligence corresponds to careful treatments of the source lifetime; mean bubble separation; going beyond the bag model approximation in solving the hydrodynamics equations and explicitly calculating the fraction of energy in the fluid from these equations rather than using a fit; and including fits for the energy lost to vorticity modes and reheating effects. The lowest level of diligence incorporates none of these effects. We compute the percolation and nucleation temperatures, the mean bubble separation, the fluid velocity, and ultimately the gravitational wave spectrum corresponding to the level of highest diligence for three explicit examples: SMEFT, a dark sector Higgs model, and the real singlet-extended Standard Model (xSM). In each model, we contrast different levels of diligence in the calculation and find that the difference in the final predicted signal can be several orders of magnitude. Our results indicate that calculating the gravitational wave spectrum for particle physics models and deducing precise constraints on the parameter space of such models continues to remain very much a work in progress and warrants care.


Introduction
Even the very early Universe is transparent to gravitational waves, making searches for the gravitational wave background of the Universe a unique probe of the cosmos before big bang nucleosynthesis. Ubiquitous in the literature is the generation of a gravitational wave background from an inhomogeneous transition of the ground state (for a review see [1][2][3]). In the standard model of particle physics, there is no mechanism for such a gravitational wave background to be produced. Specifically, both the QCD and electroweak transitions are predicted to be smooth [4][5][6][7][8][9]. This implies that any gravitational wave background resulting from a strong first order phase transition is proof that the standard model is incomplete. The electroweak transition can be made strongly first order through the introduction of new states at around the electroweak scale . The QCD transition can be catalyzed by changing the number of light fermions [44] or having a very large lepton asymmetry [45][46][47]. Additionally there are strong motivations to believe that the standard model is incomplete and additions to the standard model can also leave cosmic fingerprints. For instance, baryonic matter can only explain a fraction of the matter observed and the missing dark matter can be a part of a hidden sector that undergoes a phase transition [48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65]. Second, the near unification of gauge coupling constants along with conspiracy of gauge anomaly cancellation motivates grand unification which can sequentially break into the standard model gauge group and leave a gravitational wave background [66][67][68][69][70][71][72]. Finally, the generation of neutrino masses can arise through a B − L breaking transition [68,[73][74][75][76][77]. In each case, an observed signal not only sheds light on our cosmic history, but on a range of energy scales spanning from sub-GeV to the PeV scale [78] (even higher scales have been proposed, though technology needs to improve to make the sensitivity cosmologically relevant [79] with the possible exception of NEMO [80]).
Any strong first order transition produces three contributions to a stochastic background (see, e.g., [1][2][3]). For a transition in the thermal plasma, only a negligibly small fraction of the energy released will remain in the scalar field when the bubble wall reaches a constant velocity [81], or when it "runs away" [82]. So the bubble collision contribution is negligible when considering transitions in a thermal plasma, though it can become dominant for transitions in vacuum, e.g., those in a dark sector without a thermal bath. In addition, the contribution from magnetohydrodynamic turbulence is subdominant compared with that from sound waves [81,83,84], with an efficiency factor that is roughly (5 ∼ 10)% of that for sound waves [85], though its spectrum is highly uncertain as of now [86][87][88][89][90][91][92]. Since most transitions in the early universe are highly likely to be proceeding in the thermal plasma, we will focus here on the gravitational wave production from sound waves. This acoustic contribution has been studied both in simulations and a combination of analytic and numerical techniques and there has been much recent progress.
Given the enormous opportunity to shed light on both cosmology and particle physics, it is worth examining in detail the theoretical underpinnings of any given model in order to enumerate both theoretical uncertainties in basic methods and the degree of benefit in more accurate calculations or, equivalently the cost of various approximations. Approximations can arise in two steps in predicting an observable from a given model as shown in Fig. 1. First the calculation of macroscopic thermal parameters, including the latent heat and the time scale of the transition, are often calculated using perturbative techniques which can introduce large errors [93] in particular when ℒ Macroscopic Thermal parameters Thermal evolution of the effective potential Observables Figure 1. The uncertainty in linking a particular model with a set of observables is conceptually presented above. The break down of perturbation theory at finite temperature is the dominant error in the prediction of the evolution of the effective potential and ultimately non-perturbative methods might be required to predict macroscopic thermal parameters. The macroscopic thermal parameters of interest are often taken to be the latent heat, the time scale of the transition (usually approximated), the bubble wall velocity and the temperature of percolation, but if one desires to have an accurate prediction one needs the fluid velocity, the wall velocity, the mean bubble separation, the percolation temperature and the lifetime of the acoustic source (see also Fig. 1 of [96]).
long wavelength modes are not resummed carefully enough [94][95][96][97]. 1 The second step, which we focus on in this paper, converts macroscopic thermal parameters into a prediction for the spectrum -in particular the peak frequency and amplitude. Ultimately, both steps will likely require simulations to truly perform precision cosmology on a future hypothetical observation. 2 However, this is impractical for the analysis of large numbers of parameter sets for large numbers of models. We therefore examine several layers of improvement in the prediction of the peak amplitude that have recently arisen in several models involving physics beyond the standard model • The finite lifetime of the source first estimated in Ref. [102] and derived in the sound shell model in an expanding background in Ref. [103].
• Going beyond the bag model approximation in solving the hydrodynamic equations [104,105].
• Calculating the mean bubble separation from the evolution of the bubble number density.
• Calculating the fraction of energy in the fluid from solving the hydrodynamic equations rather than using a fit [106].
• Including fits for the energy lost to vorticity modes [106].
In this paper we will enumerate the error in a number of models in order to get a broad understanding of the numerical importance of diligence. This avoids model-specific effects where accidental cancellations between different improvements could in principle occur. The models we consider include a toy model introduced for pedagogical purposes, the Standard Model Effective Field Theory (SMEFT), a dark sector Higgs and a real scalar singlet extension (xSM) of the Standard Model. For the benefit of the reader, a demonstration of the importance of diligence is 1 Other important problems in common calculations are gauge dependence [98] and the inhomogeneous background [99,100]. 2 Infrared divergences of the dynamical mode for instance remain even after NLO resummation. As a result perturbation theory even at two loops disagrees substantially with montecarlo simulations very close to the critical temperature [101]. 30  provided at the outset in Fig. 2. Here, the relative error in the predicted peak amplitude is shown for SMEFT, the dark sector Higgs model (which we label throughout as "Dark RG") and xSM. Our paper will be devoted to fully explaining Fig. 2; for now, we provide a feel for the comparative importance of these errors. For the Dark RG, for example, the relative error is far more manageable than what it is for SMEFT. However, even for that model, the relative error is larger 3 than the error from gauge dependence that is introduced in SMEFT in some commonly used methods. Thus, even this case, which may present an unrealistically optimistic picture, still motivates diligence in the calculation. The structure of this paper is as follows. In Section 2 we outline three methods of various levels of diligence that we find used in the literature, including a level of diligence motivated by its use in the recent review [3]. In Section 3 we define the models we will use to demonstrate increasing levels of diligence, and in Section 4 we will present our results. We will end with our Conclusion.

Phase Transition Dynamics
Gravitational waves produced from first order phase transitions is a finite temperature tunneling process, from some false vacuum to the true vacuum. When calculating this transition with perturbation theory, one needs to track the minima of the effective potential from the temperature at which the energy in each vacuum is degenerate -that is the critical temperature. Below the critical temperature, bubbles of the true phase begin to form at some critical radius where the pressure is strong enough to cause expansion. The probability of such bubbles forming increases as the Universe cools, until the nucleation temperature at which there is an average of one bubble per Hubble volume. Slightly below this temperature is the percolation temperature at which bubble collisions are occurring and the final temperature when the phase transition ends. There are simple analytical expressions for these temperature scales which are the result of approximations used in the equations. However, the gravitational wave spectrum is sensitive to the level of diligence that goes into the computations and reducing the error is paramount to probing phase transitions at future gravitational wave detectors. We will now proceed to analyze the different level of diligence used in the literature.

Lowest diligence
Here we will describe the level of lowest diligence in computing the gravitational wave spectrum. At this stage, we will only introduce the various parameter definitions and wait until the highest diligence section for a more in depth look at the numerical procedure. This level will involve computing all relevant parameters at the nucleation temperature.
The tunneling rate per unit time per unit volume will have the general form wherep 0 is a dimensionless number that we will assume is O(1) and S 3 can be found by solving the bounce solutions that minimize the action given by The nucleation temperature is defined as the temperature at which the probability of a single bubble being nucleated within a Hubble volume is O(1): where M pl is the Planck mass, ξ ∼ 3 × 10 −2 , and V H (t) is the horizon volume. This equation will lead to the simple definition of the nucleation temperature [3,[107][108][109] It is important to note that the above calculation assumes that the phase transition occurs in a radiation dominated era which is not guaranteed. The strength of the gravitational wave spectrum will depend on hydrodynamic parameters such as the amount of vacuum energy released during the phase transitions, the inverse time duration of the phase transition, and the fraction of latent heat that goes into the bulk motion of the plasma (referred to as the kinetic efficiency coefficient). We discuss each of these quantities in turn.
The strength of the phase transition is characterized as where V eff is the finite temperature effective potential and the symbol ∆ signifies the difference in the symmetric phase (false vacuum) and the broken phase (true vacuum). The energy density of radiation is given by ρ rad = π 2 /30 g * T 4 n where g * is the number of effective degrees of freedom at T n .
The inverse time duration of the phase transition evaluated at the nucleation temperature can be approximated as where H 2 n = 8πGρ rad (T n )/3 is the Hubble parameter at the nucleation temperature. A smaller β/H and larger α will result in stronger gravitational waves.
The gravitational wave spectrum observed today has a simple broken power law fit [81] in terms of the aforementioned parameters given by where g n is the number of degrees of freedom at the nucleation temperature and κ is the efficiency coefficient that represents the fraction of the bulk kinetic energy in the plasma relative to the available vacuum energy. The numerical fits for the kinetic efficiency coefficient, κ, were derived in [110] for the different velocity profile types which we give in Appendix A. The spectral shape, S SW , and the peak frequency, f SW are given by Hz. (2.9) The gravitational wave spectrum may be rewritten in terms of the R.M.S velocity, U 2 f = 3 4 κα, with the replacement where Γ ∼ 3/4 is the adiabatic index which is defined as the ratio of the enthalpy and energy density in the symmetric phase. The term in the denominator on the left hand side, (1 + α), is the result of the energy density in the symmetric phase.

Moderate diligence
The level of modest diligence is the approach most frequently used in the recent literature (including the recent LISA review [3]). It closely resembles the lowest diligence with the exception that the thermal parameters are defined at the percolation temperature rather than the nucleation temperature and the finite lifetime of the source is taken into account with an ansatz correction to the peak amplitude. The percolation temperature is here approximated by solving the equation where log(A/T 4 ) ∼ 14 for an electroweak phase transition. Note that the derivative of the left hand side in Eq. 2.11 appears on the right hand side, as can be seen from Eq. 2.6. The percolation temperature is always below the nucleation temperature and hence closer to the final temperature when the phase transition ends. This makes using the percolation temperature a better approximation to estimate the thermal parameters. However, if the percolation temperature is significantly far away from the nucleation temperature, one should check if the phase transition can even reach completion at all for cases of strong supercooling, since the universe may become vacuum dominated. The strength of the phase transition and the inverse time duration of the phase transition take on the same form as in Eq 2.5-2.6 but with the replacement T n → T p such that and where H p is now the value of the Hubble parameter at the percolation temperature. The gravitational wave spectrum in Eq. 2.7 assumed that the lifetime of the source is approximately one Hubble time, Hτ sw = 1 [85] . It was later pointed out in [74], that a better approximation to the lifetime of the source is where R * is the mean bubble separation and U f is the root mean squared velocity defined at α(T p ). The mean bubble separation is related to the inverse time duration using R * = (8π) 1/3 v w /β. We then take into account the finite lifetime of the source in the gravitational wave spectrum through and calculate all temperature dependent quantities at the percolation temperature T p defined in Eq. 2.11.

High diligence
The highest diligence with which one can calculate the gravitational wave spectrum involves a number of improvements to the predictions of the peak frequency and amplitude: 1 Improving on the bag model approximation for the fluid velocity and fraction of energy that is in gravitational waves; 2 Calculating the fluid velocity and efficiency from solving the hydrodynamic equations rather than using fits (related to the first); 3 Calculating the mean bubble separation from the number density of the bubbles; 4 Taking into account the finite lifetime of the soundwave source, derived in an expanding universe [103]; 5 Calculating the suppression due to reheated droplets creating friction that slows collisions.
Note that Ref. [111] used the bag model in their simulations, so we assume that the suppression factor arising from kinetic energy lost in the fluid is independent of the change in the amplitude from improving on the bag model. Also in the last case, and only the last case, we use fits to estimate this degree of suppression as it relies on a full numerical simulation -methods to approximate this effect we leave to future work. In this section we outline in detail each of the other improvements.

Calculation of the percolation temperature
The false vacuum fraction at t > t c in an expanding universe is defined as where I(t) represents the volume of nucleated bubbles per comoving volume, double counting the overlapped space between bubbles and virtual bubbles within others [112]. The comoving radius of a bubble nucleated at t and measured at t is for FLRW space where η is the conformal time. This form of the comoving radius of the bubble assumes that the initial size of the bubble is negligible and the bubble wall velocity is constant. The above equations can be recast in terms of integrals over temperature in order to render them more suitable for numerical computations. The scale factor of an expanding universe will have the form a = c a t n where c a is a constant and n is determined by the form of energy that is dominating the expansion, i.e. n = 1/2 and n = 2/3 for radiation and matter respectively. The expansion of the universe may be treated as an adiabatic process so that the entropy s R in the radiation sector is conserved per comoving volume: (2.18) The radiation sector will have s R ∝ T 3 which will typically give T ∝ 1/a ∝ t −n . This is valid for both a radiation dominated universe and a matter dominated universe. However if non-relativistic matter can decay into radiation, the dependence is T ∝ a −3/8 [113]. In general, we may assume where γ will depend on the amount of entropy injection. The Hubble parameter is related to the rate of change of the scale factor, H =ȧ/a, and we then have which allows for rates of change with respect to time to be exchanged with that of temperature. We can also relate the scale factor at the time of the critical temperature, a c , to the scale factor at a later time by For an expanding universe in a radiation dominated era, the total energy density is which consists of the relativistic matter density ρ rad , and the vacuum energy released during the phase transition. Accordingly the Hubble parameter in an expanding universe in this notation yields with y = a/a c . The above definitions can be used to convert any integral over time to that of temperature so that Eq. (2.16-2.17) may be written as With the false vacuum fraction now defined as a function of temperature, the percolation temperature occurs when the false vacuum fraction is 70% of the total volume, i.e., when g(T p ) = 0.7. (2.26) Similarly, the temperature when the phase transition ends occurs at the time when the volume of nucleated bubbles equals the comoving volume, i.e. I(t f ) = 1. This translates into In most cases, the percolation temperature calculated from Eq. 2.26 should roughly coincide with the final temperature calculated from Eq. 2.27 and depend on v w and ∆V .

Calculation of the nucleation temperature
The evolution of the mean bubble density per proper volume is determined by which begins at n b (t c ) = 0 and includes all the bubbles ever formed. We then transform the above result into an integral over temperature: The nucleation temperature can be numerically calculated as the temperature when there is one bubble per Hubble volume: The mean bubble density depends on the false vacuum fraction and therefore should change with v w , κ m , and ∆V as well.

Calculation of the mean bubble separation
The mean bubble separation R * is related to the mean bubble number density and is given by The gravitational wave spectrum is usually written in terms of the inverse time duration of the phase transition where R * and β are interchangeable. This form of the mean bubble separation can similarly be found in an expanding universe by Taylor expanding Eq. 2.1 about η 0 with respect to conformal time for exponential nucleation where β c is the comoving inverse time duration of the phase transition: The above Taylor expansion is valid in the regime where the false vacuum fraction decreases faster than the rate p increases which occurs when I(t) becomes order 1 and most of the bubbles are nucleated. The exponent of the false vacuum can be transformed into an integration over conformal time, dt = adη: which can then be used to calculate the false vacuum fraction for η > η c . If we define the comoving bubble number density as n b,c = n b a 3 , we can rewrite Eq. 2.28 in terms of conformal time: We now integrate over η to have where we used the relation that I(η f ) = 1 for η 0 ∼ η f . The relation between the comoving mean bubble density and the mean bubble separation is now where β c and R * ,c are both functions of the bubble wall velocity. The mean bubble separation at conformal time η is now related to the inverse time duration when the phase transition ends: where R * ,c = R * /a(η) and β c = a(η f )β. The above result can be transformed into a function of temperature using a relation similar to Eq. 2.21 to give The procedure to calculate the inverse time duration of the phase transition at T f for a fixed v w is

Going beyond the Bag model
The free energy density f of a model with g * degrees of freedom consists of a zero temperature scalar potential and a thermal potential that involves the summation over all relativistic species that interact with the scalar φ. In the Standard Model, this involves the standard electroweak Higgs field and degrees of freedom g SM = 106.75. The free energy density gives where the summations run over all relativistic bosons B and fermions F at temperature T . In the high temperature expansion, the free energy density can be written as where V T (φ) is the thermal effective Higgs potential and we explicitly separate out the scalar independent terms that go as T 4 . The hydrodynamics of the plasma can be described as a perfect fluid in terms of the energy density e, pressure p, and enthalpy ω. These thermodynamic quantities can be explicitly calculated from the free energy density through the relation p = −f (φ, T ) for the pressure. The energy density and enthalpy are then related to the pressure through These quantities together, along with the velocity of the fluid v, will give the energy momentum tensor of the plasma T µν = u µ u ν ω + g µν p, (2.43) where g µν is the inverse Minkowski spacetime metric and u µ (v) is the four velocity of the fluid. During the period of gravitational wave production, bubbles of the new phase will collide and generate sound waves in the perturbed plasma. The velocity profile around a single bubble is determined by the hydrodynamic quantities and the junction conditions across the bubble wall. We will denote the broken phase inside the bubbles by subscripts "b, −" and the symmetric phase outside the bubbles as "s, +". The signs "−, +" are used to describe quantities behind and in front of the bubble wall. The continuity equations ∂ µ T µν = 0 can be integrated in the wall frame across the bubble wall to give the junction conditions where v ± and T ± corresponds to the velocity and temperature of the fluid behind and in front of the bubble wall. We may assume that T + ≈ T − and expand the thermodynamic quantities about the symmetric phase so that the ratio of velocities about the junction give where c s,b is the speed of sound in the broken phase and αθ + is the phase transition strength [105]. The speed of sound can be calculated in both the symmetric and broken phase and is defined in terms of the pressure as where (...) denotes derivatives with respect to temperature. In terms of the free energy density, the pressure in the symmetric phase is φ independent with p s = −f (0, T ) and the pressure in the broken phase is p b = −f (φ min , T ) which is evaluated at true vacuum φ min . It is important to keep every term in the free energy density when calculating the speed of sound in order to properly account for the full temperature dependence of the model. This includes keeping all light degrees of freedom that do not acquire field dependent masses that affect the Higgs effective potential.
The new phase transition strength parameter αθ + is dependent on the speed of sound in the broken phase and has a similar form to the bag parameter built up from the trace anomaly, [105]. Going back to the free energy density, we can define the phase transition strength as a function of temperature where V T is the Higgs effective potential defined in Eq. 2.40-2.41 and ∆ denotes the difference between the symmetric and broken phase. The phase transition strength has the same form as the bag model α θ defined in Eq. 2.5 when c 2 s = 1/3 and ω s = 4/3ρ rad . The junction conditions to obtain v ± will serve as boundary conditions for solving the velocity profile obtained from projecting the continuity equations into the parallel and perpendicular motions of the fluid flow. The hydrodynamic equations become where ξ = R/t is a self similar coordinate defined as the ratio between the distance to the bubble center and the time since nucleation. These hydrodynamic equations may be combined to give the enthalpy profile and the velocity profile: where ∂ ξ = ∂ ξ e + ∂ ξ p and c 2 s = dp/de are used to connect the equations. The Lorentz boost transformation used in the equations is defined in terms of the self similar coordinate µ(ξ, v) = (ξ − v) / (1 − ξv). In detonations, the fluid velocity ahead of the bubble wall is always zero so that the hydrodynamic profiles are independent of the speed of sound in the symmetric phase. Deflagrations have a non-zero bubble wall velocity ahead of the bubble wall and the equations will then depend on the speed of sound in the symmetric phase. Both profile types will always depend on the speed of sound in the broken phase through the junction conditions. In detonation, it is sufficient to use αθ + = αθ(T ) with T usually taken as the nucleation temperature or the percolation temperature. Deflagrations and hybrid modes take αθ as input and require a shooting method by varying αθ until αθ is reached far away from the bubble.
The quantity of interest for the peak gravitational wave energy density is the kinetic energy fraction K which can be solved by the hydrodynamic equations: where ρ fl = 3/v 3 w dξξ 2 v 2 γ 2 ω is the fluid's kinetic energy. We use the publicly available code in [105] to numerically compute the kinetic energy efficiency κθ for a given set of c 2 s and αθ when comparing to calculations in the bag model. The kinetic energy fraction is related to the efficiency parameter through where κ = 4ρ fl 3αθω s (2.56) [105]. The quantities c s,s , c s,b , e s ,ω s , αθ, and κ in determining the kinetic energy fraction K are all calculated at the final temperature T f when the phase transition ends. The enthalpy-weight rootmean-square fluid velocity around a single bubble may be found from the kinetic energy fraction, which is evaluated at T f .  . Suppression factor with respect to the strength of the phase transition due to vorticity and reheating effects in the plasma.Ω is the reduced peak gravitational wave energy density and Ω exp is the expected peak gravitational wave energy density. The data is taken from [106].

Gravitational Wave Spectrum
The first numerical simulations of strong first order phase transitions for α 0.1 were undertaken in [106] which showed that previous results overestimated the gravitational wave energy density. The rotational component of the fluid velocity, particularly for deflagrations, increases in relative size as the transition strength grows. This reduces the amount of available kinetic energy that can be transferred to the fluid. The rotational component for detonations, however, remains small and constant. The probable explanation for this effect is the formation of reheated droplets of the metastable phase that are produced during the collisions of the bubble walls. These droplets can then slow down the bubble walls and reheat the surrounding regions. The simulations considered a simple bag equation of state where the results only depend on input parameters of v w and α. For a given v w , we use an interpolation of their results to estimate the corresponding suppression factor to the Ω GW h 2 as a function of α. We show the suppression factor in Fig. 3 for two representative bubble wall velocities. We utilize extrapolation for when α is out of range. Although the results were performed using a bag equation of state, we numerically compute α in the beyond the bag model and assume that the suppression from vorticity and reheating derived with a bag model applies without modification. We will test this assumption in future work. Furthermore, the simulations suggest that the RMS fluid velocityŪ f reaches a maximum that is under-approximated by calculating the expectedŪ f around a single bubble. We use the results ofŪ f,max /Ū f,exp to estimate the maximum fluid velocity in the highest diligence after calculatingŪ f in the beyond the bag model. A careful calculation of the gravitational wave production in an expanding universe will result in a suppression factor of the form for a radiation dominated era where y = a/a s is a dimensionless scale factor ratio normalized by when the source of production becomes active. For a radiation dominated era, the time elapsed since some reference time, t s , is Due to the presence of shocks and turbulence in the plasma, the time elapsed is unlikely to last an arbitrarily long time. The effective lifetime of the source is then given by the timescale of turbulence, τ sw = R * /Ū f . This was the basis of the suppression factor used in Eq. 2.14-2.15. We use the estimated maximum of the fluid velocity in the beyond the bag model when calculating the lifetime of the source. The gravitational wave energy density will then be suppressed by which approaches the asymptotic value Υ = 1, the lowest diligence, when τ sw H s → ∞. When τ sw H s 1, the suppression factor is approximately Υ = τ sw H s , the modest diligence. The peak gravitational wave energy density after taking into consideration the suppressions arising from vorticity and reheating effects in the plasma as well as the lifetime of source becomes where K is calculated in the beyond the bag, β/H * is calculated from the mean bubble separation, and the last factor arises from the vorticity and reheating effects in the plasma.

Test models
In this section we examine the numerical difference in predictions arising from different levels of diligence in several models

SMEFT
SMEFT is a model independent method of examining many extensions of the Standard Model by augmenting it with a tower of high dimensional operators, each suppressed by higher and higher powers of the cutoff scale corresponding to the scale of new physics. Unfortunately, the Standard Model requires such a large change to its potential that the scale of new physics needs to be quite low to augment a strong first order phase transition [37]. In such a case the SMEFT only provides a qualitative description of the UV complete scalar sector, and then only in special circumstances [37]. In the SMEFT the tree level potential is augmented by a single higher dimensional operator where M characterizes the cut off scale, φ is the SM Higgs doublet and ∆V h is chosen such that the zero-temperature minimum is shifted to the origin. We will consider the full free energy density at one-loop given by where V CW is the Coleman-Weinberg contribution and V T is the finite-temperature correction. These are given by and 5) where N α counts the number of degrees of freedom of each particle and C α is a constant that is 5/6 for gauge bosons and 3/2 for all others. We note that the daisy terms in Eq. 3.5 are the result of a high temperature expansion which may cause an IR-divergence in the speed of sound for low temperatures [105]. We explicitly check this by including a Boltzmann suppression term when M α 2.2T . The sums run over the top quark, W and Z bosons, and the Higgs boson h. The total degrees of freedom in SMEFT is the Standard Model value g SM = 106.75. The calculation of the speed of sound requires including all the relativistic particle species in the free energy density. We will account for the remaining light particles that were neglected in V T by including the term: to the free energy density where g * = 345/4. However, in the bag model, the speed of sound is taken to be c 2 s = 1/3 and the light species can be ignored as they do not affect the phase transition dynamics. The last term in V T corresponds to the resummation of the daisy terms of the scalar bosons. To calculate the effective potential and the counter-terms at zero-temperature, we fix the zero-temperature MS-parameters by matching the physical observables at the Z boson pole mass m Z using the full self energies. To go beyond the bag model, we need the absolute pressure in each phase, and not just the relative pressure. We therefore add an overall constant in the potential such that the pressure in the broken phase at zero temperature vanishes at one loop. The scale of the Coleman-Weinberg potential is taken to be at µ ∼ T for the dynamics of the phase transition and we run the parameters to this scale.

Dark Renormalizable Models
Here we will consider a dark Higgs model [48,50,62,63] of the type SU (N )/SU (N − 1) with renormalizable operators following the conventions in [50]. The overall scale Λ and the zero temperature vacuum expectation value v are the only inputs of the model. We can then define the zero temperature parameters such as the tachyonic mass and self coupling as where factors of v/Λ will control the thermal parameters. The tree level potential is then defined as where h D is the dark Higgs of the doublet H and ∆V shifts that potential at the minimum to zero. The degrees of freedom of the full dark sector in consideration are where n G is the number of Goldstone bosons, n GB are the gauge bosons, N f is the number of fermions,and N is the rank of the group. The free energy density in consideration is where V CW and V T are defined in Eq. (3.3-3.5). Here the summations now only run over the dark sector particles. This requires us to add on the additional relativistic particles not included in the sum which now include the full degrees of freedom of the Standard Model: where g * = 106.75. We add the term ∆V so that the minimum of the tree-level potential is shifted to zero. We choose the scale of the one-loop potential to be at Λ. The inclusion of the CW-term will shift the zero-temperature vacuum expectation value away from v. We add the counter-term potential where δµ 2 , δλ, and δ∆V h are chosen such that 14) to maintain the tree-level structure of the potential. We work in the Landau gauge where the Goldstone bosons have zero mass at the h D = v which causes an IR-divergence in the one-loop potential. This causes an issue in the evaluation of the counter-terms. One prescription is to remove the Goldstone bosons from the sum in the CW potential. For the purpose of this work, we will follow the procedure in [115] to evaluate the counterterms. This shifts the Goldstone mass by its mass at the one-loop level, i.e is the field dependent mass of the Goldstone bosons.

xSM
The singlet extended SM, known as "xSM", consists of the standard SM Higgs doublet H T = G + , v EW + h + iG 0 / √ 2 and a real gauge singlet S = v s + s where the electroweak vacuum is (v EW , v s ) [28,38,[116][117][118][119][120][121][122]. The tree level potential in this setup is defined as where ∆V shifts the minimum of the potential to zero. The mass parameters µ 2 and b 2 are related to the other model parameters through the minimization conditions around the electroweak vacuum, The parameters λ, a 1 , and a 2 are related to the physical parameters θ, m h 1 , and m h 2 through the mass matrix diagonalization as whereb 3 ≡ b 3 /v s and θ is the mixing angle of the physical fields h 1 and h 2 defined as with s θ ≡ sin(θ) and c θ ≡ cos(θ). Here we associate h 1 as the SM Higgs and h 2 is some heavier scalar. The free energy density we consider in the xSM presented here contains only the high temperature expansion approximation of the full finite temperature one loop effective potential since the phase transition is primarily driven by the cubic terms. The free energy is then where we take g * = 106.75. The field dependent thermal mass Π h (T ) and Π s (T ) are where the physical masses of the W , Z, and t-quark are used to define the gauge and Yukawa couplings to h.

Results
The resulting gravitational wave spectrum is dependent on the level of precision of the thermal parameters. Until recently, the bag model was assumed to compute the phase transition strength and the kinetic energy of the fluid. Going beyond the bag model will require the calculation of the speed of sound in the plasma for both the symmetric and broken phase which should result in a more accurate treatment of the thermal parameters. These quantities are temperature dependent and will change depending on the temperature at which they are computed. Furthermore, the temperature scales of the phase transition such as the nucleation and percolation temperature are also sensitive to level of diligence in the calculations. We use the publicly available codes CosmoTransitions [108] and BubbleProfiler [123] to compute the actions in order to find the relevant transition temperatures. The lowest diligence level will compute the thermal parameters at the estimated nucleation temperature defined in Eq. 2.4. The strength of the phase transition is calculated using Eq. 2.5, the inverse time duration of the phase transition is calculated using Eq.2.6, and the peak gravitational wave energy density is calculated using Eq. 2.7.
The modest diligence level will compute the thermal parameters at the estimated percolation temperature in Eq. 2.11. We will use Eq. 2.12, Eq. 2.13, and Eq. 2.15 to estimate the strength of the phase transition, inverse time duration, and peak gravitational wave spectrum respectively. The lifetime of the source is estimated using Eq. 2.14. There is an ambiguity as to which temperature to  employ in the calculation of the thermal parameters: the nucleation temperature or the percolation temperature. The true percolation temperature T p , Eq. 2.26, should lie close to the temperature at which the phase transition ends T f . The inverse time duration should be computed from the mean bubble separation, Eq. 2.31 and Eq. 2.31, which is evaluated at T f . The highest diligence will evaluate all thermal parameters at T f to ensure that all quantities are evaluated at the same temperature as the inverse time duration. From here on out, we will associate T p with Eq. 2.11 when referring to modest diligence. The highest diligence will also utilize the beyond the bag model to calculate the strength of the phase transition Eq. 2.49 and the kinetic energy fraction Eq. 2.55 which requires the numerically calculated speed of sound in Eq. 2.47. The lowest and modest diligence calculations will assume c 2 s = 1/3 as is done in the bag model. The peak gravitational wave spectrum is found from Eq. 2.61 which accounts from the lifetime of the source, Eq. 2.60, as well as vorticity and reheating effects in the plasma. We note that the suppression factor due to the finite lifetime used in Eq. 2.14 is a valid approximation of Eq. 2.60 only when τ SW H 1. In the following sections, we will compare the different levels of diligence in SMEFT, the dark  renormalizable model, and xSM. The error in the peak gravitational wave energy density, where j = (low, mod) refers to lowest, modest, and highest diligence respectively, will be used to compare the different levels of diligence.

SMEFT
The SMEFT model we consider has the scale of the zero-temperature one loop potential set to µ = T as well as temperature dependence in the running of the couplings. This will contribute additional temperature dependence to the speed of sound calculations in the broken phase. We note that Ref. [124] also considered the beyond the bag model in SMEFT using a high temperature expansion of the effective potential. The speed of sound will never reach the bag model with c 2 s = 1/3 as seen in the top figure left of Fig. 4. The green curves show the different levels of diligence for the speed of sound in the broken phase and the dashed gray curve represents c 2 s = 1/3. The magenta curve is the speed of sound calculated in symmetric phase which is approximately the same in each level and does not deviate far from the bag model. We do not consider any additional relativistic degrees of freedom and thus expect little deviations between the speed of sound in the symmetric phase. As the scale M grows large, the speed of sound in the broken phase approaches a constant value of c 2 s ∼ 0.32. There is noticeable disagreement between the different levels below M = 600 where there is mild supercooling. For a given M , the speed of sound is only a function of temperature. The differences in c 2 s in the broken phase is the result of these different temperatures at which the speed of sound is set to when calculating the strength of the phase transition αθ(c 2 s ). The large difference in T p and T f is due to the approximations of T p in Eq. 2.11 which is less accurate when S 3 /T acquires a minimum for smaller M .  On the top right panel of Fig. 4 we show the strength of the phase transition computed at the different levels of diligence. Both the lowest and modest diligence curves have c 2 s = 1/3 whereas the highest diligence curve corresponds to the beyond the bag calculation with c 2 s shown in the top left panel. Although each level is computed at different temperatures, the lowest diligence is a better approximation of the strength of the phase transition compared to level 2 which over approximates α. This is a result of T p computed in the modest diligence placing far below T f which results in a higher estimated α. The difference between the different levels on the strength of the phase transition is negligible for large M as a result of the asymptotic behavior observed in c 2 s and the better approximation of T p when there is no barrier present. There is also a dependence on the bubble wall velocity in both c 2 s and α for the modest and highest diligence curves in computing T p and T f but we only show detonation with v w = 0.92 because the difference is minor. The kinetic energy fraction is shown in the bottom panel of Fig. 4 and should depend on the speed of sound and phase transition strength. Similar to what was seen in α, the lowest and highest diligence curves are closer together for large M while the modest diligence curve is the least accurate. For each of the levels, the largest error in both α and K occurs for smaller M where the speed of sound is significantly lower than c 2 s = 1/3 and T p is far from T f . In the left panel of Fig. 5, we show the inverse time duration β/H of the phase transition for detonation. The largest difference between modest diligence and highest diligence occurs for small M and is due to the following reason: the minimum formed in the action where T p calculated in Eq. 2.11 along with β/H in Eq. 2.13 are inaccurate when there is a minimum present. The lowest diligence is a better approximation than the modest diligence in this regime. The modest and the highest diligence become indistinguishable for large M when there is no minimum in the action. For small M , the lowest diligence curve appears to be a good approximation for modest diligence. Although β/H estimated from the action is not accurate when there is a minimum, the error using  T n appears to do better than using the approximation of T p . Contrary to the α and K, β in the lowest diligence never approaches the highest diligence for large M where the error appears to get worse. This is due to the inaccuracy in using the approximate T n . In this regime, T p is a better approximation of the inverse time duration as there is no minimum present in the action. The right panel of Fig. 5 shows the suppression factor due to the lifetime of the source in the modest and highest diligence. The error between the two levels gets worse for small M which is the result of the error in T p . The error approaches a constant as M gets large.
In Fig. 6 we show the relative error in the peak gravitational wave spectrum for both deflagration and detonation. The error with respect to the highest diligence is estimated using Eq. 4.1. For both deflagration and detonation, the modest diligence spikes to an error of ∆Ω/Ω ∼ 10 2 for small M. This is correlated with the large error observed in α, K, β, and Υ seen in Fig. 4 and Fig. 5. The modest diligence error for both profile types slowly approach a constant value as M grows large which is the result of the minimal error in α, β/H, K. The error in Υ appears to become a constant for large M . The suppression factors due to the lifetime of the source grow to zero as M grows large which results in the increasing behavior of the peak error in the lowest diligence which does not include any suppression factors. Overall we notice an error in the peak gravitational wave energy density of 10 1 − 10 2 for lowest diligence and 10 0 − 10 2 for modest diligence.

Dark Renormalizable Models
The dark renormalizable model considered in the analysis does not couple to the Standard Model and will consist of a N = 10 group, and 2N − 1 gauge bosons with charge g D = 0.8. The scale of the one-loop potential is also T independent. These will result in a speed of sound in the symmetric phase that differs from the one seen in SMEFT.
We show the speed of sound calculated using Eq. 2.47 on the left panel of Fig. 7 for the different levels of diligence. The differences between the levels of diligence in the speed of sound are only minor. We show only the highest diligence curve for the speed of sound in the symmetric phase. For small v/Λ, the speed of sound in the symmetric phase remains constant with a value slightly above the one given in the bag model. This is attributed to the additional degrees of freedom arising from the dark sector. The speed of sound above v/Λ = 2.6 begins to decrease until it reaches a discontinuity near v/Λ = 2.8. It then jumps to c 2 s = 0.336 where it begins to monotonically increase. This discontinuity is a result of the daisy terms in the effective potential. With out the daisy terms, the speed of sound in the symmetric phase would be smoothly connected and monotonically increasing.
The strength of the phase transition is plotted on the right panel Fig. 7. The different levels appear to agree very well with each other with the lowest diligence becoming slightly worse at high v/Λ. For most of the parameter space, the highest diligence has the greatest α because it is computed at the numerically calculated values for c 2 s in the broken phase which results in a amplification compared to the other two levels. This is due to the factor (1 + c −2 s ) in αθ. The error between the modest and highest diligence begins to decrease as v/Λ increases which is related to the speed of sound approaching c 2 s = 1/3. Despite the differences between the different levels, the speed of sound in the broken phase lies between c 2 s ∼ 0.325 − 0.330 and does not contribute a significant source of error to the strength of the phase transition.
The inverse duration of the phase transition is plotted in the left panel of Fig. 8 for detonation. The lowest diligence calculated using Eq. 2.6 consistently over approximates β/H while modest diligence calculated using Eq. 2.13 agrees well with the highest diligence found from the mean bubble separation. There were no minima found in the action for any of the parameters in consideration so the difference between T p calculated using Eq. 2.11 in the modest diligence and T f calculated using Eq. 2.27 in the highest diligence is only minor. The dips near v/Λ > 2.8 in β/H are the result of the shape of S(T )/T which causes the highest error between the modest and highest diligence. This dip also effects the suppression factor due to the lifetime of the source as seen in the right panel of Fig. 8 v/Λ regime has a small β/H and large α which results in a small lifetime of the source τ sw . The modest diligence suppression factor is then an appropriate approximation in this regime. The error in the gravitational wave spectrum is shown in Fig. 9 for deflagration and detonation. For both the lowest and modest diligence, the error remains roughly constant until v/Λ ∼ 2.8 where it exhibits some oscillations. This behavior is related to the dips in Fig. 8. Past v/Λ ∼ 2.8, both the lowest and modest diligence begin to increase. The error in the lowest diligence past this point is dominated by the lack of suppression factor due to the lifetime of the source. The suppression factor remains roughly constant until v/Λ ∼ 2.8 where it begins to approach zero and as a result increases the error. The increasing behavior in the modest diligence is likely due to the separation in β/H between the modest and highest diligence in Fig. 8 and the suppression factor from vorticity and reheating effects in the plasma which are stronger for larger α. The values of the speed of sound in the symmetric and broken phase calculated at T f are not far from the bag model of c 2 s = 1/3 and we do not consider it a strong source of error in the peak gravitational wave energy density spectrum. Overall we notice an error in the peak gravitational wave energy density of 10 1 − 10 3 for lowest diligence and 10 −1 − 10 1 for modest diligence.

xSM
We show in the top left panel of Fig. 10 the speed of sound in the symmetric and broken phase for a scan over the heavy singlet mass in the xSM model while holding all other parameters constant. The speed of sound in the symmetric phase is approximately c 2 s = 1/3 as in the bag model. The speed of sound in broken phase deviates far from the bag model where it approaches zero as m h 2 → 0. The speed of is strongly correlated with the cubic term that arises from the extra scalar who also acquires a tree level vacuum expectation value. The speed of sound can then be suppressed by increasing the b 3 parameter. This strong suppression in the broken phase speed of sound will lead to an amplification in the strength of the phase transition as seen in the top right panel of Fig. 10.
The strength of the phase transition in the highest diligence grows larger compared to the other levels as the singlet gets heavier. This is directly related to the suppression in the speed of sound in the broken phase. There is a minor difference in the lower singlet mass range. The kinetic energy fraction is shown in the bottom panel of Fig. 10. The lowest and modest diligence both overestimate K for the entire range of the parameter space which is not observed in α. This can be attributed to the approximations used in the kinetic energy fraction pre-factor α/(1 + α) used in the peak gravitational wave energy density in Eq. 2.7 and the speed of sound dependence in solving the beyond the bag model hydrodynamic equations. The inverse time duration of the phase transition is plotted in Fig. 11 for the different levels of diligence. The modest diligence is a better approximation than that of the first level for β/H but slightly under-approximates the spectrum for lower mass ranges. The lowest diligence is a poor approximation for β/H for the entire parameter space.
The error in the gravitational wave spectrum compared to the highest diligence for deflagration and detonation is given in Fig. 12. The largest error in the spectrum occurs for the lowest diligence and this is due to the lack of suppression factor for the finite lifetime of the source and the larger uncertainty in β/H. The suppression factor for the modest diligence case is an under-approximation to the finite lifetime of the source particularly in the higher singlet mass regions. Both the lowest and modest diligence receive significant errors from neglecting the beyond the bag contributions to the kinetic energy which over estimates the peak spectrum which also gets worse for higher singlet masses. Overall the range of error in the peak gravitational wave energy density is between [10 2 ∼ 10 3 ] and [10 0 ∼ 10 2 ] for the different levels of diligence. All of the points above the range in m h 2 shown are certainly viable points and may even reach higher levels of error. However, this range in m h 2 is chosen such that all the points remain in either deflagration or detonation for both consistency and the lack of numerical simulations for hybrids.

Mean Bubble Separation vs Inverse Time Duration
The gravitational wave energy density is dependent on determining the mean bubble separation when the phase transition ends at temperature T f . An approximation to the mean bubble separation can be determined by calculating the inverse time duration, β/H, from the first derivative of the action. This calculation is typically only valid when there is a negligible barrier at zero temperature. However, if there is a barrier at tree level, a minimum in the action will develop near T f and the second derivative, β 2 will become relevant while the first derivative will vanish. The bubble nucleation rate can then take on the form where t * is the time when the temperature is near T f and S * = S 3 (T * )/T * . The above result will ultimately lead to a new relation between the mean bubble separation R * and the inverse time duration of the phase transition β.  This subtlety is not usually taken into account and the relation between R * and β that is useful for gravitational wave calculations is simply given by the approximate formula where β is related to the first derivative of the action regardless of the presence of a barrier. Out of the models we consider, SMEFT and xSM can acquire tree level barriers that result in a minimum in the action. The lowest and modest diligence results presented here assume Eq. 4.3 always hold which can result in significant errors for these two models. Furthermore, the percolation temperature at which β/H is estimated is a function of β/H, Eq. 2.13, which can also acquire error if the barrier is not sufficiently taken care of. The highest diligence results can side-step these issues by numerically calculating the mean bubble separation from the number bubble density which is independent of any assumptions about the curvature of the action, i.e  Fig. 13 where the left figure corresponds to SMEFT and right figure corresponds to xSM. The solid lines represent the proper mean bubble separation calculated at T f . The dotted and dashed lines correspond to the mean bubble separation calculated first from β at T f and T p respectively. We denote T p to refer to the estimation given in Eq. 2.11. Below M = 600 GeV in SMEFT, the action acquires a minimum as a result of the tree level barrier at zero temperature which causes β(T p ) to significantly overapproximate HR * . As mentioned previously, this is a result of the underlying assumptions in approximating both T p and β which ignore the barrier. The mean bubble separation calculated from β(T f ) performs better than β(T p ) in this regime with nearly identical HR * predictions compared to n . This is largely due to T f being independent of any assumptions on the action. The xSM model consists of a second scalar and several parameters which when varied may induce either first step or second step phase transitions. The bench marks chosen involve scanning of the heavy singlet mass while holding the other model parameters fixed. All of the points resulted in a one step phase transition along with no minimum in the action. On the right of Fig. 13, we see that all three methods resulted in a roughly consistent approximation of HR * with slightly better performance from β(T p ) for large m h 2 . This can be attributed to the lack of minimum in the action observed in the parameter space. We reserve a further analysis of the mean bubble separation in xSM for future work.

Summary of Results
The previous results involved fixing certain characteristics associated with each of the outlined levels of diligence. In this section, we will fix all of the quantities as high diligence while varying the level of a single quantity to determine its impact on the error of ∆Ω/Ω. Table 1 shows the range of error we observe associated with varying the level of diligence in the calculation of the transition temperature, mean bubble separation, fluid velocity, finite lifetime of the source, and vorticity effects. The base level of comparison will use Ω GW calculated using Eq. 2.61 which assumes the transition temperature is at T f and includes beyond the bag effects and the suppression factors due to the finite lifetime of the source and reheating effects in the plasma. We will now proceed to describe how the range of errors are calculated.
The transition temperatures used for the different levels were T n (2.4) and T p (2.11. The frequency independent Ω GW is now calculated at high diligence using Eq. 2.61 for both the lowest and modest diligence. This is to show how Ω GW can change purely by the temperature at which  the transition is assumed to take place. The lowest diligence will use T n to calculate ∆Ω/Ω while the modest diligence will use T p . The base level comparison is Ω GW at T f which corresponds to the previously defined high diligence. Varying the transition temperature leads to an error of (10 −1 − 10 0 ) and (10 −4 − 10 1 ) for lowest and modest diligence respectively. The modest diligence can experience a larger error than the lowest diligence and this is due to the result of the strong super-cooling observed in SMEFT when M 600. The approximations used in calculating T p break down when a minimum develops in the action and results in the 10 1 peak in the error for modest diligence. The error in the lowest diligence results in an enhancement in the spectrum which is attributed to T n > T f . The modest diligence experienced both enhancement and suppression which is due to T p being much closer to T f . The lowest diligence primarily had Ω low GW < Ω high GW and modest diligence had Ω med GW > Ω high GW . For these reasons, we conclude that the type of error due to the transition temperature is random and dependent on the underlying model.
The estimation of the error due to the mean bubble separation will involve calculating R * H from the β/H at T n for the lowest diligence and R * H from n b at T f for modest diligence. We use T f for the modest diligence in determining the relevant quantities in Ω GW in this case to minimize error which may arise from using the T p approximation. All quantities in Ω GW are calculated in high diligence at T n and T f for lowest and modest diligence respectively. The lowest diligence exhibits the largest error with a range of (10 −1 − 10 0 ) while modest diligence has the range (10 −3 − 10 −1 ). The error in modest diligence observed in the table is only due to the approximation of the mean bubble separation from the inverse time duration but it is expected to be higher if T p is used as opposed to T f which helps to correct the error. Both the lowest and modest diligence had mostly Ω low,med GW < Ω high GW with modest diligence having a couple benchmarks with Ω med GW > Ω high GW . We denote this type of error as predominately suppression. The error estimate from the fluid velocity involves comparing the fits for kappa given in Appendix A to solving the hydrodynamic profiles numerically. The fluid velocity is related to the kappa through the kinetic energy fraction K in Eq. (2.10,2.55,2.57). The lowest diligence calculates Ω GW at T f in the highest diligence with the replacement that K and U f are now calculated using the fits to κ and the bag calculation for α. The modest diligence is the same as the lowest diligence except that κ is calculated using the hydrodynamic profiles with c 2 s = 1/3 in the bag model. The error associated with the different treatments of the fluid velocity is (10 −2 − 10 0 ) for lowest diligence and (10 −3 − 10 0 ) for modest diligence. This represents the amount of error that one might expect in these models when performing precise calculations of Ω GW but without taking into consideration the beyond the bag treatment of the speed of sound in the plasma. The type of error we observe for the fluid velocity is random.
To determine the impact of the suppression factor due to the finite lifetime of the source has on the error, we compare Ω GW calculated in Eq. 2.61 with out Υ for the lowest diligence and with the replacement Υ → τ SW H corresponding to Eq. 2.14 for the modest diligence. All quantities are evaluated at T f . Modest diligence will also contain the suppression to U f due to vorticity and reheating effects in the plasma. Note that this suppression is less dramatic than what one might naively expect from ref [106], as the suppression in the fluid velocity results in a longer lifetime for the soundwaves. For the range of models we consider, the error for modest diligence is in the range (10 −3 − 10 −1 ) and represents the validity of the approximation used in Υ. The error in the lowest diligence experiences the highest error with a range of (10 1 − 10 3 ). For all of the models, Ω low,med GW > Ω high GW . This type of error is an enhancement. The last row in Table 1 corresponds to the error in Ω GW calculated using Eq. 2.61 without suppression factors arising from vorticity and reheating effects in the plasma. This is compared to the full suppression in highest diligence which uses U f,max in the lifetime of the source as well. The range of error we observe is in the range of 10 −1 − 10 1 . Neglecting U f,max in the suppression factor will contribute at most an error of 0.62. The lowest diligence experienced Ω low GW < Ω high GW for the all of the models. The modest diligence experienced mostly random error. The primary focus should be on modest diligence so we denote this type of error as random.

Conclusion
In this work we have examined the cost of various short-cuts and approximations used in the literature when predicting the gravitational wave spectrum generated by a cosmological first order phase transition. Even in the case where some modest diligence has been used in the calculation, we found the cost to often be comparable to problems in finite temperature QFT such as the scale dependence that arises from the break down of perturbation theory as well as the gauge dependence. Assuming detonations, the dominant cost in cases where there is a fair amount of super-cooling is from poor estimates of the percolation temperature in Eq. 2.11. The poor estimate of the percolation temperature has a knock on effect in enhancing the errors that arise from using the bag model and an estimate for the suppression factor. In the case where there is no tree-level barrier delaying the nucleation of the phase transition, the dominant error is due to the bag model approximation.
Although the errors are often as large as finite temperature QFT errors, they are arguably easier to reduce. At present, all of these errors can be handled except for the reheating and vorticity effects where we had to rely on interpolations. High diligence calculations for multiple models were considerably more tractable than the two loop calculations required to bring scale dependence at finite temperature under control [93]. We recommend future phenomenological calculations of gravitational wave signals from primordial phase transitions at the very least take the level of theoretical uncertainties into consideration.

Acknowledgements
HG and KS are supported by DOE Grant desc0009956. GW was supported by World Premier International Research Center Initiative (WPI), MEXT, Japan.

A Kinetic Energy Efficiency Coefficient
The kinetic energy efficiency coefficient may be solved by integrating over the enthalpy and velocity profiles around a single bubble, where is the bag constant and ξ = r/t is a self similar coordinate in terms of the distance form the bubble center r and the time since nucleation t. The fits to κ are provided in [110] and are valid in the range 10 −3 < α < 10 to a precision of 10%. The fits are found by splitting the parameter space of v w into three regions and four boundary conditions. The boundary conditions are where v J is the Jouguet velocity and c s is the speed of sound. Subsonic deflagrations in the region v w c s have a kinetic energy cofficient approximated by and detonations in the region v w > v J by Supersonic deflagrations, hybrid, in the region c s v w v J can be approximated by B Toy Model

B.1 Toy Effective Potential
A general free energy density of a single scalar field, φ, under a high temperature expansion can be written in the form where ∆V is added to the potential to cancel out the zero temperature minimum such that f (φ min , 0) = 0. The Standard Model effective potential was considered in Ref. [125]. We require D > 0, E > 0, λ > 0 to ensure symmetry is broken at low temperature and generate a barrier between the symmetric and broken phase. The vacuum terms are not necessary for determining the phase transition structure of the model, however, they are necessary for determining the temperature dependence of the speed of sound. The structure of the potential along with the constraints on the parameters allows for simple analytical forms for the minima as a function of temperature. The minimum is found by minimizing Eq. B.1 with respect to the scalar field which results in where the '+' sign gives the local minimum. When T is large, the global minimum will sit at the origin with φ min = 0. As T decreases, a second minimum will develop at This minimum will eventually become degenerate with minimum at the origin at the critical temperature when The Euclidean action of a bounce configuration, S 3 , will start from infinity at T = T c and decrease with temperature. There is an analytical form for the action given by where σ = λM 2 /(2E 2 T 2 ) controls the overall shape of the potential [125]. The critical temperature and the action are necessary to determine the dynamics of the phase transition and calculate the relevant transition temperatures such as T n ,T p , and T f and the mean bubble separation R * (v w , β). where the pressure in the broken phase has additional dependence on temperature arising from φ min . The speed of sound may be found from the pressure using Eq. 2.47 in both the symmetric and broken phase. The temperature dependence from the minimum of the scalar field will result in a speed sound that is function of the model parameters and its form will depend on the overall shape of the potential.

B.2 Results for toy model
Here we show the different levels of diligence in calculating the thermal parameters and the gravitation wave spectrum in the toy model. The analysis involves individual scans over the different model parameters (E, D, λ, T 0 ) while holding the others fixed. A full analysis of the toy model should involve a randomized scan over all of the parameters but we perform the scan this way in hopes to see any trends in varying the different model parameters. In Eq. B.4, the critical temperature is a function of all four model parameters. For this reason, T C will be used as a basis for each scan. The first step in the beyond the bag calculations is to compute the speed of sound in the symmetric and broken phase. For the toy model, we only consider detonation, v w = 0.92, where the speed of sound in the symmetric phase is set to c 2 s = 1/3 and the degrees of freedom consist of only the standard model sector. The speed of sound in the broken phase may be found through Eq. 2.47. The transition temperatures for the different levels of diligence are T n (2.4), T p (2.11), and T f (2.27). Example calculations for various phase transition quantities used in the high diligence calculations such as the false vacuum fraction, mean bubble separation, lifetime of the source, and T f in the toy model may be found in Ref. [103].
In Fig. 14, we calculate the speed of sound in the broken phase for each level of diligence. The gray dashed line corresponds to c 2 s = 1/3. This involves first calculating the speed of sound as a function temperature using Eq. 2.47 and then evaluating it at (T n , T p , T f ) computed in the different levels of diligence. We note that in computing the strength of the phase transition only the highest diligence level will involve this calculation. This is merely to show the level of variance in computing the speed of sound at different temperature stages. For specific values and range chosen, there is only minor change to the speed of sound computed in the different levels however how much variance is present is strongly model dependent. We do notice that the speed of sound can have a significant deviation away from c 2 s = 1/3 in the bag model. The strongest deviation is caused by varying the barrier term, E, and the quadratic multiplicity term, D as seen in the green and purple curves. The speed of sound can go as low as c 2 s ∼ 0.22 and as high as c 2 s ∼ 0.36. Varying the zero temperature mass term, T 0 , did not have any noticeable impact on the speed of sound while the quartic coupling term, λ, had a mild impact on the speed of sound. This is likely due to the temperature independence of the terms that involve T 0 and λ. The parameters D and E on the other hand, multiply T 2 and T respectively and will result in a change in the temperature dependence. The speed of sound in the broken phase is related to the temperature derivatives of the pressure which is evaluated at p b (T ) = V eff (φ min (T ), T ) and hence D and E will impact the minimum at finite temperature. The smallest speed of sound in the broken phase corresponds to small E and large D.
We show in Fig. 15 the phase transition strength computed in the different level of diligence (left) and the comparison between α θ computed in the bag model versus αθ computed in the beyond the bag model at T f for both quantities (right). We see in the left figure that going higher in the level of diligence results in an increase in the phase transition strength compared to lowest diligence. This can be attributed to more vacuum energy being released at T p compared to T n . On the right, to better compare the difference between the bag model and the beyond the bag model, we compute the ratio of α θ and αθ computed at the same temperature, T f . For T c < 100, α θ is less than αθ which is the result of c 2 s < 1/3 as seen in Fig. 14. This has to do with the (1 + c −2 s ) factor in αθ. When T c > 100, we see that the opposite is true when c 2 s > 1/3. Similarly, the largest deviations are due to the parameters D and E.
The error in the gravitational wave spectrum of the toy model for different scans in the model parameters is shown in the left of Fig. 16. The lowest and modest diligence peak gravitational wave energy density Ω GW is calculated using Eq. 2.7 and Eq. 2.15 respectively. The comparison in error is computed with respect to the highest diligence in Eq. 2.61. The lowest diligence level has error in the range ∆Ω/Ω ∼ [10 1 , 10 3 ] for all parameter scans. The modest diligence level is closest to the highest diligence with error in the range ∆Ω/Ω ∼ [10 0 , 10 1 ] for the different scans.
The highest error occurs for T C ∼ 100. This is related to beyond the bag effects which exhibited suppression for the scans in (E, D).