Gravitational-wave signatures from reheating

We initiate a study of the gravitational-wave signatures of a phase transition that occurs as the Universe's temperature increases during reheating. The gravitational-wave signatures of such a heating phase transition are different from those of a cooling phase transition, and their detection could allow us to probe reheating. In the lucky case that the gravitational-wave signatures from both the heating and cooling phase transitions were to be observed, information about reheating could in principle be obtained utilizing the correlations between the two transitions. Frictional effects, leading to a constant bubble-wall speed in one case, will instead behave as an ``antifriction'' force in the other and accelerate the bubble wall. This antifriction will often take the bubble into a runaway regime, significantly enhancing the amplitude of the heating phase transition gravitational-wave signal. The efficiency, strength, and duration of the phase transitions will be similarly correlated in a reheating-dependent way.


I. INTRODUCTION
The remarkable transparency of the Universe to light allows us to look far back in time and learn about the early Universe.Using this fact, we can observe the clumping of matter as a function of redshift, as well as infer early Universe properties from the cosmic microwave background (CMB) power spectrum.At around a redshift of z ∼ 1100, however, this treasure trove of information ceases as the Universe becomes opaque to light.
Gravitational waves (GWs) offer a unique opportunity to look even further back in time.Unlike light, the Universe is never opaque to GWs, and thus they allow us to observe the Universe at its very youngest.This unique opportunity comes at the cost of it being much harder to observe GWs than it is to observe light.It is only recently that they have been observed for the first time by the LIGO-Virgo Collaboration [1].With future detectors such as LISA [2][3][4][5], BBO [6][7][8], and DECIGO [9][10][11][12] on the horizon, the prospects of GW detection are bright.
One of the most interesting early Universe events is the process called reheating.Inflation cools the Universe to a temperature T of zero due to its exponential expansion.Meanwhile the late-time Universe is well described by the standard model of cosmology, where the Universe is a hot thermal bath cooling due to the expansion of the Universe.Clearly sometime in between these two events the Universe must have gone from T = 0 to T > 0, a process referred to as reheating (RH).This makes RH a particularly special era in the history of the Universe, since in almost all models the temperature only ever decreases afterwards.Subsequent events such as entropy dumps only serve to cool the Universe more slowly rather than cause new RH events.Because RH likely only occurred once in our Universe, it is a unique and interesting event to study, and it is the target that we aim to elucidate.
In this article, we wish to probe how RH can be tested experimentally.Because RH is the only time when the temperature of the Universe increases, we are led to study signatures that arise from a period of increasing temperature.We therefore study the GW signature resulting from a heating phase transition (hPT) as opposed to the more commonly studied cooling phase transition (cPT). 1he GWs of a cPT (which we denote as cGWs) are chiefly generated by the colliding bubble walls, plasma sound waves, and plasma turbulence.The resulting signature is commonly characterized by the efficiencies (κ cPT ), the strength of the phase transition (α cPT ), the velocity of the bubble wall (v w,cPT ), the duration (β −1 cPT ), and the temperature of the phase transition (T cPT ).The same quantities, mutatis mutandis, characterize the GW signature of an hPT (hGW for short), with the addition of two more that parametrize reheating.In this work we take RH to begin at a Hubble scale H i corresponding to a time when the reheaton, the particle whose decay reheats the Universe, starts decaying with a rate Γ χ .The hPT parameters, e.g., β hPT , are related to their corresponding cPT parameters, e.g., β cPT , and they can be expressed in terms of each other up to O(1) numbers and RH parameters.
In general, the difference between an hPT and a cPT will be more than just a change in the parameters, as the dependence of the GW signal on both the PT parameters and frequency will change.For example, in a cPT plasma sound waves last for an entire Hubble time before the expansion of the Universe damps them.The result of this prolonged emission time is that the power in GWs is enhanced by a factor of β cPT /H cPT .In a heating PT, sound waves will instead be damped by the heating process, which adds total energy to the plasma damping the sound waves and giving a smaller enhancement of β hPT /Γ χ .While in this paper we focus more on the amplitude of the GWs, these same processes could also change the frequency dependence of the GW signal in interesting ways.In addition, novel plasma effects related to the restoration of symmetry that accompanies an hPT can enhance the amplitude of the GW signal coming from bubble collisions.
In this article, we study how an hPT percolates and describe how hPT and cPT parameters are related in a particularly simple model.Additionally, we study the special case when the heating and cooling PTs are both observable at future GW detectors.In these lucky scenarios, information about reheating can likely be obtained.Even if only an hPT were observed, it is possible that much could be learned from its frequency distribution.
In Sec.II we present a toy model with a phase transition, the object of our study.In Sec.III we investigate the details of how a heating phase transition completes, and review the cooling case.Section IV describes how GWs are generated by a heating and cooling phase transition, highlighting their differences.Section V details what is needed for an hPT to be found at future GW detectors.Finally, we conclude in Sec.VI, and supplement our results with four appendices.For the reader's convenience, in Table I we provide a list of our notation and some of the parameters most commonly used in the literature.

Variable Meaning
Hi Hubble expansion rate when reheating starts.Γχ Decay rate of the reheaton.ρχ, ρr, T Energy densities in the reheaton and radiation, and the temperature of the latter.Φ Higgs-like scalar field, whose spontaneous-symmetry-breaking potential drives the phase transition.⟨Φ⟩ b , ⟨Φ⟩ s Broken and symmetric minima of the temperature-dependent Φ potential.{µ, A, λ} Coefficients of the quadratic, cubic, and quartic terms in the Φ potential.0 < ∆ < 1 The useful combination ∆ = 4A 2 /(3λµ 2 ), which controls much of the physics of the phase transition.Tmax(tmax) Maximum temperature during reheating, and the time at which it is reached.
Tc(tc) Critical temperature and time: when the broken and symmetric phases are degenerate.T0(t0) Binodal temperature and time: when the symmetric phase becomes a maximum of the Φ potential.

T1(t1)
Spinodal temperature and time: when the broken phase no longer exists.

Tn(tn)
Nucleation temperature and time: when one bubble per Hubble patch is formed.TcPT(tcPT) Percolation temperature and time: when the cooling phase transition completes.T hPT (t hPT ) The same as above, but for a heating phase transition.
Γ/V Bubble nucleation rate per unit volume.

S
Euclidean bounce action of the bubble nucleation rate.

h(t)
Metastable volume fraction: the fraction of the volume of the Universe in the false vacuum.n b (t) Bubble number density.

RPT
Average distance between bubbles at percolation time.β Inverse duration of the phase transition.
Table I: Summary of the notation used in this paper to denote various quantities of interest.

II. A TOY MODEL WITH A PHASE TRANSITION
There are two main ingredients in the toy model that we consider.The first one deals with how RH takes place.
While there is a plethora of ways to achieve RH, as a representative toy model we consider a Universe whose energy density is entirely contained in a reheaton field χ, with RH proceeding via χ decays.For simplicity we assume that the daughter particles constitute an interact-ing dark sector (DS) with g * degrees of freedom, which quickly form a thermal bath. 2 Eventually this DS plasma reheats the visible sector (VS) containing the Standard Model (SM) via some portal interactions, whose form is irrelevant to our purposes and we thus leave unspecified.This scenario is characterized by two parameters: the Hubble scale H i at which χ starts to decay (which depends on the initial energy density ρ χ,i of the reheaton) and the decay rate Γ χ . 3Roughly speaking, the reheatondominated era lasts for a time ∆t χ ∼ 1/Γ χ equal to its lifetime, after the onset of its decay.
The other ingredient of our toy model is the field Φ responsible for the first-order phase transition.Φ is a component of a thermal bath after reheating, and will eventually generate GWs.As is standard, we base our model of Φ on the Higgs boson, where thermal corrections give rise to a cubic term in its potential and, consequently, to a first-order phase transition.We consider a finite temperature T potential [62][63][64] The parameters A, µ, λ, and T 0 have model-dependent values but we allow them to vary freely, in order to ensure the wide applicability of our results.These parameters come about from the coupling of Φ to other particles and the subsequent mass difference between the two phases ⟨Φ⟩ s ≡ ⟨Φ⟩ symmetric = 0 and ⟨Φ⟩ b ≡ ⟨Φ⟩ broken ̸ = 0.The particles and their interactions with Φ also play a crucial role in the dynamics driving bubble expansion.Note that the term −µ 2 T 2 0 Φ 2 /2 clearly corresponds to a tachyonic tree-level mass for Φ, which leads to the usual spontaneous-symmetry-breaking mechanism at zero temperature, with ⟨Φ⟩ 0 ≡ ⟨Φ⟩ b,T =0 = 3/λµT 0 and V 0 ≡ V (⟨Φ⟩ 0 ) = −3µ 4 T 4 0 /(8λ).There are several temperature values that are important in our model, as illustrated in Fig. 1.These temperatures are determined by the parameter combination ∆ ≡ 4A 2 /(3λµ 2 ).The first is T 0 , the temperature below which the symmetric phase ⟨Φ⟩ s ceases to be a minimum: for T < T 0 only the broken phase, with ⟨Φ⟩ b , is a minimum.The second temperature is the critical temperature T c = T 0 / √ 1 − ∆ at which the broken and symmetric phases have equal energies.Since degenerate minima are a requirement for there to be a first-order PT, we demand our potential parameters to satisfy ∆ < 1, which allows for the critical temperature to exist.For subcritical temperatures (T < T c ) the broken phase is 2 For the O(TeV) temperatures we consider, an interaction rate Γ th ∼ g 2 T is efficient enough to thermalize the DS plasma. 3The reheaton may itself be the inflaton in simple models such as m 2 χ 2 inflation.In this case we expect Γχ ∼ H i ∼ mχ for O(1) couplings, since inflation terminates when χ ∼ m Pl , where m Pl is the Planck mass [60].In other models such as hybrid inflation [61], the reheaton and inflaton are separate particles, and H i and Γχ are not necessarily related.
Figure 1: Finite-temperature potential V (Φ), with its symmetric (⟨Φ⟩ s = 0) and broken (⟨Φ⟩ b ̸ = 0) minima.For T < T0 (T > T1) the minimum corresponding to the symmetric (broken) phase disappears.At the critical temperature Tc both minima are degenerate.A cooling phase transition (blue arrow) from the metastable symmetric phase to the stable broken phase can take place for subcritical temperatures T < Tc, whereas a heating phase transition (red arrow) from the metastable broken phase to the stable symmetric phase can occur for supercritical temperatures T > Tc.
energetically preferred by the system, whereas for supercritical temperatures (T > T c ) the symmetric phase is more energetically favorable.Finally there is the temperature T 1 = T 0 8/(8 − 9∆) above which the broken phase ceases to exist: for T > T 1 , the only minimum is the symmetric phase.Note that if ∆ ≥ 8/9 the broken phase always exists.In the literature T 0 and T 1 are sometimes called the binodal and spinodal temperatures, respectively.If, however, A = 0 (∆ = 0 and T 0 = T 1 = T c ) there is no potential barrier separating the symmetric and broken phases and the phase transition is second order.Furthermore, thermal corrections to the self-energy of particles (what is commonly called "daisy resummation" [64][65][66][67][68]) may lower or erase the Φ potential barrier at high temperatures, thus weakening the strength of the first-order phase transition, or even negating it altogether.This puts a more stringent lower bound on A and therefore on ∆.As an estimate of this bound, we demand that these corrections be smaller than 50% at T c , which means ∆ > ∼ 0.27λ/µ 2 or A > ∼ 0.45λ.See Fig. 10 and Appendix B 3 for more details.

III. PHASE TRANSITIONS DURING REHEATING
In this section we detail the dynamics of phase transitions during reheating.While the cases of heating and cooling phase transitions are very similar, there are important differences.Because of this, we review some of the previous literature on first-order phase transitions, occasionally highlighting our new results.Of these, our discussion of heating phase transitions during reheating Figure 2: Thermal history of the Universe in our toy model, in which a reheaton energy density ρχ (blue) reheats a radiation energy density ρr (red) via decays.The reheating history is entirely determined by two parameters: the reheaton decay rate Γχ and the Hubble expansion rate Hi at the time ti = 0 at which these decays begin to take place.The curves in the figure only depend on the ratio of the two parameters; we chose Γχ = Hi as our reference point.The dashed black line represents the value ρr, c, where the radiation is at its critical temperature, here taken to be 80% of the maximum temperature, Tc = 0.8 Tmax.tc1 and tc2 are the two times when T = Tc.
takes the center stage.Although mentioned in passing in Refs.[69,70], a detailed study of the properties of hPTs taking place during RH has not, to the best of our knowledge, been published in the literature.
The thermal history that we consider is shown in Fig. 2. Initially all of the energy density is in χ, the temperature is zero and Φ is in the broken phase.As χ decays, a DS thermal bath develops and the temperature of the Universe increases as a function of time.During this "heating" era the energy density in the radiation grows linearly with time, ρ r ≈ Γ χ t ρ χ,i .Eventually the temperature grows larger than T c (i.e. the radiation density is ρ r > ρ r (T c ) ≡ ρ r,c ) and the broken phase in which the Universe finds itself becomes metastable.This means that a first-order hPT can now occur via the nucleation of bubbles of the stable, symmetric phase, which subsequently expand.At some point, corresponding to a temperature T hPT , these bubbles fill the entire Universe, which has now fully transitioned to the symmetric phase, and we can say that the hPT is completed.
Once the Universe reaches its maximum temperature T max at a time t max , it begins to cool down due to Hubble expansion.Famously, T max is larger than what is commonly known as the reheating temperature, the temperature of the plasma after the energy transfer from the reheaton [71][72][73][74][75].It depends on how much of the reheaton energy could be transformed into radiation before one Hubble time, roughly ρ r,max ∼ ρ χ,i min[1, Γ χ /H i ].After this maximum the temperature eventually falls below T c , and the symmetric phase of the Universe is now metastable.The previous process repeats but in reverse, with bubbles of the broken phase forming and growing, eventually filling up the Universe at some time when it is at a temperature T cPT , at which point it can be said that the first-order cPT is finished.
There are a few necessary conditions for PTs to take place.The first one is T c < T max , namely that the Universe should reheat above the critical temperature of the system.Otherwise the broken phase is always stable and the Universe remains in it throughout reheating, which means no PT takes place.It is also important that the PT completes before the metastable minimum disappears, i.e., T cPT > T 0 and T hPT < T 1 .If this is not satisfied, Φ will simply roll down to the stable minimum before a significant number of bubbles are formed, thereby preventing the production of a sizable amount of GWs.There are two times t c1 and t c2 when T = T c , which take place while the plasma is heating and cooling respectively; see Fig. 2. In addition, we denote by t 1 the time at which T = T 1 and the broken phase disappears, and by t 0 the time at which T = T 0 and the symmetric phase disappears.Therefore, the conditions stated above can be understood as follows: the hPT must complete at a time t hPT within the interval (t c1 , min[t c2 , t 1 ]), whereas the cPT must complete at a time t cPT within (t c2 , t 0 ).
The similarities and differences between cPTs and hPTs force us to make a brief housekeeping comment on notation: we use the generic subindex "PT" to denote that a quantity is evaluated at the end of a PT (t PT ), for statements that apply equally to the heating or cooling cases.However, if it is imperative to specify whether the quantity in question is evaluated at the end of an hPT or a cPT (t hPT or t cPT , respectively), we make this explicit.Throughout the rest of this section we review the dynamics of PTs, highlighting the differences between the well-studied cPTs and the novel hPTs.
The most common definition in the literature of the end of a first-order PT is the time t PT at which the fraction h(t PT ) of the volume of the Universe found in the metastable phase (metastable volume fraction for short) has been reduced to 1/e [37,63].This is often called the percolation time.It can be shown [37,63,76,77] that h(t) is given by4 where v w is the bubble-wall velocity, and Γ/V is the bubble nucleation rate per unit volume [37,62,78] Here S is the Euclidean bounce action associated with nucleating a Φ critical bubble of the lower-energy phase (see Appendix B 2 for its precise definition).Because at T c the broken and symmetric phases are degenerate, there is no free energy available to generate a phase transition from one to the other.As such, S diverges and Γ/V vanishes.Since S is a monotonically decreasing function of |T − T c |, the rate Γ/V grows exponentially from 0 as the temperature moves away from T c in either direction, i.e., both during the hPT and the cPT; this behavior can be seen in Fig. 3.As long as the parameter A controlling the height of the potential barrier that separates both phases is not too large, this exponential growth will generally guarantee a time at which the probability of nucleating one bubble in a Hubble space-time patch approaches 1.The time at which this happens is called the nucleation time t n [4,5,76,79] and it is roughly given by The exponential sensitivity of Γ/V to S means that in most cases the phase transition will occur in what is called the exponential nucleation regime [32,[80][81][82], which can be observed in Fig. 3. 5 In this regime, we can approximate the bounce action as S ∼ S 0 + S 1 t so that the nucleation rate is changing exponentially quickly.As a result the times t c < ∼ t n < ∼ t PT all take place in quick succession, and most of the bubbles of the new phase are nucleated towards the tail end of the process.In this regime, since t n ≈ t c , we can write and where in the last equality we have solved for S(t n ) from Eq. ( 5) by keeping the dominant exponential behavior in Eq. ( 3) and ignoring the prefactor.Once percolation is achieved and the PT ends at t PT , the bubbles of the new phase are large enough that they 5 The exponential nucleation regime is the norm in most of our parameter space, for both hPTs and cPTs.Nevertheless, there is another.Because Γ/V vanishes at both tc (since S(Tc) = ∞) and t 0 or t 1 (since S(T 0 ) = S(T 1 ) = 0), Rolle's theorem guarantees that the function Γ/V has a maximum as a function of time.In the region of parameter space for which this maximum is comparable to the Hubble rate (i.e., Γ/V < ∼ H 4 ), the so-called simultaneous nucleation regime [81,82] takes place.In this regime most bubbles of the new phase are nucleated at the time when Γ/V is at its largest.This occurs more easily when Tmax is reached during an hPT, which happens for large bounce actions.For more details about this regime, we direct the reader to Appendix C 4.
Figure 3: Example of the bubble nucleation rate Γ/V divided by H 4 as a function of time t.The red curve corresponds to nucleation during an hPT (t ∈ (tc1, min[tc2, t1]), whereas the blue curve corresponds to nucleation during a cPT (t ∈ (tc2, t0)).The dashed lines indicate the unrealized evolution of Γ/V, since they correspond to times after the phase transitions have actually finished.The horizontal line marks Γ/V = H 4 .The vertical lines correspond to different times of interest.For this plot we chose the parameters corresponding to the starred benchmark point in Fig. 6, namely Tc = 1 TeV, Hi = 2 × 10 −15 TeV = Γχ (i.e., Tc = 0.8 Tmax), and g * = 10, the potential coefficients {µ, A, λ} = {1, 0.72, 1} (∆ = 0.7), and a bubble-wall speed in the cPT of vw,cPT = 0.05.
are very close to each other and begin to collide.It is these collisions and the subsequent behavior of the system that give rise to gravitational waves.The quantity of interest is the mean bubble separation scale at percolation R PT , defined in terms of the bubble number density n b (t) as follows: with From R PT and the bubble-wall velocity v w one can obtain the characteristic time scale β −1 of the PT 6β ≡ (8π) It can be shown that in the exponential nucleation regime [32,[80][81][82] which is the definition of β more commonly found in the literature.However Eq. ( 10) has a wider range of applicability and it is more closely related to the peak frequency of the GW spectrum [81,83].Because of this, we use Eq.(10) in our results, which are numerically calculated.
There is one last significant difference between cooling and heating phase transitions, regarding the manner in which their respective bubbles expand.Indeed, while during a cPT bubbles generally reach a constant subluminal velocity, in an hPT they instead typically enter a runaway regime, in which their wall velocity quickly approaches the speed of light (v w → 1).Below we justify this claim in a more or less quantitative manner, leaving a more detailed discussion and a description of the runaway parameter space to Appendix C 2. Finally, we would like to caution the reader that the growth of bubbles in the presence of a plasma is governed by very complex dynamics, and it is the subject of ongoing research [5,[84][85][86][87][88].
Once nucleated, the bubbles of the new stable minimum grow due to the free energy difference inside and outside of their wall.We can determine whether these bubbles run away by considering a relativistic bubble wall (where friction is maximized and antifriction is minimized) and asking if the net pressure acting on it is pushing outwards, driving the wall ever faster.Mathematically, written in terms of the total force per unit area P tot acting on the bubble wall, this runaway condition is given by where ∆V 0 ≡ V 0, out − V 0, in = sign(T c − T ) |V 0 | is the zero-temperature potential difference between the outside and the inside of the bubble (for cPTs (hPTs) the zero-temperature broken minimum V 0 < 0 is inside (outside) the bubble), and ∆P T is the pressure difference produced by the plasma, which is the same as the leading mean-field contribution to the thermal potential [84]: , in is the i-th particle's mass difference between the outside and the inside of the bubbles, N i accounts for its degrees of freedom, and c i = 1 (1/2) for bosons (fermions).The second equality in Eq. ( 13) stems from the definition of µ in terms of the particle interactions (see Eq. (B6)) and from the fact that the masses of the particles depend on their Yukawa couplings to Φ.The sign takes into account that Φ = ⟨Φ⟩ b inside cPT bubbles (for which T < T c ) , while the opposite is true for hPTs.
This sign of ∆P T can be readily interpreted in terms of momentum conservation in the bubble rest frame [84,86] with the help of Fig. 4. In this frame the bubble wall is static, while an incoming flux of plasma particles is moving towards it with velocity −v w .For cPTs (with T < T c ) the massless plasma particles transmitted through the wall gain a mass m 2 i, in > 0, thus lowering their momentum.Momentum conservation then implies that the bubble wall needs to make up for the missing momentum by moving along the direction of the incoming flux of plasma particles, which corresponds to a slowing down of the bubble expansion.For hPTs (with T > T c ) the process is the opposite: the transmitted particles lose their outside mass m 2 i, out > 0 and thus gain momentum, which the bubble wall needs to balance out by accelerating in the direction opposite to the incoming particle flux.
As a result, this pressure difference acts as a plasma friction on the bubbles of a cPT, or as a plasma antifriction on the bubbles of an hPT.The former case has been much discussed in the previous literature, while the latter is presented in detail, to the best of our knowledge, for the first time in this work. 7Intuitively, the massless plasma particles outside the bubbles of a cPT experience a Φ potential barrier at the wall.This means that these particles can bounce off of the bubble walls or lose momentum upon entering the bubble, thereby exerting a friction on the bubbles and slowing down their growth.This friction typically ends up balancing the force driving the bubble expansion, and a constant v w , often subluminal, is reached [84][85][86]88].Heating phase transitions, on the other hand, have almost exactly the opposite behavior.During an hPT the plasma particles are massless in the interior of the bubbles, where ⟨Φ⟩ s = 0.As a result, it becomes energetically favorable for these particles to get inside the bubbles: the plasma is "sucked in" by them.As the plasma particles pass through, they transfer their energy to the bubble walls, accelerating them, and leading to a bubble runaway regime.
Owing to the fact that T cPT ∼ T c ∼ T hPT for generic values of the potential parameters, we can see that the bubble runaway conditions for cPTs and for hPTs described in Eq. ( 12) are reflections of each other: when the condition is satisfied in one case, it will typically not be satisfied in the other.This makes sense because for runaways to take place in cPTs, the friction of the plasma acting on the bubble has to be small, whereas in hPTs the antifriction has to be large; ultimately it is µ, which parametrizes the strength of the interaction between the plasma particles and the Φ field, that determines the size of both friction and antifriction.Therefore, if a given µ produces enough plasma friction to cause the bubbles of a cPT to expand at a constant wall speed, it will also cause that very same plasma to exert instead an antifriction on the bubbles of the hPT, accelerating them into a runaway.Because the plasma is transferring its energy into the bubble walls in order to accelerate them, this means that the energy available in an hPT for GW production 7 Phase transitions during reheating were briefly mentioned in Refs.[69,70] in the contexts of both the Standard Model and inflaton models with dynamical decay rates respectively.the case of a cooling phase transition, where the plasma particles gain mass upon entering the bubble (blue region), thereby slowing down (short green arrows to the left of the bubble wall).Momentum conservation means the plasma exerts a friction force on the wall (blue arrows), which moves it to the left; in the plasma rest frame the bubble wall is slowed down.Bottom: the case of a heating phase transition, where the particles lose mass inside the bubbles (red region) and thus accelerate (long green arrows), thereby accelerating the bubble wall (red arrows).
in bubble collisions can be very large.

IV. GRAVITATIONAL WAVES DURING REHEATING
The frequency spectrum of a SGWB is typically described in terms of the fraction of the total energy density of the Universe found in GWs, per frequency e-fold: (dρ gw /d ln f )/ρ tot .Denoting with an asterisk all those quantities evaluated at the time t * at which the GWs are produced, we can find the SGWB spectrum today by accounting for the difference in ρ tot and the redshift as follows [5,33,37]: with In the first equality the prefactor F * accounts for the radiation-like redshifting of the GWs with the scale factor a * and for the ratio of the total energy densities at t * and today,8 ρ tot, * /ρ crit = H 2 * /H 2 0 .In the last equality we have written the energy density in GWs in terms of the energy density ρ s, * of their source: Here G ∼ H 2 /ρ tot is Newton's constant, τ ∼ β −1 is the typical time scale of GW production during a PT, and κ is an efficiency factor that quantifies how much of the energy ρ s, * in the source goes into GWs.The factors N and S(f ) account for an overall normalization and spectrum, respectively, and they may depend on other phase transition parameters, such as the bubble-wall velocity v w .
The sources of GWs during a PT are typically of three kinds: bubble collisions, sound waves, and magnetohydrodynamic turbulence.The first one involves the energy stored in the Φ bubble walls, while the last two come from the response of the plasma to the nucleation and percolation of the bubbles of the new phase.Each of them has different frequency spectra and dependences on the PT parameters.Their precise form and relative contribution to the overall GW signal is the subject of ongoing research (see Refs. [4,5,37] for reviews).In cPTs, the case most commonly studied in the literature, runaway bubbles are not typically expected and the contribution from sound waves tends to be the largest [37,89]. 9On the other hand, as discussed in the previous section, the same plasma exerting a friction on cPT bubbles will instead exert an antifriction on hPT bubbles, accelerating them.Because of this, most of the energy is stored in the bubble walls, and we expect that in hPTs the dominant GW contribution comes from the collisions of runaway bubbles.Since we are interested in the detectability of GWs from reheating as a proof of concept, and since the contribution from turbulence is the most uncertain [4,5,37], we will not consider it throughout the rest of this paper, and we focus instead on GWs coming from bubble collisions for hPTs and from plasma sound waves for cPTs.

A. Bubble-wall collisions
In this subsection we briefly discuss the gravitational waves generated by bubble-wall collisions.We study this signature in the context of hPTs, where it is the dominant source of GWs, whereas they are typically subdominant in cPTs [37,89].The GW spectrum of a bubble-wall collision is typically calculated numerically with the addition of the envelope approximation [30,31,36,81,[90][91][92][93][94], which approximates the bubbles as an expanding set of infinitely thin shells that disappear when the transition completes.In these numerical calculations, the Hubble expansion is typically neglected as the PTs being studied are assumed to complete very quickly.In the hPTs that we consider, we make a similar assumption.However, as the injection of energy is dictated by Γ χ , we instead assume that the phase transition completes quickly relative to both Γ χ and H hPT .Under this assumption, the numerical results apply equally well to cPTs and hPTs.
The GW spectrum is found numerically to be10 Below we briefly explain the various parameters describing the GW spectrum from bubble collisions, with the aid of Eqs.(A30)-(A34).For illustrative purposes we focus on the case of a short era of reheaton domination (∆t for which these parameters take simple forms.For more details, as well as the case of an arbitrary duration of reheaton domination, we refer the reader to Appendix A 3. 11t * : The time at which the GWs are generated.It corresponds roughly to when the bubble collisions take place, which in turn is very close to the percolation time t hPT .We therefore take t * ≈ t hPT .
T * : The temperature of the plasma at the time when the GWs are generated.From the previous paragraph, T * ≈ T hPT > T c .Except for very fine-tuned Φ-potential parameters, the hPT and critical temperatures are similar in scale, T hPT ∼ T c .
H * : The hPT typically finishes at a time t hPT during the reheaton-dominated era, so R χ : The energy density of the reheaton over the radiation density, at t * : R χ ≡ ρ χ, * /ρ r, * ≈ ρ χ,i /ρ r,hPT = 3H 2 i m 2 Pl /ρ r,hPT .a * : The scale factor at which the hPT takes place.If the reheaton-domination era is shorter than one Hubble time, it is given roughly by a * ∼ 8 × 10 −17 (1 TeV/T max ).A longer reheaton domination corrects this expression with a factor that depends on Γ χ /H i .There is also a mild dependence on the degrees of freedom of the dark and visible sectors.
F * : This redshift factor depends on the photon and total energy densities today, as well as the duration of the reheaton-dominated era.It is approximately given by F * ∼ 4 × 10 −5 for short reheaton domination.
α, κ Φ : The vacuum energy in the scalar Φ compared to the energy in the radiation; and the corresponding efficiency factor that quantifies how much of it goes into the bubble walls (essentially gradient energy of the Φ field), determined by Eq. ( 12) in the runaway bubble regime: where In our work we never consider vacuum domination but only reheaton domination and radiation domination, i.e., α < 1, R χ .Note that, up to a minus sign, Eq. ( 21) is the same expression quantifying how much energy goes into the bubble wall for runaway cPTs [4, 5, 84-86, 95, 96].The sign difference, seen in Eq. ( 12), stems from the fact that the direction of the PT is the opposite in hPTs than in cPTs.As discussed in the previous section this sign means that, while in most of their parameter space cPT bubbles do not run away (reaching a subluminal terminal bubble-wall velocity) and κ Φ is consequently a number much smaller than 1, for hPTs we have instead α ∞ > α for most of their parameter space, because the plasma antifriction makes the hPT bubbles run away, greatly increasing the energy stored in their walls.Note that our definition of α here is directly borrowed from previous literature, which has only dealt with cPTs.It is therefore not the most natural way to parametrize the energy density available in an hPT, which is not As a result of this definition, α ∞ > α in runaway hPTs means that κ Φ can be much larger than 1.See Eq. (C6) and Fig. 10 as well as Appendix C 2 for a more detailed discussion on the runaway condition.
v w : As discussed in the previous section, for most of the region of parameter space in which a strongly firstorder phase transition takes place, the plasma exerts an antifriction on the bubble walls during an hPT, leading to a runaway regime and therefore v w → 1.The parameter space leading to a runaway hPT can be found in Fig. 10 in Appendix C 2.
β : The inverse of the characteristic time scale τ of the PT, given by Eq. ( 11) in the exponential regime.Up to O(1) factors S(t hPT ) < ∼ S(t n ), while d ln T /dt < ∼ Γ χ ρ χ,i /(4ρ r,hPT ) from the equations governing RH (Eq.(A2)).From this and our discussions in Appendices B 4 and C 4, we find where B. Plasma sound waves For the cPTs that we consider, the main source of GWs are the sound waves.The SGWB created by sound waves is found to be Note that, unlike the GWs sourced by bubble collisions, the spectrum from sound waves scales like one power of H * /β rather than two.This is because the fluid bulk motion sourcing the GWs lasts for about a Hubble time, longer than the PT duration 1/β [4,5]. 12We continue focusing on the simple case of a short reheaton-dominated era, in which case the cPT typically occurs during radiation domination.t * , T * , H * , R χ , a * , F * : The GWs are generated at a time t * shortly after the cPT is completed at t cPT ; we therefore simply assume t * ≈ t cPT .For typical potential parameters T * ≈ T cPT < ∼ T c .The values of H * and R χ at this time will be smaller than their hPT counterparts due to both the expansion of the Universe and the reheaton decays.If the cPT completes firmly during the postreheating radiation-dominated era then R χ ≈ 0, while cPT /(90m Pl ), a * ≈ 8 × 10 −17 (1 TeV/T cPT ), and F * ≈ 4 × 10 −5 .However, the cPT may also take place during reheaton domination.Estimates of these quantities for both cases are listed in Eqs.(A30)-(A34) in Appendix A 3.
v w : In a typical cPT, the plasma exerts a friction that balances out the expansion force of the new phase, which leads to the bubble walls moving at a constant, typically subluminal speed.Complex model-dependent dynamics govern the expansion of bubbles in a thermal plasma, and their terminal velocity can in principle be derived from these; in this work we simply bypass the issue by taking v w in a cPT to be a free parameter.
α, κ sw : These quantify the energy released during a cPT and how much of it goes into the bulk motion of the fluid.There are several ways to compute these quantities in the literature.The simplest one relies on the bag model, where the vacuum energy is used, and thus α = |V 0 | /ρ r [85].Recent theoretical and numerical developments advocate instead for the use of the trace of the stress-energy-momentum tensor [5,[99][100][101].However, in this latter model a more thorough knowledge of the plasma fluid is necessary in order to compute κ sw , (e.g., the plasma sound speeds of both the symmetric and broken phases), knowledge that we do not have.We simply use the bag model and the corresponding fits to κ sw (α, v w ), conveniently provided in Ref. [85].For the subluminal cPT bubble-wall velocities we consider in this work, κ sw ≈ 5αv 6/5 w .β : Also given by Eq. ( 11) in the exponential regime.Since the temperature is cooling due to the Hubble expansion of the Universe, d ln T /dt = −H * .It can then be shown that during radiation domination where the estimate for the general case is shown in Eq. (C17) (see Appendices B 4 and C 4 for more details).Finally, we say a few words about GWs from sound waves in hPTs.Since the runaway regime is common in hPTs, the plasma puts energy into accelerating the bubble walls and therefore bubble collisions dominate GW production, with sound waves most likely playing a subdominant role.It is difficult to tell, without dedicated numerical simulations, whether the GWs from sound waves in hPTs have a different frequency spectrum or amplitude than their cPT counterparts.This certainly seems possible: a first difference between both scenarios is that while in cPTs the bubble walls inject energy into the plasma, leading to radially outward fluid bulk motion and subsequently sound waves, in hPTs energy is removed from the bath and put into the bubble walls, with the sound waves moving in the opposite direction to the walls.Yet another difference is the duration of the sound waves sourcing the GWs.For a cPT, the Hubble expansion eventually damps the sound waves.For an hPT, the reheaton deposits its energy in the thermal bath at a rate ∼ Γ χ ρ χ /ρ r , thereby increasing the radiation energy density and damping the sound waves.The duration of the sound waves in hPTs is then ∼ ρ r,hPT /(Γ χ ρ χ ) As a result, we expect that the amplitude of the GWs from sound waves in an hPT will differ from that in a cPT by a factor of ∼ (H cPT /Γ χ )(ρ r,hPT /ρ χ ), accounting for their shorter relative duration.

C. Comparison between the GW spectra from cooling and heating phase transitions
The above discussion allows us to compare the GW spectra for hPTs and cPTs, which are dominated by collisions and sound waves, respectively.In the simplest case of a short reheaton-dominated era the ratios of both take their simplest forms: where we denote the cPT bubble-wall speed by v w,cPT , we take T PT ≈ T c , and the numerical benchmark values we show are typical of our parameter space.We remind the reader that these expressions were derived only for a shortly lived reheaton-domination era, that they are to be taken only as heuristic, and are to be trusted only as order-of-magnitude estimates.
From the equations above it can be seen that, barring a fine-tuning of the Φ potential parameters in order to get vastly different κ Φ and |d ln S/d ln T |, the amplitude of the SGWB from both hPT and cPT can be of the same order of magnitude for a modest coincidence between ρ r,c and ρ χ,i , of about a couple of orders of magnitude (or equivalently a coincidence between T c and T max of a factor of a few).This coincidence has to be more severe the larger the hierarchy between Γ χ and H i .The peak frequencies of the spectra can also be close to each other, with f hPT slightly smaller than f cPT for subluminal cPT bubble-wall velocities and the aforementioned coincidence between ρ r,c and ρ χ,i (since of course T c < ∼ T max ).Increasing the Γ χ /H i ratio inverts the order of the peaks, with f hPT eventually becoming bigger than f cPT .Furthermore, in the limit of a long reheaton-dominated era (Γ χ ≪ H i ), the hGWs are typically quieter than the cGWs.This is both because the ratio ρ r,c /ρ χ,i is smaller (since ρ r,c ≤ ρ r,max ∼ ρ χ,i Γ χ /H i in this limit), and because more time has elapsed between t hPT and t cPT , which means that the relative redshift suppression (F hPT /F cPT ) between both spectra becomes more significant.
Example cPT and hPT SGWBs, for a benchmark pa-rameter space point, are shown in Fig. 5.The characteristics described above are easily appreciated in this figure .As this example illustrates the same detector, in this case the future BBO probe [6][7][8], can observe both spectra.In the lucky case that both signals are loud and have sufficiently separated peak frequencies, so as to not have one of them buried under the other, one could then in principle distinguish between them, identifying which belongs to the cPT and which to the hPT.Then the correlations between them, of which Eqs. ( 31) and ( 32) are examples, would allow one to extract both the scale of reheating and the reheaton decay rate.

V. PROSPECTS AT FUTURE GRAVITATIONAL-WAVE DETECTORS
In this section we explore the visibility of GWs generated by cooling and heating PTs in future detectors.We focus on the upcoming BBO experiment [6][7][8], as it is the most relevant detector for our parameter choice.We quantify visibility in terms of the signal-to-noise ratio (SNR) of the GW spectra.
To obtain the SNR of the PT GWs we employ the peak-integrated sensitivity (PIS) curves Ω PIS (f ) introduced in Ref. [102].We compute the SNRs of GWs, originating mostly from bubble collisions in hPTs and from plasma sound waves in cPTs (abbreviated as hGWs and cGWs respectively), by simply comparing the amplitudes The red curve represents the GW spectrum arising from bubble collisions during the hPT (hGW) and the blue curve represents the GW spectrum originating from sound waves in the cPT (cGW).Both spectra are detectable by the future GW probe BBO [6][7][8], as shown by the PIS curves [102] (black lines, solid for GWs from bubble collisions, and dashed for GWs from sound waves).For this plot we chose the parameters corresponding to the starred benchmark point in Fig. 6, namely Tc = 1 TeV, Hi = 2 × 10 −15 TeV = Γχ (i.e., Tc = 0.8 Tmax), and g * = 10, the potential coefficients {µ, A, λ} = {1, 0.72, 1} (∆ = 0.7), and a bubble-wall speed in a cPT of vw,cPT = 0.05.
of the GWs at their peak with Ω PIS (f ): where t obs is the observation time, f peak is the frequency at which Ω GW is the peak, and "bc" ("sw") corresponds to the GWs from bubble collision (sound waves).
The results for hGWs and cGWs are shown in Fig. 6 in terms of the T c /T max -∆ parameter space (and the corresponding values of H i and A, respectively).The benchmark point used in previous figures of this paper, such as in Fig. 5, is marked with a star in Fig. 6.The SNR contours for t obs = 1 year for hGWs (cGWs) at BBO are shown in red (blue), where the increasing opacity indicates larger SNR.We have fixed T c = 1 TeV, Γ χ = H i , g * = 10, {µ, λ} = {1, 1}, and v w,cPT = 0.05.This specific choice of values has only a modest impact on our results, and the GW features described in this section are generic.We direct our reader to Appendix D for a more detailed study of the parameter space.
Note that ∆, which controls the height of the potential barrier separating the broken and symmetric minima of the Φ potential and thus the action S, strongly determines the strength of the GWs.For both hPTs and cPTs a larger ∆ makes d ln S/d ln T smaller, which increases the duration β −1 of the PT.Eventually, however, sufficiently large values of ∆ will kill the GW signature by making the PT impossible, as clearly seen in the region above the dotted line.Indeed, this region corresponds to those points with Γ(t max )/V < H(t max ) 4 .Since at t max Γ/V is at its largest (because the temperature is at its maximum and thus the action S is at its minimum), no bubbles are produced within a Hubble patch in this region.This means that the Universe remains in the broken phase throughout all of reheating, never transitioning to the symmetric phase (via an hPT), and therefore never coming back to the broken phase (via a cPT).Thus no PT takes place and therefore no appreciable GWs are produced.We nevertheless show the continuation of the cPT contours in this empty region for illustrative purposes.
The parameter T c /T max (which can be turned into ρ r,c /ρ χ,i for a given Γ χ /H i ratio) also has a crucial impact on the visibility of the hGWs.A large hierarchy between T c and T max means that the time elapsed between t hPT and t cPT (which are on opposite sides of the reheating curve of Fig. 2) is also large.As such, the GWs associated with the hPT are produced much earlier than those generated during the cPT, and are therefore more redshifted and correspondingly quieter.Thus their signal falls outside of the BBO sensitivity window.The combination of this effect and the one controlled by ∆ described in the previous paragraph gives the BBO-visible hPT region its characteristic crescent shape.The cGWs have a milder dependence on T c /T max .A strong coincidence between these two temperatures means ρ r represents a larger share of the total energy density, which makes the GWs louder.On the other end, the more different T c and T max are, the later the cPT takes place.For sufficiently large hierarchies the cPT occurs squarely during radiation domination and its GWs become insensitive to reheating.This is shown in Fig. 6 as the insensitivity of the cPT contours to low values of T c /T max .
The region enclosed by the green contour corresponds to points where the peaks from both cGWs and hGWs can be distinguished.This "double peaks region" is therefore defined as the parameter space where the peak amplitude of hGWs is larger than the amplitude of the cGWs at that same frequency, and vice versa.These points are potentially the best ones in terms of how much we could learn from these GWs.Indeed, observing both GW peaks in a GW detector could allow us to extract the most information about reheating by studying the correlations between cooling and heating PTs, as we have attempted to do in this paper.The region located above the dashed line has T max < T 1 , which means that the plasma never reaches T 1 and therefore ⟨Φ⟩ b never disappears.The dark gray region is the parameter space where the plasma antifriction during the hPT is not strong enough to cause the bubble walls to enter a runaway regime.Finally, the light gray region corresponds to those points where the daisy contribution becomes so large that the cubic term in Eq. ( 1) is severely suppressed and there is no strongly first-order PT (SFOPT).For more details on both grey areas see Appendices C 2 and B 3, respectively.

★ ★
Figure 6: SNR contours for 1-year observation time of the SGWB spectra generated during reheating by hPT bubble collisions (red) and cPT sound waves (blue), for the upcoming BBO detector [6][7][8], as a function of the ratio Tc/Tmax and the potential parameter ∆.The corresponding values of Hi and A are also shown.The star corresponds to the benchmark point used in the previous figures.Points within the green contour have a total SGWB with double peaks (from both hPT and cPT).The ⟨Φ⟩ b minimum never disappears in the region above the dashed line, since Tmax < T1.In the space above the dotted line Γ/V < H 4 at tmax, which means that no PT takes place and no GWs are generated.In the dark grey region there are no runaway hPT bubbles, whereas in the light grey region the daisy contributions to the thermal potential prevent an SFOPT.

VI. CONCLUSIONS
In this article, we explored phase transitions that occurred when the Universe was heating up, a process called reheating.Reheating is a unique and interesting cosmological event in the history of the Universe, about which we currently know nothing.If a phase transition occurred during reheating, then the resulting GW spectrum would carry information about the reheating process to us and teach us something about this exciting era of the early Universe.
We discussed how the GW signature of heating phase transitions depends on the phase transition and reheating parameters.In some optimistic scenarios, such as the one depicted in Fig. 5 and in the green region of Fig. 6, it would be possible to see the same phase transition in both its heating and cooling directions.If this came to pass then, by cross correlating the observed spectra and our knowledge of the phase transition in both directions, one could learn more about the process of reheating.For example, the time scales β −1 of the two phase transitions, which appear in the peak frequency of their respective GW spectra, are related by Γ χ /H cPT up to O(1) numbers.Therefore, by comparing the measured values of the peaks, one would be able to obtain the decay width of the reheaton.Many other features are correlated between heating and cooling phase transitions.Finite bubble-wall velocities in the cooling case tend to be correlated with runaway bubbles in the heating case, and vice versa if the roles were switched.These correlations offer an opportunity to extract information about reheating.
While not explored in our article, it would be extremely interesting if the frequency spectrum of the GWs produced during a heating phase transition were radically different from what has been studied in the literature.In this paper we argued that the main contribution to the GWs in this case comes from bubble collisions, and assumed that the envelope approximation can be used to describe its spectrum.This assumption, while justified (see Sec. IV), needs to be corroborated by dedicated numerical simulations.Furthermore, it is entirely possible that plasma sound waves represent a non-negligible source of GWs during a heating phase transition, and it is reasonable to expect that they would be different at the O(1) level from those in the cooling case, since cooling phase transitions tend to inject energy into the plasma while heating ones tend to remove it.Additionally, the sound waves of cooling phase transitions are damped by the Hubble expansion, whereas those from heating transitions are damped by the injection of energy from reheaton decays.If differences of this sort between both spectra were to be firmly established, it is possible that the detection of GWs generated by a heating phase transition would by itself be enough to infer detailed information about reheating.Future research along this direction is both necessary and of great interest to anyone attempting to uncover the physics of reheating.have taken a(t) = 1 today.The decomposition into three products is useful as some of these factors will be evaluated numerically.The first factor is simply given by the number of e-folds between those two times: The second term depends on the nature of the portal interaction that allows the DS to reheat the VS.For simplicity we simply assume that the VS is populated by the DS while both baths are relativistic.From energy conservation considerations, we can write where T denotes the temperature of both baths, T X = T (t X ), and g vs is the VS degrees of freedom reheated by the DS.Throughout this paper we assume that the entire SM is reheated, and therefore take g vs = 106.75.Note that t rd is indeed arbitrary: an earlier time (and consequently hotter T rd ) yields a smaller Eq.(A4) and a larger Eq.(A5) in the same proportion, and consequently the same Eq.(A3).Finally, the last factor can be derived from the usual entropy considerations.Assuming no remaining DS light relics, this factor is where g s0 = 3.94 and T γ0 = 2.7255 K ≈ 2 × 10 −4 eV are today's entropy degrees of freedom and temperature, respectively.
The final exact expression is then For an easier comparison to previous work, we can multiply and divide by g where [a * ] denotes the value of a * for cPTs taking place during RD, the case commonly studied in the literature.We can use our results to find F * , the prefactor responsible for the redshift-related suppression of the GW spectra (see Eq. ( 15)), and compare it to its value [F * ] for cPTs during RD [5]: (A13)

Analytic expressions for reheating
Depending on the reheaton decay rate, one can classify the reheating history into two relatively easy to study broad classes: those with γ ≫ 1, or those with γ ≪ 1.For each of these cases we can find analytic expressions for various quantities of interest, such as the maximum temperature reached during reheating or the scale factor at a given temperature.We present these expressions here in a quick-and-dirty fashion.A more thorough job can be done, and is included in the explanatory notebook 01 reheating.nb of graphare.
a. Fast reheaton decays: γ ≫ 1 Before discussing the consequences of extremely fast reheaton decays, γ ≫ 1, it is worth mentioning under what circumstances this can occur.The simplest example is simply low-scale hybrid inflation [61].In hybrid inflation, the inflaton drives inflation and a separate waterfall field abruptly ends inflation after the appropriate number of e-foldings.Reheating at the end of hybrid inflation is typically extremely efficient, e.g., tachyonic reheating completes within a single oscillation [103].Because the mass of the waterfall field (which in this context is the reheaton) is much larger than the Hubble rate, this results in γ ∼ m waterfall /H ≫ 1 so that we are in the limit of a fast reheaton decay.Alternatively, one can simply give the waterfall field a Yukawa coupling to a fermion and have it quickly decay in a more traditional manner.
If one finds low-scale inflation unappealing, an alternative manner in which to arrive at the fast reheaton decay scenario is to utilize either another phase transition or another slowly rolling scalar that controls the mass of the decay products of the reheaton.In the early Universe the mass of the decay products is large compared to the reheaton mass kinematically forbidding the decays.Only later will the mass of the decay products become small enough to allow for the reheaton to decay quickly.Leaving aside any specific ways to realize fast reheaton decays, we simply focus on its phenomenological consequences.
In the and the χD era is very short.ρ χ,i is quickly and efficiently transformed into radiation energy density, and the Universe enters into an RD era.Therefore, the maximum amount of radiation is ρ r,max ≈ ρ χ,i , and the time t max at which this occurs is solidly within RD and before one Hubble time has passed: i , a max ≈ a i ≡ a(t i = 0) (see Fig. 7).
Because an hPT can only take place while the temperature is increasing, t hPT < t max and a hPT ≈ a max ≈ a i .Without tuning T hPT (or T c ) to be too close to T max , the hPT will occur during χD, and ρ tot,hPT ≈ ρ χ,hPT ≈ ρ χ,i .On the other hand, since T cPT < T max necessarily, the cPT will always occur during RD.Of course, as is typical during RD, the energy density scales like and χD is long.During most of this era the decay rate, but not the Hubble expansion rate, can be approximated as negligible for the purposes of estimating the evolution of the reheaton energy density.Since the reheaton is matter-like then ρ χ ∼ a −3 and a ∼ x 2/3 , which means that and thus with a χR and ρ χR denoting the scale factor and reheaton energy density at the time of reheaton-radiation equality. 13The numerical coefficients inside the square brackets can only be obtained by solving the differential equations Eqs.(A1) and (A2) analytically in the γ ≪ 1 limit.The time of increasing temperature lasts for little less than one Hubble time (see Fig. 7); consequently, the hPT occurs squarely during χD.As discussed before, the radiation energy density grows linearly during this time (ρ r ≈ ρ χ,i γx), which allows us to estimate the maximum energy density in radiation as given roughly by where once again the number in the brackets can only be obtained by solving Eqs.(A1) and (A2) analytically and making use of the fact that ρr (t max ) = 0 by definition.Because t max < ∼ H −1 i , then a max ≈ a i as well.After one Hubble time the expansion of the Universe is felt and ρ χ ∼ a −3 , as stated in Eq. (A16).Because the Universe is in χD (ρ χ ≫ ρ r ) the reheaton decays are still very much relevant to the evolution of ρ r .However, the Hubble expansion wins and the temperature no longer rises with time, as clearly seen in Fig. 7. Equation A2 can be solved in this regime to find that ρ r ∼ a −3/2 : where we have finally dropped the O(0.1-1) numerical coefficients, irrelevant to the level of precision with which we are working.It is convenient to note that during χD ρ χ /ρ r ∼ a −3/2 : Of course, once χD ends and the Universe becomes RD, we have ρ r ∼ a −4 : The previous discussion allows us to find the important quantities R χ, * , H * , a * , F * , and d ln T /dt.The first four appear in the formulas for the SGWB from PTs (Eqs.( 16)- (19) and Eqs. ( 25)-( 28)), while the last determines the duration of the PT, β −1 (Eq.( 11)), which determines the peak frequency of the GW spectra.We can estimate these quantities both for hPTs and cPTs, in both the γ ≫ 1 and γ ≪ 1 limits.This task is most easily accomplished by defining a handful of useful parameters with a clear physical intuition.The first is what we call the reheating efficiency parameter ε rh , which quantifies how much of the initial reheaton energy density is transformed into radiation: It is clear that the first argument of min[X, Y ] is picked when γ ≫ 1, while the second when γ ≪ 1, which is as we found Eqs. (A14) and (A18).
Noting that in Eqs.(A8) and (A11) a * and F * are defined relative to an arbitrarily chosen benchmark time during RD; we now make this choice concrete.Since for γ ≫ 1 a max is well within RD we make this our benchmark for this limit.For γ ≪ 1, assuming a sudden transition between χD and RD, we can make a χR our choice.It is interesting to point out that a χR < a max for γ ≫ 1, while the opposite is true for γ ≪ 1.We thus choose where once again the first argument of max[X, Y ] occurs for γ ≫ 1 and the second for γ ≪ 1.We can then define the relative scale factor between t max and t rd as well as the redshift dilution in the radiation energy density between these two times Of course, ε rh , a rd , A rd , and D rd , can all be written in terms of γ and each other.But they represent quantities with distinct physical meanings, and thus it is more expedient to interpret them separately.
Based on these definitions, we can write the following expressions for hPTs (which only occur during χD and thus a hPT ≈ a i ≈ a max ), and cPTs (which can occur during either χD or RD): Putting everything together, we list A rd D −1/4 rd cPT; (A32) where is the usual Hubble expansion rate during RD.One can easily convince oneself that the first argument in min[X, Y ] and max[X, Y ] in the above expressions is picked when the cPT occurs during RD (both for γ ≫ 1 and γ ≪ 1), while the second argument is picked only for cPTs taking place during χD (which can only happen for γ ≪ 1).Right: bubbles for supercritical temperatures T ∈ (Tc, T1), i.e., λ ∈ (1, 9/8).

Critical bubbles and bounce action
A phase transition between the two minima of V (Φ) can be described by the nucleation and subsequent growth of "bubbles" of the new, stable minimum inside a space filled with the old, metastable minimum.These bubbles are created by thermal fluctuations and they correspond to a field configuration Φ bub that interpolates between the new minimum at the center of the bubble, and the old minimum far from it [62,78].To determine the bubble configuration it is necessary to solve the Euclidean equation of motion for Φ that minimizes the Euclidean bounce action S[Φ] of the system: 16T M 4π dr r2 1 2 where we have defined the dimensionless quantities Φ ≡ Φ/ ⟨Φ⟩ b , r ≡ M r, and V (Φ) ≡ V (⟨Φ⟩ b Φ)/(M ⟨Φ⟩ b ) 2 , and we have assumed spherical symmetry so that d 3 x = 4πr 2 dr and Φ(x) = Φ(r).The solutions to Eq. (B15) are sometimes called critical bubbles, since these are the field configurations that will not collapse upon formation due to pressure. 17e have now achieved the main purpose of all of the manipulations and definitions of this section: to show that the equation of motion for Φ depends only on λ [37,83].Indeed We find Φ bub (r) by numerically solving Eq. (B15) for various values of λ ∈ (0, 9/8), and subsequently compute the corresponding action.Note that S, being an integral over infinite space, is potentially infinite.However, only action differences, taken relative to the initial metastable false vacuum value ⟨Φ⟩ fv , ever enter the physical quantities in b /(T M ) in Eq. (B14) depends on just this parameter combination, the rest of it being a function of λ.The thin grey vertical line separates the cPT (left half) and hPT (right half) regimes.The red dashed lines are the semi-analytic expressions in Eq. (B23), while the blue dashed lines show the estimates from Eq. (B25).Note that the linear scale of the abscissa for λ > 1 has been increased for easier reading.
Our results are shown in Figs. 8 and 9.In Fig. 8 we plot the bubble configurations Φ(r) for representative values of λ; the solutions corresponding to cPTs and λ ∈ (0, 1) are shown in the left plot, and those corresponding to hPTs and λ ∈ (1, 9/8) are shown in the right plot.In Fig. 9 we plot the corresponding values of the bubble action (difference) S as a function of λ.

Daisy resummation
Finally, we comment about the cubic term of Eq. (B1) defined in Eqs.(B5) and (B6).This term is of vital importance to the order of the phase transition; were it to go away, the phase transition would cease to be first order."Daisy" resummation of "ring" diagrams effectively leads to a temperature-dependent correction to the bosonic mass stemming from its self-energy [64][65][66][67][68].This means that the cubic term heuristically becomes where g B T accounts for the model-dependent thermal mass of the B boson.From now on we will simply take g B ∼ y B , a conservative choice, since the thermal mass contributions are typically further suppressed by factors of O(10); see, for example, Ref. [64].
Evidently, for very large temperatures the Φ 3 dependence vanishes.In order to ensure that the hPT is indeed strongly first order, we exclude from our results the region of parameter space that yields sizable daisy contributions.
To that purpose we estimate the relative size of the self-energy corrections to be daisies and then conservatively demand this to be no larger than 50%.This condition is satisfied for Note that Eq. (B19) is function of λ and thus of the temperature: larger temperatures yield larger daisy contributions, and thus more values of ∆ violate Eq. (B20) and are thus excluded from our analysis.Since cooling and heating  B20) is not satisfied, for different values of λ/µ 2 .In these regions the daisy resummation corrections to the potential cubic term become relevant and the first-order phase transition ceases to be strong.We evaluate ∆ daisies at λ = 1 to obtain the representative constraint Eq. (B21), which we use when presenting the sensitivity of GW detectors in our results (e.g., Fig. 6).In green we show the region for which Eq. ( C7) is satisfied and the bubbles are in the runaway regime.The dotted black lines correspond to contours of the runaway efficiency κΦ, see Eq. (C7).Note that the linear scale of the abscissa for λ > 1 has been increased for easier reading.
PTs are sensitive to different temperatures, the restrictions on the parameter ∆ are different for each.However, since our 50% threshold for the thermal corrections is an arbitrary number, and since we are working with approximate expressions (we have taken g B ∼ y B ), we simply evaluate this condition at T c (i.e., λ = 1).This is morally equivalent to the condition ⟨Φ⟩ b /T c ≫ 1 used in other parts of the literature (e.g., Ref. [79,104]).The final constraint is: In Fig. 10 we show the regions of λ-∆ parameter space for which the daisy contributions satisfy Eq. (B20).

Analytic expressions for the bounce action
It turns out that Eq. (B15) can be solved analytically for the special cases where λ ≈ 1 (T ≈ T c ) and λ ≈ 0 (T ≈ T 0 ) or λ ≈ 9/8 (T ≈ T 1 ).In the first case the energy difference between the stable and metastable minima of the potential is very small and the thin-wall approximation can be employed.On the other hand, in the second case the stable minimum of the potential is very deep and the thick-wall approximation can be used instead.Following Ref. [62] for both cases, we find that the integral in Eq. (B14) for the bounce action scales like 18 (B22) 18 For the detailed calculation see our explanatory notebook 02 phase transition.nb.
Using Eq. (B16) and the power-law behaviors listed above, we are able to find a reasonable semianalytic fit to S(λ) for the cPT and hPT regimes (λ ∈ (0, 1) and λ ∈ (1, 9/8), respectively), which we show as dashed red lines in Fig. 9: These equations, although analytic, are still too complicated to be useful for back-of-the-envelope estimates.We can do better.
While λ is a compact variable, S spans many orders of magnitude.This means we can find a linear fit of λ to ln S to obtain an even simpler approximation that can be inverted (in order to find λ(S) analytically).For this fit it is best to use the slope d ln S/dλ at the inflection point of ln S, since that is where the slope is at its smallest, and thus a linear fit will have the broadest range of applicability.The inflection points and corresponding slopes for the cPT and hPT ranges are (B24) which lead to the following linear fits for S(λ) and their inversions: cPT regime;   temperature scales we consider in this paper.The typical value of λ at the nucleation time tn can then be read at the intersection of these lines.Note that the linear scale of the abscissa for λ > 1 has been increased for easier reading.
The bubble nucleation rate per unit volume is given by [37,62,78 where S is the finite-temperature bubble bounce action difference S[Φ bub ] − S[⟨Φ⟩ fv ] found in Eq. (B14), and M is a quantity with the dimensions of energy.Its precise form is model dependent and can only be computed numerically.Dimensional analysis, however, tells us that M ∼ T , and we simply take the remaining dimensionless coefficient to be O(1).Finally, the S 3/2 prefactor comes from the three translational zero modes that appear in the path-integral computation of the partition function, and which must be treated separately [37,62].In Fig. 12 we plot Γ(λ)/(VT 4 ) for various values of µ √ ∆/λ which, as we saw above, is the parameter combination controlling the size of the action S. We can see that the nucleation rate vanishes at λ = 1 and λ = 0 (λ = 9/8), corresponding to T = T c and T = T 0 (T = T 1 ) respectively.As long as µ √ ∆/λ is not too large, a significant number of bubbles can be nucleated within a Hubble space-time patch.Note that during reheating T 1 (and thus λ = 9/8) may never be reached, in which case Γ/(VT 4 ) will reach a maximum value at λ(T max ) and then go back down to 0 as T decreases back to T c .

Bubble expansion
In this section we fill in the details of the bubble runaway condition described in Sec.III of the main text.Recall that to test if a runaway occurs, one needs to consider the net pressure P tot acting on the wall of a relativistic bubble, which is a sum of pressure differences (between the inside and the outside of the bubble) due to the zero-temperature vacuum energy and the plasma interactions with Φ (equal to the leading mean-field contribution to the thermal potential) [84]: where in the last equality we have used Eq.(B6) and the fact that and where ⟨Φ⟩ fv and ⟨Φ⟩ tv denote the false and true vacuum Higgs vacuum expectation values, outside and inside the bubble respectively.Note that, of course, sign 1 − λ used here and sign (T c − T ) used in the main text, are equivalent.
Putting this together allows us to write the runaway condition as simply where in the last line we defined the runaway efficiency κ Φ , which can be rewritten as κ Φ = sign(1 − λ)(α − α ∞ )/α; see Eq. ( 21) in the main = text.In Fig. 10 we show the regions of λ-∆ parameter space for which the runaway condition in Eq. (C7) is satisfied.It can be seen that, up to the different values of λ explored by the transition, the runaway conditions for cPTs and hPTs are almost mirror images of each other, as described previously in the text.Indeed, since the (anti)friction is ultimately described in terms of the interaction strength y i of the particles with the Higgs, we see that runaway for cPTs (hPTs) requires sufficiently small (large) couplings.
Finally, a word about higher-order terms.Relatively recent developments attempt to compute the next-to-leadingorder contributions to the friction coming from transition radiation of gauge bosons, finding that these scale proportionally to the Lorentz factor ∼ γ [86,[105][106][107]. 19This could severely curtail the chances of a cPT entering a runaway regime.Whether such higher-order terms are present in hPTs, where the particles lose their mass upon entering the bubbles of the true vacuum, or whether they behave as a friction or antifriction, is a question left to future work.In any case, if the DS plasma does not contain such massive gauge bosons this new friction term could very well not be there [5].

Percolation, bubble number density, and mean bubble separation
The bubble nucleation rate Γ/V and the bubble-wall speed v w can be used to determine the phase in which the Universe finds itself at any given time.Indeed, the metastable volume fraction h(t) is given by [37,63,76,77] where t c is the time at which T (t c ) = T c , a(t) is the scale factor at time t, and we have assumed v w to be constant.
The PT completes at the percolation time t PT when the metastable volume fraction h(t PT ) has been reduced to 1/e.Since the bubbles can only nucleate with a probability Γ/V(t) in the metastable phase volume, quantified by h(t), we can write the bubble number density and corresponding mean bubble separation R b (t) as a function of time: We can then derive the mean bubble separation scale at percolation, R PT ≡ R b (t PT ), and from there β = (8π PT .As an interesting aside, we point out that one can use R b (t) to create an alternative definition of the percolation time.Indeed, the average bubble radius is One can then define t PT to be the time at which ⟨R(t PT )⟩ = R b (t PT ).This seems like a natural way to determine when bubble collisions take place.We have numerically corroborated that this alternative definition yields similar results to the standard one of h(t PT ) = 1/e, and thus we use the latter, more common one in our computations.Finally, we have also numerically made sure that the impact of the expansion of the Universe is negligible, and we therefore use a(t) = 1 in Eqs.(C8)-(C10) throughout this paper.

Analytic expressions for the phase transition
As discussed in the main text, PTs during reheating fall into two main categories.In cPTs and most hPTs, the transition occurs via exponential nucleation, where most bubbles are produced close to the end of the PT.In hPTs for which either T 1 > T max or T 1 simply does not exist (∆ > 8/9), the PT can take place via simultaneous nucleation, as long as t n ∼ t max or Γ/V < ∼ H 4 for all times.In this regime most of the bubbles of the new phase are produced at the time t max , where the action S reaches its minimum and the nucleation rate Γ/V reaches its maximum (since T /T c = T max /T c and therefore λ is at its largest).Ignoring the expansion of the Universe by setting a(t) = 1 in Eqs.(C8) and (C9), we can find approximate expressions for h(t) and n b (t) in each of these regimes [81,82].
Starting with the exponential nucleation regime, we can expand S(t) around t PT in a Taylor series.To first order, with S PT ≡ S(t PT ), S 1 ≡ S ′ (t PT ) < 0. Plugging this into Γ/V, we can then use the method of steepest descent to find analytic expressions for h(t) and for n b (t PT ), which we approximate as the asymptotic value of n b (t → ∞) since n b (t) changes little after that time: In the exponential nucleation regime, the condition h(t PT ) = 1/e marking the end of the PT is equivalent to 8πv  For those hPTs in the simultaneous nucleation regime, the nucleation rate Γ/V is dominated by the minimum of the action, which takes place at the temperature the farthest from T c , namely T max .Since T ′ (t max ) = 0 by definition, the expansion of S(t) around t max looks like We can then use Eq.(B7) and Eq.(A23) to find λ max ≡ λ(T max ) and from there S min = S(λ max ), Γ(t max )/V, and dS/d ln T | tmax (using Eqs.(B25) and (B26)).An even simpler estimate relies on noting that for simultaneous nucleation to take place Γ(t max )/V cannot be too different from H(t max ) 4 (i.e., t n ≈ t max ); otherwise, Γ/V would be rising exponentially and the system would be in the exponential nucleation regime.Since H max ≈ H i we simply have Γ(t max )/V ∼ H 4 i and S min ∼ 4 ln(T max /H i ).It can be shown that d 2 ln T /dt 2 ∼ −H 2 i max[1, γ] at t max , and therefore As a final word on the topic, we would like to mention that a recent study [81] showed that the GW spectrum coming from PTs is essentially the same regardless of whether the PT took place in a simultaneous or exponential regime, although others [91,108]  ).This is roughly accurate to O(10%). Figure 13: SNR contours for 1-year observation time of the SGWB spectra generated during reheating by hPT bubble collisions (red) and cPT sound waves (blue), for the upcoming BBO detector [6][7][8], as a function of Tc/Tmax and ∆, with γ = 0.1 (left) and γ = 10 (right).The corresponding values of Hi and A are also shown.The notation is the same as in Fig. 6, namely, the points inside the green contour have double peaks (from both the hGWs and the cGWs), the region above the dashed line has Tmax < T1, the region above the dotted line has Γ/V < H 4 at tmax, the dark grey region has no runaway hPT bubbles, and the light grey region has daisy contributions to the thermal potential that prevent a SFOPT.In the left panel, γ = 0.1 and we have added a purple line separating the parameter regions in which the cPT takes place during RD and during χD; for γ = 10 in the right panel and γ = 1 in Fig. 6, the cPT always takes place during RD.For this plot we chose Tc = 1 TeV, g * = 10, {µ, λ} = {1, 1}, and vw,cPT = 0.05. Figure 14: SNR contours for 1-year observation time of the SGWB spectra generated during reheating by hPT bubble collisions (red) and cPT sound waves (blue), for the upcoming BBO detector [6][7][8], as a function of Tc/Tmax and ∆, with Tc = 100 GeV and vw,cPT = 0.03 (left), and Tc = 10 TeV and vw,cPT = 0.12 (right).Note that we tune vw,cPT so that the cGWs and hGWs have similar amplitudes.The corresponding values of Hi and A are also shown.The notation is the same as in Fig. 6, namely: the points inside the green contour have double peaks (from both the hGWs and the cGWs), the region above the dashed line has Tmax < T1, the region above the dotted line has Γ/V < H 4 at tmax, the dark grey region has no runaway hPT bubbles, and the light grey region has daisy contributions to the thermal potential that prevent a SFOPT.For this plot we chose γ = 1, g * = 10, and {µ, λ} = {1, 1}.

Figure 4 :
Figure4: Dynamics of the bubble wall in its rest frame.The plasma particles of the incoming flux (green arrows to the right of the bubble wall) move with a velocity of −vw.Top: the case of a cooling phase transition, where the plasma particles gain mass upon entering the bubble (blue region), thereby slowing down (short green arrows to the left of the bubble wall).Momentum conservation means the plasma exerts a friction force on the wall (blue arrows), which moves it to the left; in the plasma rest frame the bubble wall is slowed down.Bottom: the case of a heating phase transition, where the particles lose mass inside the bubbles (red region) and thus accelerate (long green arrows), thereby accelerating the bubble wall (red arrows).

Figure 5 :
Figure5: Stochastic GW background from the phase transitions occurring during reheating.The red curve represents the GW spectrum arising from bubble collisions during the hPT (hGW) and the blue curve represents the GW spectrum originating from sound waves in the cPT (cGW).Both spectra are detectable by the future GW probe BBO[6][7][8], as shown by the PIS curves[102] (black lines, solid for GWs from bubble collisions, and dashed for GWs from sound waves).For this plot we chose the parameters corresponding to the starred benchmark point in Fig.6, namely Tc = 1 TeV, Hi = 2 × 10 −15 TeV = Γχ (i.e., Tc = 0.8 Tmax), and g * = 10, the potential coefficients {µ, A, λ} = {1, 0.72, 1} (∆ = 0.7), and a bubble-wall speed in a cPT of vw,cPT = 0.05.

Figure 9 :
Figure 9: The action S = S[Φ bub ] − S[⟨Φ⟩ fv ] as a function of λ.Note that it has been rescaled by µ √ ∆/λ = 2A/ √ 3λ 3 , since ⟨Φ⟩ 2b /(T M ) in Eq. (B14) depends on just this parameter combination, the rest of it being a function of λ.The thin grey vertical line separates the cPT (left half) and hPT (right half) regimes.The red dashed lines are the semi-analytic expressions in Eq. (B23), while the blue dashed lines show the estimates from Eq. (B25).Note that the linear scale of the abscissa for λ > 1 has been increased for easier reading.

Figure 10 :
Figure10: λ-∆ parameter space.Shaded in grey are the regions for which Eq. (B20) is not satisfied, for different values of λ/µ 2 .In these regions the daisy resummation corrections to the potential cubic term become relevant and the first-order phase transition ceases to be strong.We evaluate ∆ daisies at λ = 1 to obtain the representative constraint Eq. (B21), which we use when presenting the sensitivity of GW detectors in our results (e.g., Fig.6).In green we show the region for which Eq. (C7) is satisfied and the bubbles are in the runaway regime.The dotted black lines correspond to contours of the runaway efficiency κΦ, see Eq. (C7).Note that the linear scale of the abscissa for λ > 1 has been increased for easier reading.
to the dashed blue lines in Fig.9.As an illustrative benchmark, for T ∼ 1 TeV and H ∼ T 2 /m Pl ∼ 10 −15 TeV, S ∼ 4 ln(T /H) ∼ 138 and λ ∼ 0.74 for a cPT and λ ∼ 1.09 for an hPT, taking µ √ ∆/λ = 1.Note that λ in Eq. (B25) depends logarithmically on S, so varying H and T has little impact on the resulting λ.From the above approximations and Eq.(B7) we can obtain the rate of change of S with respect to the temperature: of |d ln S/d ln T | as a function of ∆ and λ, for both the cPT and hPT regimes, are shown in Fig.11, as well as their corresponding estimates from Eq. (B26).

Figure 11 :
Figure 11: Plots of |d ln S/d ln T | as a function of λ and ∆, for both the cPT (blue) and the hPT (red) regimes.Also shown are the curves corresponding to the approximate expressions listed in Eq. (B26).Note that the linear scale of the abscissa for λ > 1 has been increased for easier reading.

Figure 12 :
Figure12: Bubble nucleation rate (Γ/V)/T 4 as a function of λ, for different values of µ √ ∆/λ.The blue (red) lines correspond to the nucleation rate during a cPT (hPT).The dashed black line shows (10 −15 ) 4 , a typical value of H 4 for the O(TeV) temperature scales we consider in this paper.The typical value of λ at the nucleation time tn can then be read at the intersection of these lines.Note that the linear scale of the abscissa for λ > 1 has been increased for easier reading.

Figure 15 :
Figure15: SNR contours for 1-year observation time of the SGWB spectra generated during reheating by hPT bubble collisions (red) and by cPT plasma sound waves (blue), for the upcoming BBO detector[6][7][8], as a function of Tc/Tmax and γ, with ∆ = 0.5 (left), and ∆ = 0.7 (right).The corresponding values of Hi are also shown.The regions between the green contour near the dotted line and above the green contour at the left-upper corner have double peaks (from both the hGWs and the cGWs), as in Fig.6.The region to the right of the dotted line has Γ/V < H 4 at tmax.For this plot we chose Tc = 1 TeV, g * = 10, {µ, λ} = {1, 1}, and vw,cPT = 0.05.
Using the chain rule for S ′ and Eqs.(A34) and (B26) we can estimate β.Taking advantage of the fact that in the exponential nucleation regime t c < ∼ t n < ∼ t PT , we note that S PT < ∼ S n , ρ χ /ρ r,hPT < ∼ ρ χ,i /ρ r,c , H PT < ∼ H n < ∼ H c < H i , and that λ(S PT ) ≈ λ(S n ) ≈ λ infl.since λ depends only logarithmically on S for a large range of values (see Eqs. (t) ≈ S min + 1 2 S 2 2 (t − t max ) 2 , (C18)with S min ≡ S(t max ) and S 2 2 ≡ S ′′ (t max ).Repeating the exercise, we find