title: Gravitational Waves From Dark Sectors, Oscillating Inflatons, and Mass Boosted Dark Matter
authors: Amit Bhoonah Joseph Bramante Simran Nerval Ningqiang Song
text: 
  Learning Representations of Endoscopic Videos to Detect Tool Presence Without Supervision
      David Z. Li 1Masaru Ishii 2 Russell H. Taylor 1 Gregory D. Hager 1Ayushi Sinha 3
      December 30, 2023
  =========================================================================================
  
  
  
  
  Gravitational wave signatures from dynamical scalar field configurations provide a compelling observational window on the early universe. Here we identify intriguing connections between dark matter and scalars fields that emit gravitational waves, either through a first order phase transition or oscillating after inflation. To study gravitational waves from first order phase transitions, we investigate a simplified model consisting of a heavy scalar coupled to a vector and fermion field. We then compute gravitational wave spectra sourced by inflaton field configurations oscillating after E-Model and T-Model inflation. Some of these gravitational wave signatures can be uncovered by the future Big Bang Observatory, although in general we find that MHz-GHz frequency gravitational wave sensitivity will be critical for discovering the heaviest dark sectors. Intriguingly, we find that scalars undergoing phase transitions, along with E-Model and T-Model potentials, can impel a late-time dark matter mass boost and generate up to Planck mass dark matter. For phase transitions and oscillating inflatons, the largest dark matter mass boosts correspond to higher amplitude stochastic gravitational wave backgrounds.
  
  
  
  
  
  § INTRODUCTION
  
  Despite a wealth of cosmological data collected in the last two decades, almost nothing is certain about the dynamics of our universe at temperatures above an MeV, when nuclei were synthesized. Primordial fluctuations evident in the cosmic microwave background (CMB) and the large scale distribution of galaxies are both consistent with a period of inflation and subsequent matter or radiation dominated expansion <cit.>. However, both of these primordial structures become apparent at temperatures below an MeV, and both are consistent with a spectrum of initial scalar perturbations. With only scalar primordial perturbations, it is not possible to pinpoint the energy density of the universe either during or after inflation. In contrast, primordial gravitational waves are eminently sensitive to the energy density of the universe at the time of their creation <cit.>. It follows that the discovery of primordial gravitational waves, either as tensor perturbations to the CMB or as a stochastic gravitational wave background (SGWB) sourced by scalar field dynamics <cit.>, could provide a fascinating new window on the dark history of our universe at energy densities above an MeV. In this work we will consider two separate sources of gravitational waves in the early universe: first order phase transitions and inflatons oscillating after inflation. 
  
  Because they would provide a stunning new view of the early universe, it is natural to consider whether primordial gravitational waves can also inform us about the formation of dark matter, which is inextricably linked to the state of the universe at temperatures above an MeV. Although we have accumulated extensive gravitational evidence for dark matter, the couplings, mass, and mechanism by which dark matter was produced in the early universe remain a mystery.  While much work has scrutinized production of thermally equilibrated weakly interacting massive particle (WIMP) dark matter in a radiation dominated universe, for dark matter masses ranging up to 10 TeV <cit.>, simple physical mechanisms that result in the observed relic abundance of dark matter with a mass ranging up to the Planck mass are still being uncovered. A canonical example of superheavy dark matter produced out of thermal equilibrium is the WIMPzilla, gravitationally produced at the end of cosmic inflation <cit.>. However, this nonthermal process is limited to producing dark matter with a mass below ∼ 10^14 GeV, since gravitational production is exponentially suppressed for masses exceeding the Hubble constant at the end of inflation, which itself has been limited to H_e ≲ 10^14 GeV by the observed tensor-to-scalar ratio of primordial perturbations <cit.>. On the other hand, dark matter heavier than 10^14 GeV is particularly compelling, since it may be connected to Grand Unified Theories (GUT) at mass scales above 10^14 GeV <cit.>, and could soon be discovered at underground experiments using multiscatter and other techniques <cit.>. A number of papers have considered how suppressed gravitational production of dark matter heavier than 10^14 GeV can be supplemented by production through the decay of the inflaton or out of equilibrium “freeze-in" processes during the radiation dominated epoch following inflation <cit.>. 
  
  In this paper we consider a simple mechanism for gravitational production of dark matter with a mass well above 10^14 GeV. The mechanism operates by coupling dark matter weakly to a scalar field, which undergoes a phase transition that boosts the mass of the dark matter after its initial abundance is set by gravitational processes at the end of inflation. In contrast, some previous work has considered how mass shifts and phase transitions can alter the relic abundance of thermally produced dark matter <cit.>. Here we find distinct features and compelling gravitational wave signatures associated with a scalar field that can boost the mass of dark matter with a nonthermal abundance fixed by the end of inflation. In particular, we find that parameters associated with a strong first order phase transition that source a substantial SGWB are associated with a large dark matter mass boost. We also find that a large dark matter mass associated with certain inflaton potentials also corresponds to a large source of gravitational waves generated by the inflaton oscillating in its potential.
  
  The remainder of this paper is organized as follows. In Section <ref>, we detail a simplified model consisting of a scalar with a quartic self-coupling generated at one-loop order through coupling to a vector and fermion field. First order phase transitions and associated gravitational wave signatures of this potential are explored in Section <ref>. Stochastic gravitational waves sourced by E-Model and T-Model inflation are computed using lattice simulations in Section <ref>. Section <ref> details the cosmology of very heavy dark matter that receives a late-time mass boost from a scalar field. In Section <ref>, we conclude with a summary of observational prospects for gravitational wave signatures. In Appendices <ref> and <ref> we detail higher order corrections and running couplings for the thermal potential studied in Section <ref>.
  
  
  
  § PHASE TRANSITIONS FROM A SIMPLIFIED MODEL
  
  
  
  In this section we introduce a simplified model that will be useful for investigating phase transitions and dark matter mass boosts in the early universe. Our model consists of a scalar field with a small tree-level quartic self-coupling, coupled to a fermion and gauge field. More specifically, for the parameter space we will find most interesting, the scalar field's quartic self-coupling is predominantly loop-induced, through coupling to fermion and gauge fields. To explore the salient features of such a model, we consider the case that the gauge group is Abelian - the non-Abelian generalization is relatively straightforward. 
  
  Then in total the simplified dark sector consists of a scalar φ charged under a U_D(1) group, a vector A_μ, and a fermion χ. For the moment, we do not consider the couplings of this sector to the Standard Model. The Lagrangian can be written as a sum of tree level terms, one loop contributions, ℒ_1, and counter-terms which we omit for clarity,
  
      ℒ = -1/4F_μνF^μν + iχγ^μ∂_μχ + 1/2D_μφ(D^μφ)^† - V_0(φ) - iy χχφ + ℒ_1,
  
  where D_μφ = ∂_μφ + igA_μφ and V_0(φ) is the tree level scalar potential. Here we assume that the fermion and the vector boson acquire masses through spontaneous symmetry breaking when the scalar develops a non-zero vacuum expectation value. In future sections when it is cosmologically relevant, we will also add an “initial” mass term for the fermion of the form -m_fiχχ, where χ will be our candidate dark matter field.
  
  Here we assume that the scalar has no self-coupling or couplings to other scalars. Our motivation for doing so is twofold. First, such scalars can arise naturally in theories containing a scalar field with a broken shift symmetry <cit.>. Second, such a scalar gains a large vacuum expectation value (vev) at the minimum of its potential. This second feature can be understood by considering a scalar potential with self coupling λ,
  
      V_0(φ)     =  λ/4!φ^4 - m^2_S/2φ^2.
  
  Assuming both m^2_S and λ are greater than zero, the minimum of this potential is given by
  
      < φ> ∼m_S/√(λ).
  
  Clearly, a large vacuum expectation value can be generated either through a large scalar mass or a small self-coupling. Of course if the scalar is coupled to other fields, its self-coupling cannot be arbitrarily small, as it will receive radiative corrections from coupling to other fields. In our setup, a small loop-induced effective quartic coupling λ_eff will arise from the scalar field's couplings to a vector boson and a fermion. In the case that the effective self-coupling is generated at loop level by the fermion coupling, λ_eff∝ y^4. In this case, if y is small, λ_eff will be very small, resulting in a very large scalar vev. 
  
  Through its Yukawa coupling to the scalar field, the fermion will obtain a mass. However, it is interesting to note that for the setup described above, the Yukawa mass for the fermion can actually increase as the Yukawa coupling shrinks, since in the absence of other couplings, the Yukawa term would scale as m_χ∼ y m_S/ √(λ)∼ m_S/y. We explore this feature in more detail in Section <ref>. One might also consider a fine-tuned scenario, where both the Yukawa coupling and final scalar vacuum expectation value are large. As we will see in Eq. (<ref>), the loop-level fermion and vector contributions to the scalar's self-coupling carry opposite signs and can be tuned to cancel, thereby generating an extremely small λ_eff. In this fine-tuned case, the size of the vacuum expectation value is then an interplay between the scalar mass, the size of the Yukawa and gauge couplings, and the degree of tuning between them. 
  
  Spontaneous symmetry breaking in this model is realised through radiative corrections, similar to the Coleman Weinberg (CW) <cit.> mechanism. An important difference to highlight is that while CW assumes the scalar is massless but self-interacting at tree level, we assume that the scalar is massive and its self-interaction parameterized by λ_eff are generated at loop-level through coupling to other fields,
  
      V(φ) =λ_eff/4!φ^4 - m^2_S/2φ^2.
  
  A comment about the potential can be made here. Since λ_eff is zero at tree level, the potential may at first appear tachyonic. Specifically, in Eq. (<ref>), if m_S^2 > 0 the tree level potential is unbounded from below. However, as we show after computing the one loop effective potential, quantum effects stabilise the potential. For the parameters we will consider hereafter, the potential is bounded from below and the physical mass of the scalar - given by the second derivative of the potential evaluated at its minimum - is positive, meaning that the system has no tachyonic degrees of freedom. Here as in the CW model the scalar φ obtains a vacuum expectation value which depends on its coupling to other fields.
  
  The method we will use for studying loop effects on spontaneous symmetry breaking is effective potentials. Here we take care to define what we mean by “effective potentials" since the term can appear in a few different contexts in the literature. The first arises when one is mainly looking to study the low energy regime of a theory by integrating out heavy degrees of freedom whose effects are unimportant at energies far below the high energy cutoff of the theory. The second context, which is the one considered in this work, consists of integrating out quantum fluctuations of fields in the presence of classical sources <cit.>. While there are different ways of doing this, the most widely used technique for one loop computations is the background field method, where all the fields in the action are split into a classical background piece (c) and a quantum piece (q). In our case this means shifting the scalar, vector, and fermionic fields as
  
      φ   →φ_c + φ_q
      
          A^μ   → A^μ_q
      χ   →χ_q
  
  respectively. Here, as is in prior literature <cit.>, we turn on a background field for the scalar. When we then integrate out the pieces quadratic in the quantum fields, we are left with an effective potential that depends only on the classical background field φ_c, which is a space-time independent constant (in general it need not be) and encodes all the one particle irreducible graphs of the original theory at one loop. We have been quite sparse in our review the background field method as there are already many helpful references on this topic in the literature. For examples of background field methods applied to potentials significant to cosmology like those considered here, we refer interested readers to <cit.> for a review. 
  
  The one-loop effective potential, which we evaluate in the Landau gauge, contains a divergent piece and an “effective quartic" interaction - effective because it is generated from loop-level corrections and so takes the form coupling multiplied by φ^4 - and logarithmic pieces. In more detail, we take the scalar in our theory to be complex and so there are actually two real scalar fields, φ_1 and φ_2. However, the effective potential in the Landau gauge only depends on φ = √(φ^2_1 + φ^2_2). This introduces factors of two in our loop calculations which we can absorb in normalizations for the scalar field. In the interest of clarity, we therefore use φ throughout our calculations. If we had a tree level quartic coupling, we would also have another divergent piece from the λφ^2 term that would require renormalization using the mass counter-term. However, in our current setup this mass renormalization occurs at second order in perturbation theory and will not be considered in this work. Here we use the following renormalization conditions: 
  
  
      V^''(φ) |_φ = 0    =-  m_S^2
      
       V^''''(φ) |_φ = M    = 9/16π^2(4y^4 - 3g^4)
  
  
  where primes denote derivatives with respect to φ. These are similar to the renormalization conditions in CW, but with the loop-generated quartic substituted for a tree-level quartic in (<ref>). 
  
  With these conditions, the renormalized one-loop effective potential at scale M = < φ> ≡ν is
  
      V(φ) = -m^2_S/2φ^2 + 17/192π^2(4y^4 - 3g^4)φ^4 - 4y^4 - 3g^4/64π^2log[φ^2/ν^2]φ^4.
  
  The vacuum expectation value at one loop is given by the zero of the first derivative of the potential at scale φ = ν,
  
      ν^2  = 96π^2/31m_S^2/(4y^4 - 3g^4),
  
  where we emphasize that this expression is valid to one loop order at φ= ν, at which field point logarithmic terms like those in Eq. (<ref>) are absent for the one loop expression.
  Imposing that spontaneous symmetry breaking happens, ν^2 > 0, requires 4y^4 > 3g^4. Of course, this spontaneous symmetry breaking will also generate a mass for the fermion m_f = yν and the vector field m_V = g ν. There are two ways in which this would result in a large vev for φ and corresponding large masses for the fermion and vector boson. First, if the Yukawa and gauge couplings are small, g,y ≪ 1, this can result in a large vacuum expectation value for φ. In lower regions of Figure <ref>, we show parameter space where a large fermion mass results from small Yukawa and gauge couplings. Secondly, even if these couplings are order unity, g,y ∼ 1, in the case that 4y^4 ≈ 3 g^4, the fermionic and vector contributions may cancel to some degree, again leading to heavy fermions and vector bosons but a comparatively light scalar. However, we note that if a large ν is generated via a cancellation between the gauge and Yukawa couplings, some additional consideration may be required. The results here have been computed at one loop level, and so determining cancellation between couplings to a certain precision would require a more precise two loop treatment. While a two loop calculation is beyond the scope of this work, we can gain some insights into the effect of higher order loops using parametric estimates. The leading two loop corrections to the effective potential from the gauge and Yukawa sectors are 𝒪(g^6) and 𝒪(y^6) respectively, and as long they are both perturbatively smaller than the leading terms 𝒪(g^4) and 𝒪(y^4), higher order effects should be subdominant. In our setup we anticipate that two loop corrections will not significantly affect our results, so long as g/y ≲ 1.07, since for the strongest couplings we consider g =1.07, y= 1, the leading terms correcting the one loop λ_eff result should scale like g^2/16 π^2∼y^2/16 π^2∼ 0.007, compared to the larger one loop term 4y^4 - 3g^4≈ 0.07. In future work it might be interesting to explicitly calculate the effective potential at two loop order in perturbation theory. 
  
  
  A few additional remarks are in order. We define the physical mass of the scalar, M_S as the second derivative of the potential at φ = ν,
  
      M^2_S  = V^''(φ) |_φ = ν = - m^2_S + 27/32π^2(4y^4 - 3g^4)ν^2
  
  
  It can be verified by substituting in the expression for ν from (<ref>) that the square of this scalar physical mass, M^2_S, is positive.  
  
  For completeness, we also comment on gauge invariance for the scalar potential. A detailed discussion of gauge dependence in potentials like the one we consider can be found in Ref. <cit.>. In our derivation, we have performed our calculations in Landau gauge where the gauge fixing parameter (usually called ξ in the literature) does not appear. For our potential we anticipate that as expected (and proven for the CW model in <cit.>), once couplings are appropriately fixed at renormalization scale M, our results are gauge invariant. 
  
  
  
  
  § FIRST ORDER PHASE TRANSITIONS & GRAVITATIONAL WAVES
  
  
  We now turn to the cosmological evolution of our simplified model in a radiation dominated background. When thermal effects are included, the effective potential (<ref>) acquires temperature dependent pieces
  
      V(φ, T) = -m^2_S/2φ^2 + 17/192π^2(4y^4 - 3g^4)φ^4 - (4y^4 - 3g^4)/64π^2log[φ^2/ν^2]φ^4
          + T^4/2π^2∫_0^∞ dx  3x^2log[1-e^-√(x^2 + g^2φ^2/T^2)] - 4 x^2log[1+e^-√(x^2 + y^2φ^2/T^2)] .
    
  The integrals above have no known analytic forms and must be evaluated numerically. However, it is possible to expand this potential in the high limit, where yφ/T≪ 1 and gφ/T≪ 1. The effective potential then takes the form
  
      V(φ, T) = D(T^2 - T_0^2)φ^2 -E T φ^3 + λ(T)/4φ^4 ,
    
  where we define the following terms 
  
  
      E     = 3 g^3/4π , 
      
      D      = (4y^2 + 3g^2)/24 , 
      
      T_0    = √(m_S^2/2D) =  √(12 m_S^2/(4y^2 + 3g^2)) ,
      λ(T) = 17/48π^2(4y^4 - 3g^4) + 1/16π^2   (4y^4log[y^2ν^2/a_F T^2] - 3g^4log[g^2ν^2/a_B T^2]) ,
  
  
  which can be compared to similar terms defined in <cit.>.
  Here a few numerical constants have been added, log a_F≈ 1.14 and log a_B≈ 3.91. These constants are obtained through matching numerical solutions of Eq. (<ref>) to Eq. (<ref>). The temperature dependent extrema of this potential, not including the one at zero, are given by
  
      ν_± (T) = 3ET/2λ±1/2λ√((9E^2 - 8λ D)T^2 + 8λ D T_0^2) .
  
  We note that the cubic term ∝ ET is responsible for creating a barrier between the two extrema, which is typical of a first order phase transition. 
  
  Before detailing the phase transitions facilitated by this potential, let us first review the qualitative features of a first order phase transition. An analogy can be made to boiling water in a container. As the temperature inside the container is raised to the boiling point of water, the vapour phase is energetically favoured; the liquid phase becomes the false vacuum and the vapour phase becomes the true vacuum. A growing vapour bubble converts regions of liquid into vapour (in other words, converts regions of the false vacuum state into the true vacuum state), and soon the entire container is filled with vapour, the new true vacuum. For cosmological phase transitions where the universe is undergoing expansion and dilution, matters proceed in the opposite temperature direction. 
  
  In Figure <ref>, we show the shape our potential takes during a first order phase transition. At very early times (high temperatures) the global minimum of the scalar field undergoing the phase transition is at a field value of zero; the universe is said to be in the symmetric phase. A local inflection point forms away from zero at a temperature T_1. As the temperature drops below T_1, a local maximum and local minimum form, given by the ν_-,+ solution of (<ref>) respectively. At this temperature the energy of the ν_+ vacuum state is higher than the energy of the minimum at the origin - the ν_+ minimum is energetically disfavoured. However, as the temperature decreases further the local minimum at ν_+ descends until we reach a temperature T_C at which it is energetically degenerate with the minimum at the origin. For temperatures T < T_C, the ν_+ minimum is energetically favoured. 
  
  
  By definition, the first order phase transition begins at T < T_C. Early on during the phase transition when T ≲ T_C, there is a minimum at ν_+. At this point the universe is still in the vacuum state located at the origin. As the phase transition progresses, bubbles of the (now) lower energy ν_+ phase nucleate and grow. These growing bubbles expand, percolate, and collide to create larger and larger regions in the ν_+ phase. If we imagine for a moment that the universe is static, then the phase transition completes when the entire universe is in the ν_+ phase - all the false vacuum regions have been converted to the true vacuum state. However in an expanding universe two important factors need to be considered: (1) the bubble nucleation and growth process needs to be fast compared to Hubble expansion and (2) an expanding universe cools down, which affects the nature of the phase transition. For the phase transition to be first order and proceed through bubble nucleation, the potential needs to have two minima: in this case at the origin and at ν_+. However, at some temperature T_0, the minimum at the origin vanishes, and so the phase transition is no longer first order for temperatures below T_0. Hence for our purposes, it is useful to determine whether the phase transition completes quickly relative to the rate of Hubble expansion, before the universe reaches a temperature T_0. An important temperature in this regard is the nucleation temperature T_N, at which the bubble nucleation rate per unit volume per unit time is comparable to the Hubble rate of the universe. At T_N the phase transition is efficient and the entire universe is converted to the ν_+ phase. If T_N > T_0, it is safe to assume the entire transition is first order – this criterion for first order phase transitions was first set out in <cit.>. The nucleation temperature T_N will also be the characteristic temperature for the production of gravitational waves via bubble collisions and associated acoustic waves generated in the primordial plasma.
  
  We now begin a quantitative treatment of first order phase transitions caused by the scalar potential in Eq. (<ref>). First, we analyze the potential at the origin (φ = 0) to validate some of our prior description of the potential for temperatures above and below T_0. Note that as depicted in Figure <ref>, T=T_0 occurs well after any first order phase transition has completed. Because the extrema of the potential at various temperatures determine the characteristics of its phase transition, we consider the second derivative of Eq. (<ref>),
  
      ∂^2V(φ,T)/∂φ^2 = 2D(T^2 - T_0^2) - 6E T φ + 3λφ^2.
  
  It can be seen from (<ref>) that the extremum at φ = 0 is a minimum for temperatures T > T_0, a saddle point at T = T_0, and a maximum when T < T_0. At T = T_0, there is only one non-zero extremum, the minumum at ν_+(T_0) = 3ET_0/λ(T_0)≡ν_T_0.
  
  Now we turn to the dynamics of the potential around temperature T_1, which is defined by the formation of an inflection point at ν_+. Until the universe reaches this temperature, there is only one minimum at the origin. At T_1 the term under the square root in (<ref>) vanishes and the potential forms an extremum at φ = ν_+.
  The value of this extremum is found using (<ref>), ν_T_1≡3ET_1/2λ(T_1). Provided the temperature dependent quartic interaction term λ(T_1) is positive, this is an inflection point that turns into a minimum at temperatures below T_1. In our case, λ(T_1) < 0 only occurs at extremely high temperatures above the Plank scale, well beyond the region of validity for our calculations. Initially, V(ν(T_1), T_1) > 0 and the minimum at zero is energetically favoured since V(0, T_1) = 0. However, as the universe expands and cools, we reach some temperature T = T_C where V(ν(T_C), T_C) = 0. This condition simplifies to
  
      T^2[1 - E^2/λ D]  |_T=T_C = T_0^2.
  
  Bubbles of true vacuum form at T_C, since at this temperature a new true minimum appears at ν_+. Below T_C, bubbles may nucleate. To compute the volumetric rate of bubble nucleation we use a Euclidean action treatment <cit.>,
  
      Γ = A  exp (-S_E),
  
  where A is a dimensionful 𝒪(T^4) term and S_E is the Euclidean action for tunneling from φ = 0 to the new minimum. Near the time of bubble nucleation t_N, which corresponds to a nucleation temperature T_N, the Euclidean action S_E in (<ref>) can be Taylor expanded around time t_N, 
  
      S_E   ≃ S_E(t_N) - (t - t_N) β,
  
  where β = - dS_E/dt|_t = t_N can be understood as the inverse duration of the phase transition. Furthermore, the Euclidean action can be factorized into a spatial and temporal piece, S_E = S_3/T, such that
  
      S_3 = 4π∫_0^∞ dr  r^2[ 1/2(dφ/dr)^2 - V(φ, T) ],
  
  is the three dimensional Euclidean action. By solving the Euclidean action for (<ref>), an expression can be found for S_3/T of the form <cit.>
  
      S_3/T   ≃4.85/E^2[2D(T^2 -T_0^2)]^3/2/T^3 f(x)
       ,
      
      
      f(x)     = 1 + x/4(1 + 2.4/1 - x + 0.26/(1 - x)^2)  , 
      
      x      = λ D(T^2 -T_0^2)/E^2T^2 .
  
  With these expressions, we can now give a more precise definition of the nucleation temperature T_N.  Because we are most interested in the onset of the phase transition, we want to find the temperature at which one bubble is nucleated on average per Hubble volume (V_H≡ H^-3); this nucleation temperature T_N is defined by
  
      ∫_t_C^t_N dt Γ/H^3 = ∫_T_N^T_C dT/TΓ/H^4 = 1,
  
  where for the first equality we have assumed that the universe is expanding adiabatically. To simplify this expression, we can use the fact that in a radiation dominated universe the Hubble parameter is given by H^2 = ρ/3M_Pl^2, where ρ is the energy density and M_Pl is the reduced Plank mass. We assume that at the time of the phase transition, only standard model particles in the plasma are relativistic, yielding ρ = π^2/30× g_SMT^4 where g_SM≈ 106.75 is the total number of Standard Model effective relativistic degrees of freedom. Putting everything together, the definition of the nucleation temperature becomes
  
      (90/π^2 g_SM)^2∫_T_N^T_C dT M^4_Pl/T^5exp[- S_3/T] = 1.
  
  Since we expect that for strong first order phase transitions this integral is peaked around T_N, for practical purposes the approximate expression
  
      S_3/T|_T = T_N∼ - 4log[T_N/M_Pl]
   
  
  is often used to evaluate T_N. We note that we solve for T_N explicitly with Eq (<ref>) in our numerical implementation.
  
  Given the strong first order phase transition predicted by our model, one of its most interesting consequences and phenomenological signatures is the generation of gravitational waves (GW) from the well studied mechanism of bubble nucleation and collisions. The study of such a cosmological SGWB requires complex numerical simulations that are beyond the scope of this work. Instead, as in much of the literature about GWs from first order phase transitions, we present the most salient features of this phenomenon. The interested reader is invited to consult <cit.> and references therein for a broader discussion. As discussed in these references, the formulae we quote rely on simulations showing that a certain fraction of the latent heat released from the first order phase transition is deposited into graviational waves during the propagation of bubble wavefronts. For the regime we consider, GWs are primarily sourced by sound waves formed in the plasma after bubble collisions, but there are other secondary contributions and subtleties to consider: the bubble collisions themselves and in medium magnetohydrodynamic turbulence can become important, especially in regimes where the bubble walls are very relativistic. Other factors that can alter GW production is whether the background energy density is radiation or matter dominated and the lifetime of acoustic waves <cit.>. Contributions to produced GWs combine approximately linearly <cit.>:
  
      Ω_GWh^2≃Ω_Bh^2 + Ω_Sh^2 + Ω_Th^2 ,
  
  where Ω_S is the contribution from sound waves as the colliding bubbles dump their energy into the plasma, Ω_B the contribution from bubble collisions, and Ω_T the contribution from magnetohydrodynamic turbulence effects caused by the colliding bubbles. In this work we will focus mainly on Ω_S, since this is expected to produce the dominant contribution to gravitational waves for the model we consider. 
  
  In the treatment above, we have assumed that in a thermal background, the expansion of the universe is adiabatic, dT/dt = - TH, where H is the Hubble parameter. This gives us an important parameter that determines the properties of gravitational waves emitted from first order phase transitions,
  
      β/H|_H = H_N =  T  d/dT(S_3/T) |_T = T_N ,
  
  where H_N is the Hubble parameter at nucleation.
  In order to compute the resulting energy density of gravitational waves today, we must first compute the latent heat, or free energy density difference between the true and false vacuum at the time of the phase transition,
  
      ϵ = - V (ν_+,T_N) + T∂ V( ν_+,T)/∂ T|_T = T_N .
  
  Another important parameter for gravitational wave generation is the ratio of the latent heat released during the phase transition to the energy density in the high energy phase
  
      α = ϵ(T)/ρ_vac(T)|_T = T_N .
  
  Using α, the redshifted value of the gravitational wave energy density produced by acoustic waves has been shown
  <cit.> to follow a power law spectrum
  
      Ω_Sh^2 (f) = 2.6× 10^6 κ^2 v_b (α/1 + α)^2(100/g_*)^1/3(H/β) S_ac(f)
  
  where κ is the efficiency with which GW are generated by acoustic waves and S_ac(f) is a power law frequency spectrum obtained numerically,
  
      S_ac(f) = (f/f_ac)^3(7/4+3(f/f_ac)^2)^7/2.
  
  The peak frequency of gravitational wave production f_ac and average bubble wall velocity v_b will depend on a number of factors, including φ's coupling to the plasma, the equation of state of the universe, and the relative velocity of bubbles as compared to the plasma's sound speed <cit.>. Here we make some standard choices, that are supported by acoustic GW production simulations undertaken in <cit.>. We fix the average bubble velocity to v_b =0.9. For this velocity,
  
      κ≈α/0.73 +0.083 √(α)+α
  
  is the efficiency for transferring acoustic wave energy into gravitational energy <cit.>, and the peak frequency is f_ac≈2 β/√(3) v_b which redshifts to 
  
      f_ac=1.9× 10^-5Hz1v_b(βH)(T_N100GeV)(g_*100)^1/6 .
  
  
  
  
  
  
  With our model and gravitational wave observables defined, we are ready to investigate the generation of gravitational waves from first order phase transitions in the early universe. In Figures <ref> - <ref> we have plotted the critical temperature T_C, peak gravitational wave frequency f_ac, and gravitational wave energy density Ω_gw h^2, for our scalar potential, in terms of its three model parameters, the Yukawa coupling y, gauge coupling g, and scalar mass m_S. In these plots we have detailed phase transition signatures for a range of scalar masses m_S, from 10^6 GeV to 10^15 GeV. In producing these plots, we do not use the high temperature approximation, but rather solve numerially for the full thermal potential given in Eq. (<ref>). It is worth noting that higher order corrections and renormalization group improvements to the thermal potential, e.g. daisy corrections and taking into account the scale dependence of the effective potential do not have a significant impact on our results, as we discuss in Appendices (<ref>) and (<ref>). 
  
  
  
  To be specific, we compute the critical temperature T_C numerically by solving
  
      V(ν_+,T_C)=0 ,
  
  where the minimum ν_+ is defined by the value of the potential where dV/dφ=0 and d^2V/dφ^2>0. The nucleation temperature T_N is obtained by solving Eq. (<ref>), where S_3/T is taken from Eq. (<ref>). Note that at temperatures close to T_C, x in Eq. (<ref>) may exceed one, in which case f(x) in Eq. (<ref>) is not defined. However, this does impede our computations since at T∼ T_C, the energetic difference between the true and false vacuum is very small, and bubble nucleation is not yet efficient. We therefore set exp(-S_3/T)=0 for x>1 in solving Eq. (<ref>). The latent heat ϵ is attained numerically using Eq. (<ref>), which is then combined with Eq. (<ref>) to find the energy density sources as a gravitational wave background.
  
  In Figure <ref>, we show the final mass of the fermion field, m_χ, at the end of the phase transition. The parameter space where χ recieves a large mass from φ, which occurs at y ∼1, g/y ∼ 1, can be compared to the critical temperature, peak frequency, and SGWB densities for the same parameter space shown in Figures <ref>-<ref>. Generally a larger Yukawa coupling indicates a lower phase transition temperature, as can be understood from Eq. (<ref>) since both 1/D and T_0 decrease with a larger Yukawa coupling. As g/y approaches a particular value 3g^4 = 4y^4, the critical temperature grows drastically due to an exceptionally small effective scalar quartic, λ(T) ≪ 1. In Figure <ref>, we show the frequency at which the produced gravitational wave spectrum is peaked. As m_S varies from 10^15 GeV to 10^6 GeV, the peak frequency drops from more than a GHz to less than 100 Hz, which falls within the frequency reach of LIGO. Compared with T_C and the peak frequency f_ac, and the final fermion mass m_χ, the gravitational wave energy density today at peak frequencies is most sensitive to the gauge-Yukawa ratio g/y as can be seen from Figure <ref>. Intriguingly, several features, including a low gravitational wave peak frequency, high GW energy density, and a large fermion mass (which can be dark matter) are all attained when 3g^4∼ 4y^4. 
  
  
  
  
  
  
  § GRAVITATIONAL WAVES FROM INFLATION
  
  
  
  We now turn to a different potential source of gravitational waves: scalar field dynamics at the end of inflation. The inflaton models we will consider here attain large vacuum expectation values after inflation, much like the simplified thermal scalar potential studied in the previous section. Such potentials could therefore lead to mass boosted dark matter, as will be discussed in Section <ref>.
  
  If an inflaton field ϕ oscillates in its potential after inflation, it can self-interact in a way that causes large inhomogeneities to develop in its field density. These field configuration are sometimes called “oscillons” <cit.>. The formation of these inhomogeneities in the inflaton field density will source a gravitational wave background, which may be detected at future gravitational wave observatories. The possibility that inflatons could be detected through post-inflation gravitational wave production has been studied previously in <cit.>. This phenomenon is also interesting from the standpoint of detecting generic non-inflationary superheavy dark sectors, because this source of gravitational waves does not require a special thermal potential that results in a strong first order phase transition, like the one we investigated in Section <ref>. Compared to prior work on inflatons generating gravitational waves at the end of inflation, we will be considering both E- and T-model potentials, whereas previously only T-model potentials were studied in <cit.>. We will obtain a broad range of results for both E- and T-Model features, and find that E-model potentials result in different gravitational wave signatures than T-model potentials. 
  
  Here we will be interested in computing gravitational wave background spectra produced by two classes of inflationary potential, the E-Model
  
      V(ϕ) = Λ^4(1 - e^-√(2/3α)ϕ/M_Pl)^2
  
  and T-Model potentials
  
      V(ϕ) = Λ^4tanh^2(ϕ/√(6α)M_Pl),
  
  where Λ has dimensions of energy, ϕ is the inflaton field, and α is a constant that determines the shape of the potential. These models were proposed in <cit.>, where they were shown to arise from Kahler and superpotentials terms <cit.>. They are useful as generic inflation models, since for certain α values, E-Model and T-Model inflation are identical to a range of inflationary models which match present CMB observations, including Starobinsky <cit.> and Higgs inflation <cit.>. 
  Previously, gravitational wave background spectra from post-inflatonary field oscillations were computed for a T-Model inflationary potential in <cit.>. Here we extend this treatment to E-Model inflation, and compute results for a broader range of T-Model potentials.
  
  Another useful feature of the E-Model and T-Model classes of inflationary potentials, is that the free parameters of these potentials (Λ, α) are fixed by certain inflationary observables, including the primordial scalar power spectrum and the primordial tensor power spectrum. Specifically, the ratio of these primordial fluctuations, i.e. the tensor-to-scalar ratio (r) will have a fixed relationship to (Λ, α).  Here we define the tensor power spectrum and the scalar power spectrum as is customary,
  
      P_T(k) = A_T(k_*)(k/k_*)^n_T ,
  
  
      P_s(k) = A_s(k_*)(k/k_*)^n_s - 1
  
  where k is the wavenumber, k_* is the pivot tensor/scalar, A_T(k_*),A_s(k_*) are the amplitudes of the tensor and scalar power spectra which depend on the choice of k_*, and n_T,n_s are the tensor and scalar spectral indices. We have defined these using the usual conventions given in <cit.>; for a review on inflationary models see e.g. <cit.>. For slow-roll inflation, r, n_s, and A_s are expressible in terms of the slow-roll parameters ϵ and η <cit.>,
  
      ϵ = 1/2(M_PlV^'/V)^2,         η = M_Pl^2V^''/V,
  
  
      r = 16ϵ,          n_s = 1 - 6ϵ + 2η,          A_s = V/24π^2M_Pl^4ϵ,
  
  where M_Pl is the reduced Planck mass, V is the inflaton potential, and V^' and V^'' are the first and second derivative of the inflation with respect to ϕ. With these quantities defined, it can be shown that for E-Model and T-Model potentials, r, n_s, and A_s take the following form 
  
      r = 12α/N_e^2,      n_s = -2/N_e + 1,      A_s  = Λ^4N_e^2/ 18 π^2 α M_Pl^4 ,
  
  where these expressions are accurate so long as r ≲ 10^-3, and we have defined  as the number of e-folds of inflation lasting from time t_0 to t_e.
  
  In our study of gravitational waves from oscillating fields after inflation, we fix the parameters of our inflaton potentials using current results from Planck's 2018 data release <cit.>. Currently, Planck has observed n_s = 0.9649 ± 0.0042 at the 68% CL <cit.> and A_s = (2.100 ± 0.030) × 10^-9 at the 68% CL <cit.>. For simplicity, we fix the number of inflationary e-folds to N_e = 60. The r values explored in this paper are 10^-4, 10^-6, 10^-8, 10^-10 which are well below the current upper bound r < 0.061, set by combining Planck 2018 and BICEP2/Keck Array data <cit.>. Upcoming experiments such as LiteBIRD <cit.>, Simons Observatory <cit.>, and CMB-S4 <cit.> will be able to resolve an r value on the order of 10^-3. It is therefore interesting to note that this work explores r values much lower than those discoverable with any planned CMB experiments. The SGWBs we compute are given in terms of the fraction of the universe's energy density they take up,
  
      Ω_gw(f) ≡1/ρ_critdρ_gw/dlnf,
  
  where ρ_crit is the critical energy density and ρ_gw
  the energy density of gravitational waves.
  
  
  
  
  
  
  
  
  
  In order to determine the SGWB from E-Model and T-Model inflatons oscillating in their potentials at the end of inflation, we simulated the dynamics and self-interactions of these inflatons using the lattice field theory code HLattice <cit.>. Before the lattice simulation starts, HLattice solves for the inflaton field value at the end of inflation by evolving the field according to the Klein-Gordon equation, ϕ̈ + 3Hϕ̇ + V' = 0, until the following condition is met: ϕ̇ + Hϕ <0. Because we are studying low scale inflation models, we replaced this HLattice condition with ϵ≥ 1, where ϵ is the slow-roll parameter given in Eq. (<ref>). For low scale inflation models, this resulted in significantly smaller initial field values for ϕ, compared to the original HLattice condition.
  
  We ran simulations that lasted for just over an e-fold after the end of inflation, which corresponds to the scale factor of the universe, a, expanding by a factor of up to three for the T-Models. For the E-Model cases, the simulations lasted for just over an e-fold after inflation for r = 10^-4, with a expanding by a factor of three and less than an e-fold after the end of inflation, corresponding to the scale factor expanding by approximately a factor of a half for r =10^-6. The r = 10^-6 case was only run to this value of a, because it did not evolve significantly past this value and was therefore not shown. Given technical limits and finite computational resources, the r = 10^-8, 10^-10 cases were not computed for the E-Model. For all cases shown, we found that the gravitational wave spectra converged on a fixed result well before the end of the simulation, as can be verified in Figures <ref> and <ref>. Note that the spectra given in these figures for each a value, represent the energy density of gravitational waves at that a value redshifted to present day. The lattice box size (L) and box resolution (n) were chosen in order to ensure the simulation smoothly resolved SGWB spectral features. We checked that the results were robust to minor variations in the lattice spacing and simulation box size. For the E-Models, L = 0.003H^-1 and n = 256 was used for both r values. For the T-Models, L = 1.0H^-1 and n = 128 were used for all r values. The frequencies that the simulation probed ranged up to nπ/L with a spacing of π/L.
  
  
  
  
  
  
  
  
  
  Table <ref> gives simulation parameters used for the E-Model and T-Models respectively. Given a value of α, which determines the flatness of the inflationary potential, ϕ_0 was chosen such that inflation lasted for approximately 60 e-folds, i.e. N_e ≃ 60. The initial kinetic energy of the inflaton ϕ̇_0 was chosen to be very small, and is the same for every case. (For the sake of being completely reproducible, we note that ϕ̇_0 = -10^-10 in HLattice code, where HLattice uses the convention M_Pl=1024.) For comparison, in Table <ref> we have also given ϕ_e, the field value where the slow-roll conditions are no longer met, which we computed using standard slow-roll formulae. 
  
  
  In Figure <ref>, the stochastic gravitational wave background spectra are shown for E-Model inflationary potentials for tensor-to-scalar ratios r = 10^-4, 10^-6, where as usual, smaller r corresponds to a smaller energy density at the end of inflation. As a consequence of this decrease in energy density, the overall amplitude of the SGWB produced from inflation decreases with a decreasing r value. Contrary to prior work <cit.>, we do not observe growth in the SGWB consistent with the formation of oscillons for α≲ 10^-3 for the E-Models. We believe this is because our HLattice simulations include the effect of metric fluctuations and associated gravitational dynamics, which can affect the formation of oscillons, since as we observe, there is a large amount of background gravitational waves generated very soon after the end of inflation. This initial large inhomogeneity in the E-Model potential, as compared to the T-Model, may be the result of the E-Model potential not being symmetric around its minimum. We see in Figure <ref> that the SGWB amplitude decreases by approximately three orders of magnitude. Based on preliminary investigation, we do not expect oscillon configurations to form, and these SGWB to increase at later times.
  
  
  In Figure <ref>, the T-Model SGWB spectra were computed for r = 10^-4, 10^-6, 10^-8, 10^-10. Here the formation of oscillon field configurations is evident for the r = 10^-6, 10^-8, 10^-10 SGWB spectra, whose SGWB amplitudes increase markedly after a≈ 2, when oscillons form. This can be contrasted with r=10^-4, for which the SGWB amplitude is relatively static, because oscillons do not form. On the other hand, before the formation of these oscillons (at a=1), it can be seen that the initial amplitude of the stochastic gravitational background decreases with a decreasing r value, as does the frequency of peak SGWB amplitude. In comparison with prior results <cit.>, as expected we observe that T-Model inflatons only form oscillons for α≲ 10^-4. Here we find that, more precisely, oscillon configurations can still form for α≤3e-4. It is interesting to note that the peak amplitude of SGWB spectra decreases by a few orders of magnitude between r=10^-6 and r=10^-8, and by four orders of magnitude between  r=10^-8 and r=10^-10. 
  
  
  
  
  
  
  § HEAVY NONTHERMAL DARK MATTER FROM A MASS BOOST
  
  
  
  In this section we explore how scalar fields detailed in Sections <ref> and <ref>, can give a mass boost to nonthermally produced dark matter, resulting in dark matter produced out of thermal equilibrium, with a final dark matter mass of up to m_χ∼ 10^19 GeV. If dark matter is produced out-of-equilibrium at a temperature well above the temperature at which the universe reheats T_rh, then assuming that the universe remains radiation dominated from T_rh until the era of matter-radiation equality in the early universe, the relic abundance of dark matter can be computed by comparing it to the presently observed radiation energy density,
  
      Ω_DM h^2 = Ω_R h^2 ( T_rh/T_0) ρ_DM,e^i/ρ_e.
  
  Here ρ_e is the energy density of the universe at the end of inflation, ρ_DM,e is the dark matter energy density at the end of inflation, Ω_DM h^2 ≃ 0.12 is the physical dark matter density parameter,  Ω_R h^2 ≃ 4.3 × 10^-5 is the physical radiation density parameter, and T_0 ≃ 2.34 × 10^-13 GeV is the present temperature of the universe. In an expanding universe, radiation energy density dilutes as T^4, while matter dilutes as T^3, so the factor of T_rh/T_0 accounts for shifts in the density of matter relative to radiation as the universe expands from T_rh to T_0. Note that there is an implicit assumption in Eq. (<ref>) that the ratio ρ_DM,e/ρ_e does not shift significantly between the end of inflation and the time at which the inflaton decays – in practice this may introduce a small correction to the final dark matter relic abundance.
  
  The standard method to compute the abundance of gravitationally created particles at the end of inflation ρ_DM, e, is to integrate the particle's Bogoliubov operator to the end of inflation <cit.>. The precise abundance of particles created will depend on cosmological dynamics following inflation <cit.>, however, to good approximation the gravitationally produced relic abundance for a fermion with mass m_fi is
  
      ρ_DM,e≃ 10^-2 H_e^4 e^-2m_fi/H_e,
  
  where this expression holds so long as m_fi≳ H_e, where H_e is the Hubble constant at the end of inflation.
  
  Here we will be most interested in how the relic abundance is altered, if after a density ρ_DM,e of dark matter particles is created at the end of inflation, these particles undergo a mass shift m_fi→ m_χ during a dark sector phase transition. This mass shift would boost the dark matter mass, and corresponding relic abundance, by a factor m_χ/m_fi. For simplicity, we can assume that the universe remains radiation-dominated after reheating, so that there are not additional corrections to the dark matter relic abundance, from e.g. a post-reheating period of matter-dominated expansion. Then substituting the Friedmann relation ρ_e =3 M_Pl^2 H_e^2 and Eq. (<ref>) into (<ref>), we arrive at the total relic abundance for dark matter that is gravitationally produced,
  
      Ω_DM h^2 ≃ 0.003 Ω_R h^2 ( T_rh/T_0) H_e^2/M_Pl^2 e^-2m_fi/H_e( m_χ/m_fi).
  
  The final term in this expression accounts for dark matter undergoing a mass boost from m_fi→ m_χ, at some time after the end of inflation. Note that this expression for the relic abundance is only correct so long as m_fi≳ H_e, which is required by the cosmology detailed above. Finally, we note also that Eq. (<ref>) does not assume instantaneous reheating, allowing for a Hubble constant at reheating that is much smaller than the Hubble constant at the end of inflation.  This would occur, for example, in models of inflation where the decay of the inflaton is slow compared to the Hubble rate at the end of inflation.  In particular, in Eq. (<ref>) T_rh has been fixed to achieve the correct relic abundance of dark matter.
  
  In Figure <ref> we show parameter space for mass-boosted dark matter, in terms of the Hubble constant at the end of inflation H_e, the temperature of the universe at reheating T_rh, the initial dark matter mass m_fi, and the final dark matter mass m_χ. The present dark matter relic density has been fixed to match observations for a number of initial dark matter masses m_fi. For simplicity, we have restricted the parameter space to the regime of validity of Eq. (<ref>), so that the initial mass of the dark matter m_fi > H_e exceeds the Hubble constant at the end of inflation. In Figure <ref> it must also be required for consistency that the reheat temperature is low enough relative to the energy density associated with H_e. Thus, we have truncated H_e y-axis value so that H_e ≫ H_rh throughout this plot for consistency.o that H_e ≫ H_rh throughout this plot for consistency.
  
  
  
  
  We now provide a discussion of the relic abundance of dark matter, as pertains to Eq. (<ref>). There are a number of conditions which can be specified to ensure that the relic abundance of dark matter given by Eq. (<ref>) represents the total relic abundance of dark matter, particularly so that additional dark matter is not produced by other processes, for example dark matter coming into thermal equilibrium with the radiation bath after the universe reheats at T_rh. These conditions will depend on the particulars of the mass-boosting cosmology after the end of inflation. 
  
  First we discuss when a mass boost will occur for the simple scalar potentials we will consider. In order to provide a compact treatment of mass-boosting scalars in this section, we will be considering simplified versions of the scalar potentials detailed in Section <ref> and <ref>. In particular, we will focus on the fact that these scalar potentials V ( φ) are sufficiently flat in the early universe, that V”( φ_i) ≪ H^2. We assume that these scalars attained initial field values close to φ = 0 prior to inflation. During inflation, such a scalar will obtain a vacuum expectation value from Hubble fluctuations of order φ≈ H_e. As previously stated, for all cases we consider H_e ≪ m_fi, so this Hubble-induced vacuum expectation value φ∼ H_e will not alter the mass of χ during inflation, since from inflationary fluctuations the induced mass will be y H_e≲ H_e. In this scenario, the scalar acts just like a curvaton field <cit.>, or like any scalar/
  axion/modulus with a flat potential in the early universe, and the initial field value φ∼ H_e remains fixed by Hubble friction so long as V”( φ_i) ≲ H^2, see e.g. <cit.>. Therefore we can safely assume that during and after inflation φ maintains its initial field value because of Hubble friction, until the Hubble constant drops to a smaller value than the scalar mass, which is when H^2 ≲ V”(φ) ≡ m_eff^2. This provides the first restriction on the dark matter mass-boosting scalar considered here: the effective scalar mass should be smaller than the Hubble constant at the end of inflation V”(φ) ≡ m_eff^2 < H_e^2. With this simple scalar framework laid out, we are now ready to address the additional conditions required to ensure Eq. (<ref>) is the majority contribution to the abundance of dark matter. This will depend on whether the mass boost occurs before or after reheating.
  
  
      
    (A) m_χ > T_rh.  First we consider the case that a dark matter mass boost occurs after the end of inflation, but before reheating. For the scalar fields we consider in this section, this can be enforced by the requirement that the effective mass of the scalar field m_eff follows H_e > m_eff > H_rh, so that the scalar will roll down its potential between the end of inflation and T_rh. In this case, it is safe to assume that dark matter will not be produced in the thermal bath after reheating if m_χ≫ T_rh .
      
      
    (B) m_fi > T_rh. Next we consider that case that a mass boost occurs only after reheating, but the initial dark matter mass m_fi exceeds the reheating temperature. In this case again, it is usually safe to assume a negligible abundance χ is produced after reheating. In this case, the parameter space below the m_fi = 10^9 and 10^11 GeV lines in the left and right panels of Figure <ref> should be excluded.
      
      
    (C) m_fi < T_rh. Lastly for completeness and because it has some bearing on (A) and (B) we consider the possibility that a mass boost occurs after reheating, with the dark matter fermion mass smaller than the reheating temperature. This can be accommodated so long as χ and ϕ are weakly coupled to the radiation bath (because χ is gravitationally produced, no non-gravitational coupling need be assumed). In this case the requirement that χ is not substantially produced depends on the coupling to Standard Model particles in the radiation bath, for a complete treatment of “freeze-in” production and how small a coupling is required, we defer to <cit.>. In the following, for simplicity we will only consider cases (A) and (B).
  
  
  
  So far our treatment has been general, since we have not specified what scalar field provides the mass boost indicated by the ratio m_χ/m_fi. We now consider two possibilities: this scalar field may either be a scalar field that goes through a phase transition, much like that detailed in Section <ref>, but here different because we will set the gauge coupling to zero for simplicity, or a scalar field that rolls to a large field value after inflation, as explored in Section <ref>. 
  
  First we consider a mass-boosting scalar field, with a quartic self-coupling generated at loop-level from its Yukawa coupling to a fermion as laid out in Section <ref>. Setting the gauge coupling to zero (g=0), and omitting the logarithmic term which will yield a small correction when φ has settled at its minimum, Equation (<ref>) becomes
  
      V ( φ ) = -m_S^2/2φ^2 + 17y^4/48 π^2φ^4,
  
  where we remind ourselves that y is also the Yukawa coupling between the scalar and the dark matter field χ. 
  
   After the mass-boosting scalar rolls to its minimum, the dark matter mass will be boosted through its Yukawa coupling to the scalar. In the case that the Yukawa mass is much larger than the dark matter's initial mass, the final mass of the dark matter will be m_χ≈ y ν, where ν is the final vacuum expectation value of φ. Then assuming that after φ has settled at its minimum, the Yukawa term i y χ̅χφ provides the majority of the dark matter mass, the dark matter mass will be
  
      m_χ≈ y ν = √(24/31)π m_S/y,
  
  where here we have used the exact vacuum expectation value as given in Eq. (<ref>). 
  
  It will be instructive to consider some benchmark parameter using Eqs. (<ref>) – (<ref>). 
  
  
    Benchmark I: For m_fi = 10^14 GeV, m_S = 10^12 GeV, H_e = 10^13 GeV, and T_rh=10^11 GeV (cf. the right panel of Figure <ref>) the correct relic abundance is obtained for m_χ = 10^15 GeV, corresponding to y = 2.8 × 10^-2 and ν = 4 × 10^17 GeV. For this benchmark point, φ will roll to the bottom of its potential and boost the mass of χ when H ∼ m_S∼ 5 × 10^12 GeV, well before reheating at H_rh = 2300 GeV.
  
  
    Benchmark II: For m_fi = 10^11 GeV, m_S = 10^9 GeV, H_e = 10^13 GeV, and T_rh=10^9 GeV (cf. the left panel of Figure <ref>) the correct relic abundance is obtained for m_χ = 10^13 GeV, corresponding to y = 2.8 × 10^-3 and ν = 4 × 10^16 GeV. For Benchmark II, φ will roll to the bottom of its potential and boost the mass of χ when H ∼ m_S∼ 10^9 GeV, well before reheating at H_rh = 0.2 GeV. In this case since m_S ∼ T_rh, φ could undergo a first order phase transition as in Section <ref>.
  
  
  
  In principle, arbitrarily large dark matter masses can be obtained in the limit that y → 0 in Eq. (<ref>), because for the model detailed in Section <ref>, the Yukawa coupling scales inversely with the scalar's vacuum expectation value, y ∝√(m_S / ν), a large-enough dark matter mass implies a super-Planckian vacuum expectation value for φ (where here we reiterate that the gauge coupling g=0 for simplicity), 
  
      m_χ≈ 4 × 10^16  GeV( m_S/6 × 10^13  GeV) ( ν/10^19  GeV)^1/2
  
  where here we normalized the scalar mass (which must satisfy m_S ≲ H_e) to the maximum value allowed by the Planck 2018 bound <cit.> on the Hubble constant at the end of inflation H_e^ max≲ 6 × 10^13  GeV. We note that this does not necessarily restrict the maximum mass of mass-boosted dark matter, since there are a number of models that accommodate super-Planckian scalar field vacuum expectation values <cit.>.  However, there is one additional considerations in this regard. Namely, we should consider whether a vacuum energy added to cancel the contribution of V(φ) at its minimum, which would be ∼ - m_S^2 ν^2 in Eq. (<ref>), implies some departure from the cosmology we have already laid out. Because there exist mechanisms that dynamically adjust the cosmological constant (see e.g. <cit.>), such an addition is not necessarily required. However, if we assume this additional term is present, this would place an additional restriction that ν≲ M_Pl, since this would ensure that when φ begins rolling to the bottom of its potential after inflation at H ≈ m_S, the vacuum energy term would be subdominant to the universe's energy density at this time, i.e. 3 H^2 M_Pl^2.
  
  Lastly, we discuss under what conditions the inflationary potentials examined in Section <ref> can provide a mass boost to gravitationally produced dark matter.  These potentials have fields that move over Planckian field ranges (see Table <ref>), and so it is natural to consider whether they might provide a large mass boost to dark matter, which is gravitationally produced towards the end of inflation. Given that the inflaton field value itself is shifting appreciably during inflation, so too will the dark matter's mass, considerably complicating the resulting dynamical production of dark matter. The standard Bogoliubov coefficient evolution that yields Eq. (<ref>), would need to be re-computed to incorporate a dynamical mass for χ. We leave a full investigation of this to future work. 
  
  Here we will consider the simpler case, that the scalar potentials given in Section <ref> are not inflatons, but rather low-mass scalar fields, that stay fixed at some field value during and after inflation ϕ = ϕ_i, and then roll to their minima sometime after the end of inflation. For the case we consider, we will choose to set this initial field value to the “ϕ_e” field values given in Table <ref>. We choose these field as a starting field value, because for larger field values, these potentials would initiate a second period of inflation. (This is allowed <cit.>, but beyond the scope of our consideration here). In this case, we note that the scalar would couple to the fermion χ with a Yukawa coupling ℒ_yuk = - y χχ̅(ϕ - ϕ_e). The dynamics of these potentials are similar to curvaton fields <cit.>, in that they will slowly roll to their final field value sometime after the end of inflation, although in this case we do not require that they necessarily produce the observed spectrum of primordial perturbations, as for a curvaton field. 
  
  The dynamics of dark matter production in this scenario will be the same as above, but with the additional restriction that the effective mass of the late-time rolling scalars is less than the Hubble constant at the end of inflation
  
      m_eff^2 = V”(ϕ_e) ≪ H_e^2
  
  This condition ensures that the scalar will not begin oscillations until after the end of inflation, and this can be satisfied for all of the potentials studied in Section <ref>. Finally, note that so long as the scalar field is free to oscillate in its potential (and is not for example damped by scattering with background fields) it should produce the same gravitational wave signatures detailed in Section <ref>. For concreteness, let us briefly examine the E-Model for α = 3 × 10^-4. In this case, the effective mass m_eff = √( V”(ϕ_e))≈ 2 × 10^12  GeV, well below the bound H_e ≲ 10^14  GeV indicated in Figure <ref>. Looking at Table <ref>, the vev of this scalar, and the corresponding mass of a coupled dark matter field (m_χ≈ y ϕ_e) will shift by up to ∼ 0.1 M_Pl at the end of the scalar's oscillations.
  
  
  
  § CONCLUSIONS
  
  
  
  
  
  
  Gravitational waves produced by scalar fields may soon reveal details from the first moments of our universe. We have found connections between gravitational waves sourced by scalars and dark matter, including a mechanism whereby scalar fields boost the mass of non-thermally produced dark matter. For many cases of mass boosted dark matter, there will be an associated SGWB, produced through a first order phase transition or oscillations of a scalar field after inflation. In our study of first order phase transitions, we particularly examined a scalar field whose self-interactions are generated at one-loop by its Yukawa coupling to a fermion and gauge coupling to a vector field. For oscillating fields after inflation, our lattice simulations have extended prior work on oscillating scalar field configurations formed by E-Model and T-Model inflationary potentials, and their SGWBs. Both of these scalar fields tend to produce larger amplitude gravitational wave backgrounds for larger dark matter masses.
  
  
  
  
  Figure <ref> shows stochastic gravitational wave spectra from inflation, phase transitions, and sensitivity curves for upcoming and current gravitational wave experiments. The E-Model and T-Model spectra shown correspond to the results from Section <ref>, with T-Model results ranging from r=10^-4 - 10^-10 and E-Model results ranging from r=10^-4-10^-6.  Gravitational wave spectra from phase transitions are shown for g = 1.05, y = 1, and m_S = 10^3- 10^15 GeV as indicated. Sensitivity curves are also shown, both for direct GW observatories, and indirect cosmological probes. Dashed-dotted projected sensitivity curves include the contribution of astrophysical foregrounds while dotted curves do not include astrophysical foregrounds. Solid curves represent bounds on SGWB. The sensitivity curves for BBO, LISA, and aLIGO were taken from <cit.> and take into account the number of detectors, mission observation time, and the overlap reduction function of each detector pair. The BBO and LISA curves additionally take into account astrophysical foregrounds. The LIGO sensitivity curve for the first observing run (O1), the first and second observing run (O1+O2), and the design sensitivity are taken from the data supplement to <cit.>. A bound on SGWB contributions to relativistic degrees of freedom, derived from a fit to CMB power, baryon acoustic oscillations, lensing, and primordial nuclear abundances and prospects for an improvement of this cosmological SGWB bound, using projected constraints from the COrE and EUCLID satellites is taken from <cit.>. 
  
  There are a number of interesting conclusions that can be drawn from Figure <ref>. For low mass scalar phase transitions, the peak of the SGWB could possibly be detected in upcoming space-based gravitational wave observatories such as the Big Bang Observatory (BBO). If a SGWB is detected, its provenance from either acoustic waves during a first order phase transition or oscillating field dynamics could be determined using the strikingly different SGWB spectral features apparent in these two cases. While BBO appears sufficient to find dark sector phase transitions for scalar masses m_S ∼ TeV, our work has reinforced the importance of MHz-GHz searches, since the peak amplitude of many SGWB signatures lie in the MHz-GHz frequency range, motivating the development of sensitivity in this regime. Experiments have so far obtained sensitivity to MHz-range SGWBs exceeding the current critical density, e.g. Ω_gw,0 h^2 ≲ 10^13 at MHz frequency <cit.> and Ω_gw,0 h^2 ≲ 10^26 at 0.1 GHz <cit.>. There are some laboratory proposals at GHz frequencies <cit.>. However, even in the absence of laboratory tests, the cosmological impact of MHz-GHz gravitational waves, which can contribute to the number of relativistic degrees of freedom in the early universe <cit.>, or convert into photons in the presence of a magnetic field, are both promising high-frequency SGWB detection methods <cit.>. 
  
  In conclusion, we have shown that a range of gravitational waves background signatures could reveal the early universe dynamics of an inflaton or hidden sector phase transition, and might also be a signature of very heavy dark matter which received a mass boost in the early universe. In future work, it would be interesting to examine the extent to which gravitational spectra are produced, as the scale of inflaton potentials and corresponding tensor-to-scalar ratio r are decreased by many orders of magnitude. It would also be interesting to revisit our computations at higher loop order, to allow for analysis of the simplified potential we have studied in a more strongly-coupled limit. In addition, while we have studied a very simple model demonstrating links between mass boosts for heavy dark matter, gravitational waves from phase transitions, and a simplified scalar potential, these connections could be re-examined for an enlarged model, particularly since we have found observable phase transitions that could be associated with Grand Unified Theories predicting many scalars with masses around 10^15 GeV.
  
  
  
  
  
  § ACKNOWLEDGEMENTS
  
  The work of AB, JB, SN, NS is supported by the Natural
  Sciences and Engineering Research Council of Canada (NSERC). We thank, Avi Friedlander, Alan Goodman, Andrew Long, Amalia Madden, and Pauline Perrin for helpful discussions and especially Mustafa Amin, Djuna Croon, Kaloian Lozanov, and Graham White for comments on the manuscript. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities. 
  
  
  
  
  
  
  § THERMAL CORRECTIONS AND HIGHER ORDER EFFECTS
   
  
  This appendix provides further details of calculations for the thermal potential discussed in Sections <ref> and <ref>.  The validity of the perturbative expansion is addressed using analytical arguments and numerical results which show that daisy contributions, which can dramatically affect the nature of phase transitions, are subdominant. The derivations are similar to those of <cit.>, with a few key differences. 
  
  We begin with the calculation of the finite temperature effective potential. It will be useful to go into some detail about the more formal aspects of this exercise to highlight how our model differs from scalar potentials with a tree level quartic term. As with the zero temperature case, we use the background field method. To compute the one loop effective potential at a finite temperature T, we will have to calculate an integral of the form <cit.>:
  
      I(φ) = T/2∑_j c_j∑_n∫d^3𝐤/(2π)^3 log[k_E^2 + M^2_j(φ)/T^2]
  
  Here the sum over j is a sum over particle species and the expression k_E^2 + M^2_j(φ) is the inverse of the square of the propagator for that particle species (except for scalars where it is just the inverse of the propagator). We will return to this shortly, but for now let us proceed with evaluating this integral. k_E is the Euclidean momentum at finite temperature, k_E^2 = ω_n^2 + k^2, with ω_n = 2π T n,  2π T (n+1) the Matsubara frequencies for bosons and fermions respectively and k the three momentum amplitude. To evaluate (<ref>), it is useful to define E^2≡ k^2 + M^2(φ) and use the following: 
  
  
      v_B(E)     = ∑_n = -∞^+∞log[4π^2T^2n^2 + E^2], ∂ v_B/∂ E = ∑_n = -∞^+∞2E/(4π^2T^2n^2 + E^2)
  
  for bosons and
  
      v_F(E)     = ∑_n = -∞^+∞log[4π^2T^2(n+1)^2 + E^2],  ∂ v_F/∂ E = ∑_n = -∞^+∞2E/(4π^2T^2(n+1)^2 + E^2)
  
  for fermions. The sums in ∂ v_B,F/∂ E can be evaluated exactly using 
  
      ∑_n = 1^+∞y/y^2 + n^2 = -1/2y + π/2(π y)
  
  With (<ref>), ∂ v_B,F/∂ E can be written in closed form and integrated with respect to E to obtain v_B,F. The first term of the series sum (<ref>) gives the zero temperature effective potential while the second one gives the thermal part. We have skipped a step by step derivation here but for further details the reader may consult in particular references <cit.> and <cit.>. The integral expressions in (<ref>) are easily obtained, bearing in mind that gauge boson contributions come with a factor of c_B = 3 and fermions a factor of c_F = -4 from the Dirac trace. Let us however maintain a more general presentation. The one loop effective potential is then written as:
  
      V_1(φ,T) = V_0(φ,0) + V_1(φ,0) + T^4/2π^2[c_BJ_B(φ,T) + c_FJ_F(φ,T)]
  
  where V_0(φ,0) and V_1(φ,0) are the tree level and one-loop effective potentials at zero temperature respectively. The thermal part of the potential reads 
  
      J_F,B(φ,T) = ∫_0^∞ dx  x^2log[1 ±exp(√(x^2 + M^2_F,B(φ)/T^2))]
  
  The high temperature expansion is obtained by using the well know expressions, equations (174) and (200) of <cit.>, for the bosonic and fermionic contributions respectively. 
  
  
  
   §.§ Daisy Resummation
  
  Let us now pay closer attention to 
  the argument inside the logarithm for expression (<ref>). The M^2 term is just the tree level field-dependent mass of the boson/fermion that contributes to the effective potential. We know, however, that this mass changes due to radiative corrections, and if we could include the effect of this “effective mass" into our potential we would have a more accurate result. This is not just a matter of accuracy. At zero temperature, we only require that the coupling constants be small to ensure the validity of perturbation theory but at finite temperature this is not sufficient. Loop corrections at finite temperature depend both on coupling constants and temperature, and a high temperature can compensate for a small coupling constant, making perturbation theory unreliable - higher order corrections to the one-loop effective potential become increasingly important. 
  
  To assess the impact of loop corrections, we need to define a loop expansion parameter l which is dimensionless and depends on the coupling constants, temperature, and cutoff scale of the theory. The contribution to the effective potential of a loop amplitude of superficial degree of divergence d is of the order l^d. If l^d≪ 1 higher order corrections, which come in powers of l^d, are negligible. At fixed order in perturbation theory daisies, loops with d=2, provide the largest contribution. Therefore, if l^2 is 𝒪(1) we cannot reliably trust perturbative methods. Fortunately, there is a relatively straightforward procedure to include their contributions to all orders in perturbation theory. Daisy resummation consists of calculating M^2 to some order in perturbation theory and substituting M^2→ M^2_eff into the original expression for the effective potential. For reasons that will soon be apparent, let us - as is conventionally done in <cit.> for example - write:
  
      M^2 = M_0^2 + Π^2(0)
  
  where Π^2 is called the polarization tensor of the particle species in question, calculated to some order in perturbation theory. With this modification, the one-loop effective potential in (<ref>) is modified to
  
      I(φ) = T/2∑_j c_j∑_n∫d^3𝐤/(2π)^3 log[k_E^2 + M^2_j + Π_j^2(0)/T^2]
  
  It is useful to seperate out this expression so that we may write the final expression for the effective potential as: 
  
      V_1(φ,T) = V_0(φ,0) + V_1(φ,0) + T^4/2π^2[c_BJ_B(φ,T) + c_FJ_F(φ,T)] + V_D(φ,T).
  
  This potential is identical to (<ref>) except for the daisy “correction" V_D(φ,T),
  
      V_D(φ,T) = -T/2∑_j c_j∑_n∫d^3𝐤/(2π)^3 log[1/T^2(1 + Π_j^2(0)/k_E^2 + M^2_j)]δ^0_n,
  
  where we have used the Kronecker delta δ^0_n to illustrate that in the Daisy resummation only the infrared (IR) divergent n=0 energy mode is important. This is equivalent to a high temperature expansion <cit.>. The interested reader is invited to consult <cit.> for an alternative but nonetheless in depth derivation of this effective potential. 
  
  Proceeding, the resulting integral is
  
      V_D(φ,T) = -T/2∑_j c_j∫d^3𝐤/(2π)^3 log[1/T^2(k^2 + M_j^2(φ) + Π_j^2(0))] - log[1/T^2(k^2 + M_j^2(φ))],
  
  where k is the three momentum integral. The integral can be evaluated using dimensional regularization,
  
      ∫d^Dk/(2π)^Dlog(k^2 + Δ^2) = Γ(-D/2)/(4π)^D/2(Δ^2)^D/2,
  
  where D = 3 + ϵ and ϵ is set to zero at the end. The final result is
  
      V_D(φ,T) = -T/12π∑_j c_j(M_j^2(φ) + Π_j^2(0))^3/2 - (M_j^2)^3/2
  
  
  This is the general form of the daisy potential. For our purposes, all we need now to calculate its exact form are polarization tensors for the scalar and vector, Π^2_S(0) and Π^2_V(0) respectively. We only consider bosonic contributions as fermionic ones are not IR divergent. Let us first consider the scalar contribution as it is the most important difference between our model and prior calculations in the literature. In the interest of generality, let us assume for the moment that we have a tree level term λ/4φ^4 in our lagrangian. We will set λ to zero at the end of the calculation. The expression for M^2_S to one loop order is,
  
      M^2_S = -m_S^2 + 3λφ^2_c + Π^2_S(0)
  
  The tree level contribution is M^2_0 = -m_S^2 + 3λφ^2_c, which has a field dependent value due to the quartic term. 
  
  Restricting ourselves only to first order contributions, we find that the following diagrams contribute at first order to the scalar polarization tensor.
  
  The contribution due to the self interaction λ 
  
      Π^2_SS(k) = 3λ T ∑^∞_n=-∞∫d^3k/(2π)^31/k^2+ m_S^2 + (2π n)^2 T^2|_T≫ m_S
  
  where the SS subscript in Π^2 indicates that we are calculating the scalar contribution to the scalar polarization tensor. We are typically interested in the high temperature behaviour of the system, so it is customary to set m_S = 0 since its contribution is negligible at that temperature. Following the same steps as before, i.e., performing the sum first and then the three momentum integral, we find
  
      Π^2_SS(0) =Π^2_SS(k) |_k=0 = λ T^2/4.
  
  The other polarization tensors are calculated in the same manner. They are <cit.>
  
      Π^2_SF(0)     = y^2/2T^2
      Π^2_SV(0)     = g^2/2T^2
  
  for the fermionic and vector contributions to the scalar polarization tensor respectively. This is everything we need to write the final form of the scalar contribution to the daisy potential 
  
      V_DS(φ,T) = -T/12π[ (-m_S^2 + 3λφ^2_c + λ T^2/4 + y^2/2T^2 + g^2/2T^2)^3/2 - (-m_S^2 + λ/2φ^2_c)^3/2].
  
  The difference between our model and effective potentials with a tree level quartic interaction is now much clearer: in our model, the scalar does not contribute to daisy contributions. This is evident by setting λ to zero - the φ dependence of V_DS(φ,T) vanishes and V_DS only contributes an irrelevant constant to the daisy resummed effective potential.  
  
  
  For our model then, the only contribution to the daisy potential comes from the gauge bosons. Here the computation is slightly tricky since only the longitudinal part of the gauge boson field A_μ contributes. Let us start by defining transverse, T_μν, and longitudinal, L_μν, projection tensors
  
      T_μν    = η_μν - k_μk_ν/k^2
      
              L_μν    = k_μk_ν/k^2,
  
  where η_μν is the Minkowski metric. It is easily verified that T_μνT^μρ = T_ν^ρ, L_μνL^μρ = L_ν^ρ, and T_μνL^μν = L_μνT^μν = 0. Then the daisy potential can be written as
  
      -T/2∑ c_j∫d^3𝐤/(2π)^3 log[1/T^2(k^2 + (gφ)^2 + (T_μν + L_μν)(Π^2)^μν(0))] - log[1/T^2(k^2 + (gφ)^2)]
  
  To proceed, we will simply use a result obtained in <cit.>, which the interested reader is invited to consult for a more in depth presentation:
  
      Π^2A_μν(0)     = - L_μνΠ^2A_00(0) 
      Π^2A_00(0)     = (gT)^2/3
  
  Here we see the purpose of introducing the polarization tensors: the key contribution to the daisy potential from gauge bosons come from their longitudinal polarizations. Putting everything together, the daisy potential for our model is therefore:
  
      V_D(φ,T) = - g^3T/12π[ (φ^2 + T^2/3)^3/2 - φ^3]
  
  Having computed the daisy resummed one loop effective potential, we now wish to assess the extent to which it affects our original calculations and preserves the nature of our phase transition. This is of course a difficult question requiring complex lattice simulation which are beyond the scope of this work. Here we will restrict ourselves to computing the critical temperature. If daisy corrections are significant they can change the nature of the phase transition from first to second order. The potential loses its valley - hill - valley shape at high temperature and there is no temperature at which the effective potential evaluated at its non-zero minimum is equal to zero, i.e., there is no solution for T_C. Therefore, determining regions of parameter space where T_C exists for the daisy improved effective potential is a good indicator of a first order phase transition. We consider a 10^15 GeV mass scalar and, as Fig. (<ref>) shows, the critical temperature does not change drastically and daisy contributions are therefore subdominant.  
  
  
  
  
  
  
  
  Retrospectively, let us examine why we should have expected this outcome from power counting arguments. Since Yukawa interactions do not contribute to the potential barrier between the two minima to leading order in perturbation theory, we neglect them. At high temperature the (non daisy resummed) effective potential when the temperature is close to T_0 can be written as:
  
      V(φ, T) = ν_eff^2φ^2 - E T φ^3 + λ(T)/4φ^4,
  
  where we have written the quadratic term in (<ref>) as ν^2_effφ^2, where ν_eff is temperature dependent. Near T_0 all three terms are roughly the same order of magnitude. This potential has a minimum φ∼g^3/λ_T_0 and an effective mass ν^2_eff∼g^6/λ_T_0T^2 obtained by equating the second and third terms and first and third terms respectively. The loop expansion parameter is given by l ≃g^2T/gφ∼λ_T_0/g^2 <cit.>. In our case, it can be verified that for the parameter space we are interested in λ_T_0≪ g^2 so higher order loop corrections to the effective potential are negligible.
  
  
  
  § RG IMPROVEMENT OF THE EFFECTIVE POTENTIAL
   
  
  Here we discuss how the effective potential depends on the scale μ at which it is evaluated, i.e. the renormalization group (RG) improvement of the model. We will find that for our model, all couplings run very weakly with this scale. In our setup for the zero temperature one-loop effective potential, we define the effective potential at a fixed scale μ = ν, where ν is given by (<ref>). However, it is known that both the field value φ and couplings g(μ) and y(μ) depend on μ and the values we used in our calculations are really just the field and coupling constant values at scale ν. Our effective potential at another scale μ≠ν will be different from the one we have used for our calculations in this work and these are to be trusted only if that difference is small. In this Appendix we show that this is indeed the case.
  
  
  
  
  Let us start with the field value φ. A shift in scale μ→μ^' shifts the field value according to <cit.>
  
      φ   →φξ(t) 
      ξ(t)     = exp[-∫_0^tdt^'γ(t^') ]
  
  where t(μ) = log(μ/ν) has been introduced for convenience. γ(t^') is the anomalous dimension of the scalar field. To one loop order it is given by
  
      16π^2γ(t) = 6y^2(t) - 3g^2(t).
  
  
  
  Equations (<ref>) and (<ref>) show that as long as the evolution of the coupling constants is small for the μ range under consideration, the effective potential does not change significantly. The couplings evolve according to their RG equations
  
      16π^2d g/dt    = g^3/3,             g(t) = g_ν/1 - g_ν^2/24π^2t
      
              16π^2d y/dt    = 5y^3,             y(t) = y_ν/1 - 5y_ν^2/8π^2t,
  
  where we again remind ourselves that g_ν,y_ν are the values of the gauge and Yukawa couplings at μ = ν respectively. Since t = log(μ/ν), there is a logarithmic dependence on scale for both couplings.
  
  In our case we started by defining the coupling values at a scale μ = ν and then considered thermal effects leading to a first order phase transition, for which the scale range is naturally the range of temperatures T_0≲ T ≲ T_C. To show that our calculations can be trusted, it is sufficient to prove that two conditions are met: (1) that the couplings do not vary significantly in the range T_0≲ T ≲ T_C and (2) that the coupling values in that range are not too different from the values at μ = ν. 
  
  We find that typically T_C ≃𝒪(ν), with ν higher by 𝒪(1) factors. Only in the special cases mentioned previously where y ≃ g/1.05, does ν differ from T_C by a few orders of magnitude. Therefore, both conditions mentioned above are met if we can show that the couplings do not vary significantly in the range T_0≤μ≤ν, which we plot in Figure <ref>.
  As we can see clearly, both the gauge and Yukawa couplings meet the two conditions established above. This can also be understood from the fact that the couplings depend logarithmically on scale, their variation over the range of interest is small even if the latter spans multiple orders of magnitude - a fractional amount for a scale range spanning a few orders of magnitude, for example. Therefore, our calculations appear to be trustworthy since the effective potential does not vary significantly over the range of interest to us. 
  
  
  
  JHEP.bst
  
  
  
  