We study the generation of gravitational waves GWs during a firstorder cosmological phase transition PT using the recently introduced Higgsless approach to numerically evaluate the fluid motion induced by the PT. We present for the first time spectra from strong firstorder PTs alpha  0.5 alongside weak alpha  0.0046 and intermediate alpha  0.05 transitions previously considered in the literature. We test the regime of applicability of the stationary source assumption characteristic of the soundshell model and show that it agrees with our numerical results when the kinetic energy sourcing GWs does not decay with time. However we find in general that for intermediate and strong PTs the kinetic energy in our simulations decays following a power law in time and provide a theoretical framework that extends the stationary assumption to one that allows to include the time evolution of the source. This decay of the kinetic energy potentially determined by nonlinear dynamics and hence related to the production of vorticity modifies the usually assumed linear growth with the source duration to an integral over time of the kinetic energy fraction effectively reducing the growth rate. We validate the novel theoretical model with the results of our simulations covering a broad range of wall velocities. We provide templates for the GW amplitude and spectral shape for a broad range of PT parameters.
['General Relativity and Quantum Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
Gravitational waves from decaying sources in strong phase transitions
{'Gravitational waves from decaying sources in strong phase transitions': "Chiara Caprini, a,b Ryusuke Jinno, c Thomas Konstandin, d Alberto Roper Pol, a, 1 Henrique Rubira, e Isak Stomberg d, 2  \n- a D'epartement de Physique Th'eorique, Universit'e de Gen'eve, CH-1211 Gen'eve, Switzerland\n- b Theoretical Physics Department, CERN, CH-1211 Gen'eve, Switzerland\n- c Department of Physics, Kobe University, Kobe 657-8501, Japan\n- d Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany\n- e Physik Department T31, Technische Universitat Munchen James-Franck-Straße 1, D-85748 Garching, Germany  \nE-mail: chiara.caprini@cern.ch, jinno@phys.sci.kobe-u.ac.jp, thomas.konstandin@desy.de, alberto.roperpol@unige.ch, henrique.rubira@tum.de,  \nisak.stomberg@desy.de  \nAbstract. We study the generation of gravitational waves (GWs) during a first-order cosmological phase transition (PT) using the recently introduced Higgsless approach to numerically evaluate the fluid motion induced by the PT. We present for the first time spectra from strong first-order PTs ( α = 0 . 5), alongside weak ( α = 0 . 0046) and intermediate ( α = 0 . 05) transitions previously considered in the literature. We test the regime of applicability of the stationary source assumption, characteristic of the sound-shell model, and show that it agrees with our numerical results when the kinetic energy, sourcing GWs, does not decay with time. However, we find in general that for intermediate and strong PTs, the kinetic energy in our simulations decays following a power law in time, and provide a theoretical framework that extends the stationary assumption to one that allows to include the time evolution of the source. This decay of the kinetic energy, potentially determined by non-linear dynamics and hence, related to the production of vorticity, modifies the usually assumed linear growth with the source duration to an integral over time of the kinetic energy fraction, effectively reducing the growth rate. We validate the novel theoretical model with the results of our simulations covering a broad range of wall velocities. We provide templates for the GW amplitude and spectral shape for a broad range of PT parameters.  \n1", '1 Introduction': 'The door to a new era with the promise of groundbreaking discovery was opened with the inaugural direct detections by the LIGO-Virgo collaboration of gravitational waves (GWs) emanating from mergers of black holes and neutron stars [1-3]. The forthcoming observing run by the LIGO-Virgo-KAGRA (LVK) collaboration are expected to accumulate more events [4]. Efforts among Pulsar Timing Array (PTA) collaborations have furthermore unveiled convincing evidence of a stochastic gravitational wave (SGWB) background at nanoHertz frequencies [5-8]. While a compelling candidate for the source of this radiation is the superposition of supermassive black hole mergers, implying an astrophysical origin, it is important to point out that primordial sources of cosmological origin can also explain the observed signal (see, e.g., Refs. [9, 10]). These breakthroughs in GW detection give us ears to astrophysical events and cosmological history inaccessible through other means of observation.  \nLooking into the 2030s, the launch of the Laser Interferometer Space Antenna (LISA) mission [11, 12], designed to probe GWs in the unexplored milli-Hertz frequency band, is setting the stage for a potential overhaul of modern cosmology [13-15]. Until the time of launch, joint strides in data analysis techniques and theoretical progress are necessary to leverage the full potential of the LISA mission once it flies. It is in this spirit that studies on cosmological SGWB production gain their motivation.  \nCosmological sources of GWs that can be explored by LISA include inflation, particle production, topological defects, and primordial black holes throughout different stages of the expansion history of the Universe (see Refs. [16, 17] and references therein). Of particular interest is the phenomenon of first-order phase transitions (PTs) [18] that could have occurred in the early Universe. While the Standard Model (SM) predicts a crossover at the electroweak scale [19-21], many theories beyond the Standard Model accommodate a first-order electroweak PT (see Ref. [22] and references therein). During such a transition, the order parameter initially becomes trapped in a false vacuum expectation value in the symmetric phase. Subsequently, vacuum or thermal fluctuations locally induce a transition to the true vacuum in the broken phase, forming tiny seeds of bubbles [23-25]. The released vacuum energy drives the expansion of these bubbles, which eventually collide with each other, generating anisotropic stresses in the energy distribution and thus sourcing GWs. This process is highly non-thermal, suggesting that the baryon asymmetry of the Universe might have its origin in it [26-30]. While bubble collisions themselves are an important source of GWs [31-39], it has been shown in Ref. [40] that compressional fluid motion induced in the primordial plasma by the scalar walls often dominates the GW production for PTs when the broken-phase bubbles do not run away, 1 which is expected to be the case unless the PT is dominated by vacuum, e.g., for supercooled PTs [22, 45-47]. This occurs when the friction exerted on the bubble walls by the fluid particles is strong enough to balance the vacuum energy released and, hence, the bubbles reach a terminal velocity [36, 40, 46, 48, 49].  \nThe production of GWs from fluid perturbations can be decomposed in two contributions: sound waves (or acoustic/compressional turbulence) [40, 50-61] and vortical turbulence [32, 53, 55, 62-73]. While analytical modeling of GW production from fluid perturbations is important [38, 51, 53, 54, 59, 64, 71, 74], numerical simulations are essential for a comprehensive understanding of the entire process. It is believed that after a first-order PT, the fluid motion initially manifests as compressional motion (well approximated by sound waves when the fluid perturbations can be assumed to be of linear order) and then develops non-linearly, allowing for the formation of shocks and vorticity, and the subsequent development of turbulence [75, 76]. The non-linear evolution is inevitable due to the large Reynolds number of the fluid in the early Universe [44]. The transition from the fully compressional to the vortical turbulence regime can be especially important in strong transitions, when non-linearities can play an important role.  \nCurrently, large-scale simulations have been performed by the Helsinki-Sussex group, which numerically solve a coupled scalar field-fluid system [40, 50, 52, 77]. First-order PTs in the early Universe exhibit a significant hierarchy: this is the hierarchy between the typical scales inherent in the order parameter field (including the thickness of the walls) and those in the cosmological fluid (including the bubble size and the sound-shell thickness). For the  \nelectroweak PT, the hierarchical separation can be as large as M P /T ∼ 10 16 . This fact naturally leads to the idea of the Higgsless scheme proposed by part of the authors of the present paper [56, 58]. In this scheme, the microphysics of the wall is introduced as a nondynamical (although space and time dependent) energy-injecting boundary condition within the bag equation of state [49, 78, 79], and the bubble walls are assumed to have reached a terminal velocity, such that the fluid perturbations reach a self-similar solution in a very short time scale (much shorter than the time scale for collisions) [49]. See also Refs. [46, 80, 81] for discussions on the bubble wall terminal velocity. As a result, the scheme is able to capture the macroscopic dynamics necessary for GW production without being required to also solve for the hierarchically smaller scales.  \nIn this paper, we explore the previously uncharted realm of GWs sourced by fluid perturbations induced in strong first-order PTs. 2 We also update results for weak and intermediate transitions and compare with other results in the literature [52, 54, 58]. We perform approximately 1000 simulations involving a parameter scan over wall velocities v w ∈ [0 . 32 , 0 . 8] in increments of 0.04 for weak ( α = 0 . 0046), intermediate ( α = 0 . 05), and strong ( α = 0 . 5) PTs, using the Higgsless approach, to clarify the dependence on the underlying physical quantities of various characteristics of the GW spectrum, especially focusing on its overall amplitude. We then assess the long-term evolution of the system, and discuss indications of non-linearities and their impact on the GW production. We also include an analysis of the longitudinal and transverse components of the velocity field spectra.  \nThe outline of the paper is as follows. In Sec. 2, we review the GW production in a PT. We first describe the bubble nucleation history and the conservation laws of the fluid perturbations in Secs. 2.1 and 2.2. Then, Sec. 2.3 presents the description used to compute the GWspectrum at present time from the numerical results, more detail is given in App. A. We review in Sec. 2.4 the GW production under the stationary assumption, found, for example, in the sound shell model (which assumes linear fluid perturbations) [51, 54, 59, 60], and used to interpret numerical results [40, 50, 52, 56, 58]. We then present a novel model in Sec. 2.5 that extends the unequal time correlator (UETC) to a locally stationary UETC, allowing us to introduce the effect of the source decay. In Sec. 2.6, we extend this model to estimate the effects of the expansion of the Universe. In Sec. 3, we discuss updates to the previous Higgsless simulations [58], and summarize the physical and numerical setup of the simulations. In Sec. 4, we discuss the numerical results. We first present a convergence analysis for the kinetic energy and the integrated GW amplitude in Sec. 4.1, which is combined with a convergence analysis that includes the potential effects of underresolving the fluid perturbations, presented in App. B. The result of these convergence analyses is an estimate of the kinetic energy fraction when the entire volume is converted to the broken phase, which is then compared to the single-bubble value commonly used in GW studies. We then analyse in Sec. 4.2 the time evolution of the kinetic energy fraction and present a decaying power-law in time fit that accurately reproduces the simulations. To investigate the origin of the observed decay, we briefly discuss the development of vorticity in App. C, but we leave a detailed treatment of the onset of vorticity for future work. The growth of the GW integrated amplitude with the source duration is presented in Sec. 4.3, together with estimates of the GW efficiency ˜ Ω GW . We investigate the GW spectral shape in Sec. 4.4 and fit the numerical results to a double broken power law, providing estimates for the positions of the relevant spectral scales. We  \nsummarize our results in Sec. 5, and present a template for the GW spectrum.', '2 Gravitational waves from a phase transition': 'In this section we describe the production of gravitational waves from the fluid perturbations induced in the primordial plasma by the nucleation of bubbles in a first-order PT. In Sec. 2.1, we review the nucleation history of broken-phase bubbles that is used in the Higgsless approach, and in Sec. 2.2, we review the the relativistic hydrodynamic equations. In Sec. 2.3, we describe how the GW generation is tackled within the Higgsless simulations, and discuss its applicability (see also App A for details). In Sec. 2.4, we review the production of GWs under the assumption that the anisotropic stress UETC is stationary. This is usually assumed in the literature for sound-wave sourcing of GWs [22, 40, 50-52, 54, 59]. In Sec. 2.5, we propose an extension of the stationary model to a locally stationary UETC, allowing to account for the dynamics of a source that is decaying in time. We also provide a proxy to extend our results to the expanding Universe in Sec. 2.6. We validate this novel proposed model in Sec. 4 with the results from our numerical simulations, described in Sec. 3.', '2.1 Bubble nucleation histories': 'In the Higgsless approach [58], a fundamental assumption is that the broken-phase bubbles reach a terminal wall velocity v w due to the friction exerted by the fluid particles [46, 49, 80, 81], such that v w can be prescribed as an input of the simulations. This enables the construction of bubble nucleation histories , encompassing nucleation times and locations, as well as the predetermined expansion of the bubbles. We assume an exponentially increasing in time probability of bubble nucleation, P ( t ) ≃ P n exp[ β ( t -t n )], where β determines the usual rate of bubble nucleation evaluated at the nucleation time t n , such that the action S has been Taylor-expanded around this time [22]. The time dependence is hereby inherited from the temperature dependence of the tunnelling action S/T and the fact that the temperature scales inversely with the scale factor  \nβ = d dt ( S T )∣ ∣ ∣ ∣ t = t n = -HT d dT ( S T )∣ ∣ ∣ ∣ T = T n , (2.1)  \nwhere H is the Hubble scale. A detailed description of how such bubble nucleation histories are constructed is found in Ref. [58], and examples of how modified bubble nucleation histories can be constructed are found in Refs. [83, 84]. We note that in previous numerical work [50, 52, 77], bubbles are nucleated simultaneously. This is expected to have an impact on the spectral peak and the amplitude, but not on the spectral shape [52, 85].', '2.2 Relativistic hydrodynamic equations': 'The relativistic hydrodynamic equations of motion are derived from the conservation of the energy-momentum tensor T µν , which, in Minkowski space-time reads  \n∂ µ T µν = 0 . (2.2)  \nThese equations of motion hold during the PT under the assumption that the duration of the transition is much shorter than a Hubble time, i.e., β/H ∗ ≫ 1, which allows to neglect the expansion of the Universe. They can also be applied to the fluid motion after the PT ends if the fluid is dominated by radiation particles: indeed, the conservation laws then become  \nconformally invariant and hence reduce to those in Minkowski space-time by a conformal transformation [42, 86]. We take T µν to be that of a perfect fluid  \nT µν = u µ u ν w -η µν p , (2.3)  \nwhere u µ = γ (1 , v i ) is the fluid four-velocity, γ = 1 / √ 1 -v 2 is the gamma factor, w and p are respectively the enthalpy and pressure in the system, and η µν = diag { 1 , -1 , -1 , -1 } is the Minkowski metric tensor. The equations of motion couple to the state of the vacuum through the bag equation of state, for which we take the sound speed to be c 2 s = 1 / 3,  \np = 1 3 aT 4 -ϵ , w = T dp dT = 4 3 aT 4 , (2.4)  \nwith T being the temperature. The bag constant ϵ [49], defined as the difference in vacuum energy density between the symmetric and broken phases, is thus promoted to a time- and space-dependent quantity  \nϵ ( t, x ) = { 0 inside bubbles , ϵ outside bubbles , (2.5)  \nwhose time evolution is uniquely determined for each bubble nucleation history by the terminal wall velocity v w . We therefore neglect the (model-dependent) possibility that heating in the broken phase can slow down the expansion of the Higgs front when the latter propagates as a deflagration [77], and also that strong PTs can lead to runaway behavior [46, 87, 88].  \nThe relevant quantity for boundary conditions of the fluid at the Higgs interface is the difference in the trace of the energy-momentum tensor, θ = w -4 p , normalized to the enthalpy [79],  \nα = ∆ θ 3 w = 4 ϵ 3 w , (2.6)  \nwhere the second equality holds within the bag equation of state, for which the trace anomaly reduces to ∆ θ = 4 ϵ . We use this quantity in the following to parameterize the strength of the PT in the system [22]. The conservation laws for a relativistic perfect fluid are  \n∂ t T 00 + ∇ i T i 0 = 0 , ∂ t T j 0 + ∇ i T ij ( T µ 0 , ϵ ) = 0 . (2.7)  \nNote that T ij ( T µ 0 , ϵ ) depends on the state of the vacuum such that, effectively, the expanding bubbles perturb the fluid as the latent heat of the vacuum is (locally) deposited. For more details on these equations, we refer the reader to the original Higgsless reference [58].', '2.3 Gravitational wave production': "The GW spectrum as a present-time observable is computed using the following relation  \nΩ GW ( k ) = 1 ρ tot dρ GW d ln k = 3 T GW ( H ∗ /β ) 2 I ( k ) , (2.8)  \nwhere ρ tot is the total energy density of the Universe at present time, T GW ≡ ( a ∗ /a 0 ) 4 ( H ∗ /H 0 ) 2 is the transfer function, with a value h 2 T GW ≃ 1 . 6 × 10 -5 ( g ∗ / 100) -1 / 3 [17], k is the comoving wave number, which can be converted to the observable frequency as f = k/ (2 πa 0 ), and g ∗ and g ∗ s are respectively the relativistic and entropic degrees of freedom at the time of GW production. The ratio of scale factors is a ∗ /a 0 ≃ 8 × 10 -16 (100 GeV /T ∗ ) ( g ∗ s / 100) -1 / 3 .  \nThe function I ( k ) represents the spectrum of the stochastic GW signal, and it is given by a double time integral of the anisotropic stress UETC, E Π (see Eq. (A.8) in App. A), multiplied by the Green's function of the GW equation. We have introduced the prefactor ( H ∗ /β ) 2 in Eq. (2.8) to express I in normalized time and wave number units ˜ t ≡ tβ and ˜ k ≡ k/β (see App. A, in particular Eq. (A.9), for its derivation from the solution to the GW equation)  \nI ( ˜ t ∗ , ˜ t fin , ˜ k ) = k 2 ∫ ˜ t fin ˜ t ∗ ∫ ˜ t fin ˜ t ∗ E Π ( ˜ t 1 , ˜ t 2 , ˜ k ) cos ˜ k ( ˜ t 1 -˜ t 2 ) d ˜ t 1 d ˜ t 2 . (2.9)  \nHere ˜ t ∗ and ˜ t fin denote the initial and final times of action of the GW source. In Eq. (2.9) the expansion of the Universe is not taken into account, and it therefore holds only for a short duration of the GW production process, t fin -t ∗ ≪ H -1 ∗ , where H ∗ is the Hubble rate at the PT time t ∗ . It is important to notice that, even though the duration of the PT, determined by the inverse of the nucleation rate β -1 , must be short β/H ∗ ≥ 1, this does not imply that the GW sourcing time is also short. In particular, for GWs generated by fluid motion, the typical time scale of dissipation of the fluid kinetic energy is very long, as it is set by the kinematic viscosity in the early universe [65]. We discuss in Sec. 2.6 how to extend our results, obtained with simulations performed in Minkowski space-time, to an expanding Universe.  \nThe anisotropic stress UETC E Π ( ˜ t 1 , ˜ t 2 , ˜ k ) can be evaluated numerically. In particular, within the Higgsless simulations, it is evaluated from the following expression (see Eq. (A.13) and Refs. [56, 58]), given in the normalized units of the simulation ˜ t , ˜ k , and ˜ V ≡ V β 3 :  \nI sim ( ˜ t ∗ , ˜ t fin , ˜ k ) = ˜ k 3 4 π 2 ˜ V ∫ Ω ˜ k dΩ ˜ k 4 π Λ ijlm ( ˆ k ) [ ˜ T ij ( ˜ t ∗ , ˜ t fin , ˜ q, ˜ k ) ˜ T ∗ lm ( ˜ t ∗ , ˜ t fin , ˜ q, ˜ k ) ] ˜ q = ˜ k , (2.10)  \nbeing Λ ijlm the transverse and traceless operator  \nΛ ijlm ( ˆ k ) = P il P jm -1 2 P ij P lm , with P ij = δ ij -ˆ k i ˆ k j . (2.11)  \nThe function ˜ T ij ( ˜ t ∗ , ˜ t fin , ˜ q, ˜ k ) is computed from the normalized stress-energy tensor ˜ T ij ( t, x ), sourcing the GWs, as [see Eq. (A.12)]  \n˜ T ij ( ˜ t ∗ , ˜ t fin , ˜ q, ˜ k ) = ∫ ˜ t fin ˜ t ∗ d ˜ t e i ˜ q ˜ t ∫ d 3 ˜ x e -i ˜ k · ˜ x ˜ T ij ( ˜ t, ˜ x ) , (2.12)  \nwhere ˜ T ij ( t, x ) = wγ 2 v i v j / ¯ ρ , and ¯ ρ = 3 4 ¯ w + ϵ = 3 4 (1+ α ) ¯ w is the average total energy density [see Eqs. (2.4) and (2.6)].  \nFrom Eq. (2.10) we can directly obtain the GW spectrum at present time [see Eq. (2.8)] using the normalized source ˜ T ij ( ˜ t init , ˜ t end , ˜ q, ˜ k ) computed in the simulations in an interval of time ˜ t ∈ ( ˜ t init , ˜ t end ), provided that the source has stopped operating by the end of the simulation ˜ t end . This does not necessarily imply that one needs to compute the GW spectrum until ˜ t fin in a simulation, as this might not be computationally affordable, especially for slowly decaying sources. We can assume that the source has stopped operating at a wave number k also if the amplitude of the GW spectrum has reached its saturation amplitude, entering its free-propagation regime. Therefore, the GW spectrum evaluated from the simulation would be accurate for any wave numbers k that have already reached their saturated amplitudes by the end of the simulation at ˜ t end , even if ˜ t end < ˜ t fin . However, as the modes in our simulations  \nhave not reached this regime, as in previous numerical work [40, 50, 52, 56, 58, 77], one needs to take into account this limitation when interpreting the numerical results, as we discuss in Sec. 4. For this purpose, we will use the stationary UETC usually assumed for sound waves (see Sec. 2.4) and its extension to a locally stationary UETC that can incorporate the decay of the source with time (see Sec. 2.5). These models will allow us to extrapolate our results from the final time of the simulations till the final time of GW production t fin to estimate the present-day GW spectrum. We also note that numerical viscosity, which in the simulations can dissipate kinetic energy over time could limit the time and k range for which the simulations can actually be accurate. This can potentially affect the transfer of energy from the fluid perturbations to GWs, especially at large frequencies.", '2.4 Gravitational waves from stationary sound waves': "In the following, we describe the theoretical understanding of the GW production from sound waves from previous work, which will be useful to describe our numerical results in terms of the quantities previously used in the literature. 3 In the context of GW production from sound waves, it has extensively been assumed that the UETC is stationary, i.e., it only depends on the difference t -= t 2 -t 1 , E Π ( t 1 , t 2 , k ) = 2 k 2 K 2 f ( t -, k ) [40, 50-54, 59, 62, 64, 89, 90], where K = ρ kin / ¯ ρ is a time-independent kinetic energy fraction with ρ kin = ⟨ wγ 2 v 2 ⟩ . Under this assumption, Eq. (2.9) becomes  \nI ( ˜ t ∗ , ˜ t fin , ˜ k ) = k 3 K 2 ∫ ˜ t fin ˜ t ∗ d ˜ t ∫ ˜ t fin -˜ t ˜ t ∗ -˜ t cos( kt -) f ( ˜ t -, ˜ k ) d ˜ t -. (2.13)  \nIn the sound-shell model of Refs. [51, 54], the limits of the integral over ˜ t -were extended to ±∞ , allowing to commute the two integrals, such that the integral over ˜ t simply becomes the source duration, ˜ τ sw = ˜ t fin -˜ t ∗ , yielding the linear growth of the GW amplitude usually assumed in the literature [22]. The duration of the sound-wave sourcing of GWs can then be taken to correspond to the time that it takes non-linearities to develop in the fluid, ˜ τ sw ≡ τ sw β ∼ ( βR ∗ ) / √ K , with βR ∗ ≡ (8 π ) 1 / 3 max( v w , c s ) being the mean separation of the bubbles at the end of the PT [22]. The aforementioned assumptions about the integration limits can only be considered for kτ sw ∼ kR ∗ / √ K ≫ 1, and when τ sw /R ∗ ∼ 1 / √ K ≫ 1, as shown in Ref. [59], where the results of the sound-shell model are extended to all wave numbers and values of R ∗ . Therefore, for small values of K , we expect this approximation to hold for all relevant wave numbers kR ∗ ≫ √ K . Finally, it can be shown that under the same assumptions, the remaining integral over ˜ t -is proportional to ( βR ∗ ) /c s (see, e.g., App. B of Ref. [59]), such that  \nI ( ˜ t ∗ , ˜ t fin , ˜ k ) = ˜ Ω GW K 2 βR ∗ ˜ τ sw S ( kR ∗ ) , (2.14)  \nwhere ˜ Ω GW corresponds to the GW production efficiency and S is a normalized spectral shape, such that ∫ dln k S ( k ) = 1. Using Eq. (2.8), the final GW spectrum can be written as [22, 40, 50-52, 54, 59, 91, 92]  \nΩ GW ( k ) = 3 T GW ˜ Ω GW K 2 H ∗ R ∗ H ∗ τ sw S ( kR ∗ ) . (2.15)  \nEquation (2.14) shows that one could divide the function I int sim , obtained in the simulations [see Eq. (2.10)] and integrated over ln k , by K 2 βR ∗ ˜ T GW to estimate ˜ Ω GW , where  \n˜ T GW = ˜ t end -˜ t init is the time interval of the simulation in which Eq. (2.10) is evaluated, as done in previous numerical studies [40, 50, 52, 54]. Alternatively, as we work with nondimensional length scales and times given in units of 1 /β , the previous work on Higgsless simulations presented the results of the GW amplitude considering the following parameterization [56, 58],  \nQ ' ( k ) = ( ¯ ρ ¯ w ) 2 4 π 2 ˜ T GW I sim ( ˜ t init , ˜ t end , ˜ k ) ≈ 9 π 2 4 ˜ T GW (1 + α ) 2 I sim ( ˜ t ∗ , ˜ t fin , ˜ k ) , (2.16)  \nwhere the prefactor (¯ ρ/ ¯ w ) 2 in Eq. (2.16) takes into account that the authors in Refs. [56, 58] used the mean enthalpy to normalize T ij in the definition of Q ' instead of the total energy density ¯ ρ , as done in our case [see Eq. (2.12)]. However, we confirm with simulations in Sec. 4 that using R ∗ in Eq. (2.14) to describe the GW amplitude, instead of using Q ' /K 2 [56, 58], allows us to find a value of ˜ Ω GW that is almost independent of v w . 4  \nThe kinetic energy fraction K is usually taken to be the one corresponding to a single bubble, which can be expressed as  \nK ξ ≡ κ ξ α 1 + α , (2.17)  \nwhere α characterizes the strength of the PT [see Eq. (2.6)], and κ ξ ≡ ρ kin /ϵ is the singlebubble efficiency factor [49]. We compare in Sec. 4 the kinetic energy fraction found in the simulations with the single-bubble result to take into account the effect of collisions and nonlinear dynamics in K and express the results in terms of K ξ . We note that in the sound-shell model, K = K ξ is time-independent but the exact value of K/K ξ ∼ O (1) depends on the PT parameters [51, 54, 59]. Similarly, Refs. [52, 77] have also reported maximum values of the kinetic energy fraction in their simulations K max different than K ξ .  \n̸  \nWe find in Sec. 4 that Eq. (2.14) holds in our simulations as long as the kinetic energy does not decay with time after the PT ends. Furthermore, we find that when the source decays, a generalization of the stationary UETC to include the decay of the source, which is described in Sec. 2.5, predicts a GW growth with the source duration that is given by the integrated K int . This model is accurately validated by the numerical results in Sec. 4.3. This allows us to still compute numerically the GW efficiency ˜ Ω GW and estimate the expected final GW amplitude at the final time of GW production, even when the GW does not grow linearly with the source duration.", '2.5 Gravitational waves from decaying sources': 'As we show in Sec. 4, when the kinetic energy starts to decay within the duration of the simulations, potentially due to the fact that the system enters the non-linear regime, we find that the GW amplitude deviates from the linear growth of Eq. (2.14), impeding its use to estimate the GW efficiency. In such cases, we propose a generalization of Eq. (2.13), by assuming a locally stationary UETC that allows us to include the time dependence of K 2 , E Π ( t 1 , t 2 , k ) = 2 k 2 K 2 ( t + ) f ( t -, k ), where t + = 1 2 ( t 1 + t 2 ). Then, under the same approximations discussed above that yield to the linear growth in the source duration τ sw (i.e., for kτ sw ≫ 1 and τ sw ≫ R ∗ ), we find that K 2 ˜ τ sw in Eq. (2.14) can be substituted by  \nK 2 int ,  \nyielding  \nI ( ˜ t ∗ , ˜ t fin , ˜ k ) = ˜ Ω GW K 2 int ( ˜ t ∗ , ˜ t fin ) ( βR ∗ ) S ( kR ∗ ) . (2.19)  \nWe note that K 2 int reduces to K 2 ( ˜ t fin -˜ t ∗ ) = K 2 ˜ τ sw when K 2 is constant and we recover Eq. (2.14) found in the stationary assumption. A similar UETC has been recently considered in Ref. [61], E Π ( t 1 , t 2 , k ) = 2 k 2 √ P v ( t 1 , k ) P v ( t 2 , k ) cos( kc s t -), where P v is the kinetic spectrum. We note that the integral over t 1 and t 2 of the latter can be reduced to the integral in Eq. (2.18) under the assumptions discussed in Sec. 2.4, such that t -terms can be usually neglected in the integral over t + due to the assumed small compact support of t -[54] (see discussion in Sec. V of Ref. [59]). Therefore, we expect both UETC to have the same impact on the integrated GW amplitude and we emphasize that when we validate the proposed model in Sec. 4.3, we are validating the overall amplitude but not necessarily the GW spectral shape.  \nAfter validating this model with the results of numerical simulations in Sec. 4.3, we can estimate the GW efficiency ˜ Ω GW even when the kinetic energy is decaying with time. In particular, we show in Sec. 4.2 that the kinetic energy evolution in the simulations can be in general fit to a decaying power law, K ( ˜ t ) = K 0 (∆ ˜ t/ ∆ ˜ t 0 ) -b , where b > 0 and K 0 are parameters to be fit using the numerical results. We note that ∆ ˜ t and ∆ ˜ t 0 are time intervals with respect to the time-coordinate origin in the simulations and, hence, we will simply use ˜ t and ˜ t 0 in the following whenever the expansion of the Universe is ignored. 5 We will take ˜ t 0 to be the time when all the simulation box is in the broken phase. If we assume that the GW production starts around this time ˜ t ∗ ≃ ˜ t 0 , we find using Eq. (2.18) that the dependence of K 2 int with the source duration, ˜ τ sw = ˜ t fin -˜ t ∗ , is  \nK 2 int = K 2 0 ˜ t 0 (1 + ˜ τ sw / ˜ t 0 ) 1 -2 b -1 1 -2 b . (2.20)  \nThis expression reduces to K 2 int → K 2 0 ˜ τ sw for any value of b when the source duration is very short ˜ τ sw / ˜ t 0 ≪ 1, while for long durations ˜ τ sw / ˜ t 0 ≫ 1, it takes the following asymptotic limits:  \nlim ˜ τ sw ≫ ˜ t 0 K 2 int = K 2 0 ˜ t 0 1 -2 b ( ˜ τ sw ˜ t 0 ) 1 -2 b , when b < 1 2 , lim ˜ τ sw ≫ ˜ t 0 K 2 int = K 2 0 ˜ t 0 2 b -1 , when b > 1 2 . (2.21)  \nHence, K 2 int grows unbounded proportional to ˜ τ 1 -2 b sw when 2 b ≤ 1, thus generalizing the linear growth obtained in the stationary assumption to any decay rate b . On the other hand, when the decay rate is larger than 0 . 5, then K int saturates to a value K 2 0 ˜ t 0 / (2 b -1).  \nK 2 int ( ˜ t ∗ , ˜ t fin ) ≡ ∫ ˜ t fin ˜ t ∗ K 2 ( ˜ t ) d ˜ t , (2.18)', '2.6 Effect of the Universe expansion': 'Based on the assumption that the expansion of the Universe can be neglected, a superposition of sound waves would emit GWs with an amplitude that increases unbounded linearly with the sourcing time τ sw (i.e., until the development of non-linearities or the kinetic energy is efficiently dissipated) if the UETC of the source is stationary and K ≪ 1 (see Eq. (2.15) and Refs. [40, 50-52, 54, 59, 62-64, 90-92]). Similarly, with the proposed locally stationary UETC, the GW amplitude would increase unbounded proportional to τ 1 -2 b sw when the kinetic energy decays with a decay rate b < 1 / 2 [see Eq. (2.21)], and it would only saturate when b > 1 / 2.  \nIn the stationary UETC case, to take into account the expansion of the Universe, the linear increase τ sw of Eq. (2.15) can be substituted by the suppression factor Υ( τ sw ) = τ sw / (1+ H ∗ τ sw ) [59, 89], where τ sw ≡ τ fin -τ ∗ now refers to an interval in conformal time, and still apply the results of our simulations to long-lasting sources. We note that it is the interval in conformal time, instead of the interval in cosmic time, that should be associated to the eddy turnover time R ∗ / √ K when evaluating the expected time to develop non-linearities, due to the conformal invariance of the fluid equations when the fluid is radiation-dominated [42, 86]. When including expansion, the results are no longer invariant under time translations, so we need to choose absolute values for conformal times. Assuming that the PT is short and occurs during radiation-domination, we can set the initial and final conformal times of GW production to be H ∗ τ ∗ = 1 and H ∗ τ fin = 1 + H ∗ τ sw , where we set a ∗ = 1, such that the conformal Hubble rate is H ∗ = H ∗ a ∗ = H ∗ .  \nIn the generalized case when we take into account the decay of K ( t ) [see Eq. (2.19)], unless K 2 decays faster than 1 /t , the integrated K 2 int would also diverge, requiring to cut off the GW growth at a final time of GW sourcing t fin = t ∗ + τ sw . Extending Eq. (2.13) to apply in an expanding Universe [59], an effective integrated K 2 that can be used in Eq. (2.19) to estimate the effect of expansion is the following  \nK 2 int , exp ≡ ( β/H ∗ ) 2 ∫ ˜ τ fin ˜ τ ∗ K 2 (˜ τ ) ˜ τ 2 d˜ τ = ∫ ˜ τ sw 0 K 2 (˜ τ ∗ + δ ˜ τ ) (1 + δ ˜ τ/ ˜ τ ∗ ) 2 d( δ ˜ τ ) , (2.22)  \n̸  \nwhere δ ˜ τ ≡ ˜ τ -˜ τ ∗ and ˜ τ ∗ = β/H ∗ . Taking into account that the power-law decay in flat space-time should be taken in conformal time K (˜ τ ) = K 0 (∆˜ τ/ ∆˜ τ 0 ) -b due to the conformal invariance of the dynamics for a radiation-dominated fluid, then we need to express the absolute times in a flat space-time as time intervals in conformal time (see footnote 5), ∆˜ τ = δ ˜ τ +∆˜ τ ∗ . If we again assume that ∆˜ τ ∗ = ∆˜ τ 0 , the resulting integral for 2 b = 1 can be expressed as  \nK 2 int , exp = K 2 0 ∆˜ τ 2 b 0 ∫ ˜ τ sw 0 ( δ ˜ τ +∆˜ τ 0 ) -2 b [ 1 + ( H ∗ /β ) δ ˜ τ ] 2 d( δ ˜ τ ) = K 2 0 Υ b ( τ sw ) ( β/H ∗ ) , (2.23)  \nF b ( H ∗ τ ) = ( ∆ τ ∗ + τ ∆ τ 0 ) 1 -2 b H ∗ ∆ τ 0 (1 -H ∗ ∆ τ ∗ ) 2 2 F 1 [ 2 , 1 -2 b, 2 -2 b, -H ∗ (∆ τ ∗ + τ ) 1 -H ∗ ∆ τ ∗ ] , (2.24)  \nwhere we have defined a suppression factor Υ b ( τ sw ) = ∆ F b / (1 -2 b ) with ∆ F b ≡ F b ( H ∗ τ sw ) -F b (0) that reduces to the one found for stationary sources when b = 0, i.e., Υ 0 ( τ sw ) ≡ Υ( τ sw ) = H ∗ τ sw / (1 + H ∗ τ sw ) [59, 89]. The function F b is the following  \nwhere 2 F 1 is the hypergeometric function. We highlight that the emergence of a hypergeometric function has no deep physical meaning, since Eq. (2.24) arises from introducing  \nthe chosen fit K ( τ ) in Eq. (2.22). The relevant physical quantity is the resulting modification Υ b with respect to Υ (i.e., with no decay of the source) obtained from the integral in Eq. (2.23). The value of ∆˜ τ 0 corresponds to the characteristic time ˜ t 0 used in the fit of K 2 in flat space-time. We note that, in principle, using ∆˜ τ ∗ = ˜ τ ∗ -˜ τ 0 +∆˜ τ 0 in the integrand of Eq. (2.23) allows to compute the GW spectrum starting at any time ˜ τ ∗ . For simplicity, we have chosen ∆˜ τ ∗ = ∆˜ τ 0 . We also find that for any values of b , the functions Υ b ( τ sw ) always reduce to the linear growth H ∗ τ sw for short source duration, H ∗ τ sw ≪ 1. Furthermore, we note that with the inclusion of the Universe expansion, the relevant integrated K 2 becomes K 2 0 Υ b ( H ∗ τ sw ) ( β/H ∗ ), where the β/H ∗ cancels with the normalization introduced in Eq. (2.8) and the resulting Υ b only depends on the source duration in units of the Hubble time, H ∗ τ sw , which becomes the relevant time scale in addition to the PT duration, β -1 . Then, the final GW spectrum becomes  \nΩ GW ( k ) = 3 T GW ˜ Ω GW K 2 0 H ∗ R ∗ Υ b ( τ sw ) S ( kR ∗ ) , (2.25)  \nwhere Υ b is obtained from the integrated K 2 int and in particular reduces to the expression determined by Eq. (2.24) when the fit K 2 = K 2 0 (∆˜ τ/ ∆˜ τ 0 ) -2 b holds at all times of GW production.  \nAs discussed above, in Sec. 4.2 we use the numerical results of the simulations to find the values of the fit parameters b and K 0 for different PTs. Then, we validate the assumption that Eq. (2.19) applies within the duration of our simulations in Sec. 4.3, and provide an estimate of the GW amplitude as a function of the source duration ˜ τ sw . We compare the resulting evolution in flat space-time (both from the analytical fit and using the numerical results) with the one obtained including the expected effect of the Universe expansion for different values of β/H ∗ . We emphasize that the suppression of the time intervals by Υ b due to the Hubble expansion works as a proxy to estimate its effect.', '3 Numerical setup': 'In this section, we focus on describing the numerical setup of the Higgsless simulations: in Sec. 3.1, we comment on the updates in the numerical scheme with respect to Ref. [58], and in Sec. 3.2, we describe the simulation suite considered for this work.', '3.1 Updates to the numerical setup': "In this section, we highlight three updates to the Higgsless simulations with respect to Ref. [58] aimed at improving (1) the time integration scheme, (2) the mapping between the discrete and the continuum momenta, and (3) the criterion for numerical stability in simulations of strong first-order PTs ( α = 0 . 5). For a complete description of the Kurganov-Tadmor (KT) numerical scheme [93] used for the Higgsless simulations, we refer to Refs. [56, 58].  \nCommencing with (1) , in practice, the integral in Eq. (2.12) must be computed numerically on the grid of space and time. For the space grid, this is accomplished through a fast Fourier transform routine [94]. For the time grid, in order to overcome the practical limitation of memory (i.e., storing a large number of 3D time slices), one needs to resort to another method. In the first iteration of the Higgsless simulation code, the discrete integral in time of Eq. (2.12) was approximated as  \n˜ T ij ( ˜ t init , ˜ t, ˜ q, ˜ k ) = ˜ t ∑ ˜ t ' = ˜ t init δ ˜ t e i ˜ q ˜ t ' ˜ T ij ( ˜ t ' , ˜ k ) , (3.1)  \ni.e., through its Riemann sum, by stacking past time slices weighted by a complex factor from ˜ t init until ˜ t ≤ ˜ t end for each time step δ ˜ t over which the GWs are sourced. In the current version, we improve upon this scheme by treating T ij as a piecewise linear function interpolating between the support points, using a similar scheme to the one proposed in Ref. [66] for solving the GW equation. Since the integrand involving an oscillating exponential as well as the linearized T ij is now analytically integrated, this modified routine allows to capture better the UV behavior at large k , alleviating the time-step δ ˜ t required to find accurate spectra in this regime (see discussion in Ref. [66]). However, no sizable discrepancies have been observed in the UV range of the GW spectra through this change for the dynamical range and choice of δ ˜ t used in our simulations.  \nContinuing with (2) , we begin by noting that the first version of the Higgsless simulations employed a sin-prescription for the mapping of discrete momenta on the grid to their correspondents in the continuum. Care must be taken that on the grid of the simulation with N points per dimension, Fourier modes with momenta -l i and N -l i (in the i th direction) are equivalent and mapped to the same momenta in the continuum. At the same time, momenta of order l i ≃ N are equivalent to l i ≃ 0 and should be considered soft. Depending on whether the observable under consideration is sensitive to the sign of the momentum, this motivated the mapping  \n˜ k i = 3 -a δ ˜ x sin ( aπl i N ) , (3.2)  \nwhere δ ˜ x = ˜ L/N , with a = 2 when the sign is relevant and a = 1 when it is not. In the current simulations, we generally use a saw description for the momenta  \n˜ k i =      2 πl i / ( Nδ ˜ x ) , l i < N/ 2 , 0 , N/ 2 , 2 π ( l i -N ) / ( Nδ ˜ x ) , l i > N/ 2 . (3.3)  \nAs such, the saw-prescription avoids different descriptions in different contexts (such as the space Fourier transforms for the GW estimate or for the numerical fluid evolution) and maintains a good map of momenta all the way to l i ≃ N/ 2, while the previous method is only accurate in the linear regime of the sine function. At the moment, we do not find substantial differences between the two implementations, but we expect this implementation to improve the results when increasing the resolution of the Higgsless simulations.  \nThe third point (3) concerns the choice of the maximal local velocity a j +1 / 2 (on a staggered cell in direction j ), appearing in Eq. (3.7) of Ref. [58]. In summary, this quantity enters the flux limiter used in the KT scheme to preserve the shock structures in the lattice by setting a minimal numerical viscosity that reduces spurious oscillations and improves the stability of the numerical scheme. In the limit of small fluid velocities, i.e., for weak and intermediate PTs, a j +1 / 2 = c s = 1 / √ 3 is a good choice. In the case of strong PTs, however, fluid velocities often supersede 1 / √ 3 and approach 1. To improve the numerical stability of the simulation, we therefore choose a j +1 / 2 = 1 for strong PTs. In the weak regime, the numerical changes due to this choice are negligible but for stronger PTs, it improves the stability of the code significantly. In rare occasions and close to shocks, the simulation can lead to unphysical fluid velocities (essentially v > 1) as a numerical artifact. In these cases, we opted to enforce the local fluid velocity to 1. This only happened in isolated points and had no measurable impact on the conservation of T 0 µ or the GW spectra.  \nIn all other regards, the current version of the Higgsless implementation is identical to the first version in Ref. [58].", '3.2 Simulations and parameter choices': "We list the parameters considered in this study in Tab. 1. We expand upon Ref. [58] by including in our parameter scan strong PTs with α = 0 . 5. We thus run reference simulations for α ∈ { 0 . 0046 , 0 . 05 , 0 . 5 } and wall velocities v w ∈ { 0 . 32 , 0 . 36 , ..., 0 . 76 , 0 . 8 } , except for strong PTs where v w = 0 . 32 is excluded due to the non-existence of deflagrations for α ≳ 1 3 (1 -v w ) -13 / 10 [49], implying a total of 3 × 13 -1 = 38 PT parameter points. To extract our main results, we run reference simulations for each simulation box size, in which a single reference bubble nucleation history is used for all wall velocities, PT strengths, and grid sizes, thus keeping the sample variance in our reference measurements identical for different values of v w , α , and N .  \nThe bubble nucleation histories result in an (asymptotic) number of bubbles of the order of N b ≃ ˜ L 3 / (8 π v 3 w ), where ˜ L ≡ Lβ is the simulation box size, nucleated following a statistical distribution that is exponential in time and uniform in space, as described in Sec. 2.1, and then removing bubbles that nucleate inside the future causal cone of previous bubbles to take into account the evolution of the broken-phase volume with time (see Ref. [58] for details). In our simulations, ˜ L/v w takes on values of 20 and 40, yielding of the order of 300 and 2500 bubbles respectively. Using the same numerical resolution, simulations with ˜ L/v w = 40 yield a reduction in the statistical variance by increasing the number of bubbles and by offering an increased resolution of the measured quantities in the IR regime, while simulations with ˜ L/v w = 20 cover a larger dynamical range in the UV regime. For comparison, the number of bubbles for ˜ L/v w = 40 ( N b ≃ 2500) in our work and previous Higgsless simulations  \nTable 1 . Summary of simulation runs with physical and numerical parameter choices. Reference indicates the simulations constructed from a single reference bubble nucleation history (for each box size), thereby eliminating statistical differences among the sample of reference simulations. Seeds refers to simulations constructed from a set of 9 additional bubble nucleation histories, allowing to infer statistical sample variance for 3 selected wall velocities v w = 0 . 32 (0.36), 0.6, and 0.8, which correspond to a deflagration, a hybrid, and a detonation for weak and intermediate PTs, while v w = 0 . 8 is still a hybrid for strong PTs. Single-bubble refers to simulations with only one single centrally nucleated bubble, allowing us to study the convergence of self-similar profiles. We take a range of v w ∈ [0 . 32 , 0 . 8] in increments of 0.04 besides for strong transitions ( α = 0 . 5) for which we take v w ∈ [0 . 36 , 0 . 8]. A total of 1028 simulations have been performed.  \n×  \n[58] is in general larger than most 6 of the previous numerical simulations of the fluid-scalar system [50, 52, 77], especially for intermediate PTs, allowing for a reduction of the statistical variance. A potential issue of small box sizes is that for small wall velocities, the shock in front of the wall of the first nucleated bubble might collide with its mirror images (due to the use of periodic boundary conditions) before the end of the PT. To avoid this issue, we take the minimum value of v w to be 0.32 in our simulations, such that the numerical domain is filled with the broken phase before the largest bubble reaches the edges of the simulation box even for the smaller box ˜ L/v w = 20.  \nFor each of the 76 parameter points { v w , α, ˜ L/v w } , we then run simulations with different number of grid points N 3 with N ∈ { 64 , 128 , 256 , 512 } , yielding a total of 76 × 4 = 304 reference simulations. Running simulations of different grid sizes allows us to test the degree of convergence of our numerical results and to estimate physical quantities in the continuum limit by extrapolation (see Sec. 4.1). To ensure the stability of our simulations, we choose the number of time steps N t = ˜ t end /δ ˜ t to satisfy the Courant-Friedrichs-Lewy (CFL) condition δ ˜ t/δ ˜ x < 1 / 4 with δ ˜ x = ˜ L/N . We have confirmed that even for strong transitions, increasing N t beyond this threshold does not change the numerical results. For each parameter point { v w , α, ˜ L/v w } , we have also run single-bubble simulations to track the convergence of the selfsimilar fluid profiles, leading to 76 'single-bubble' simulations. These results are presented in App. B and will be used to improve the extrapolated predictions of the 'reference' multiplebubble simulations in Sec. 4.1. We note that single-bubble simulations are only run until ˜ t end = ˜ Lv w / [2 max( c s , v w )], being roughly the time when the fluid shell reaches the edge of the simulation domain.  \nIn addition to the reference simulations, we also run multiple-bubble simulations based on 9 additional distinct bubble nucleation histories per box size for all strengths, resolutions, and box sizes, for v w ∈ { 0 . 32 / 0 . 36 , 0 . 6 , 0 . 8 } , where the lower v w = 0 . 32 is used for weak and intermediate transitions, and v w = 0 . 36 for strong ones. These velocities correspond to deflagrations, hybrids, and detonations, respectively, except for strong transitions for which also v w = 0 . 8 corresponds to a hybrid. This implies a total of 3 × 3 × 2 × 4 × 9 = 648 seed simulations from which the statistical variance of the results can be estimated. We will use these simulations to provide error bars in our measured quantities, corresponding to the standard deviation from the 10 different bubble nucleation histories in Sec. 4.  \nAll reference and seed simulations are run between 0 < ˜ t ≡ tβ < 32 and the GW spectrum is extracted from the time interval spanning from ˜ t init = 16 to ˜ t end = 32. We set the origin of time coordinates ˜ t = 0 at a reference value such that the first bubble nucleates at ˜ t = 0 . 5, based on the invariance of our equations on time translations when the expansion of the Universe can be ignored. For this approximation to be valid we then require β/H ∗ ≫ ˜ t end = 32. We specifically cut out the early times up to ˜ t init = 16 to extract the contributions from the fluid perturbations after the collisions of bubbles, and to reduce the realization-dependent effects on the GW production. Consequently, we also suppress contributions to the GW spectrum from the initial collisions (see also the discussion in Ref. [56]). In this regime, we then compute I sim ( ˜ t init , ˜ t fin , ˜ k ) that allows us to robustly test the scaling of Eq. (2.19) and compute the GW efficiency ˜ Ω GW and the spectral shape S ( kR ∗ ). The time ˜ t init = 16 is shortly after the time when the broken phase fills up the whole volume of  \nFigure 1 . Velocity (upper panel), enthalpy fluctuations (middle panel), and vorticity (lower panel) in an xy -plane slice of the simulation volume at z = 0 and at different times ˜ t , for a strong PT with α = 0 . 5 and wall velocity v w = 0 . 36, which corresponds to a deflagration (see Fig. 13).  \n<!-- image -->  \nthe simulation, ˜ t 0 ≃ 10, for the reference nucleation history with ˜ L/v w = 20. We will consider times ˜ t > ˜ t 0 to fit the time evolution of the kinetic energy fraction K ( ˜ t ) = K 0 ( ˜ t/ ˜ t 0 ) -b in Sec. 4.2.  \nIn total, we have performed 1028 simulations, which we summarize in Tab. 1, with an estimated time of ∼ 10 6 CPUhours. We note that each large-resolution simulation ( N = 512) takes ∼ 10 3 CPU hours, a quite modest value that indicates the numerical efficiency of the Higgsless approach.", '4 Numerical results': "Before we present a detailed account of our numerical results, we would like to put them in perspective. Overall, our results can be summarized by the following findings:  \n- · Simulations of strong first-order PTs with α = 0 . 5: We present numerical results for strong PTs covering a wide range of wall velocities and performing systematic checks of the numerical convergence of our results. For the first time, we obtain the full GW spectra for strong PTs. 7 Stronger simulations are more challenging when it comes to  \nnumerical stability and proper resolution of non-linearities. At the same time, stronger PTs lead to a larger GW signal and therefore are preferred by a potential detection with LISA. Hence the importance of developing an accurate understanding of the resulting GW spectrum. We provide in Sec. 5 a template based on the expected GW spectrum from compressional fluid perturbations extended to decaying sources (see Secs. 2.4 and 2.5) to incorporate information from our simulations that can be used for phenomenological studies. We show in Fig. 1 an example of a simulation for a strong PT that corresponds to a deflagration with v w = 0 . 36.  \n- · Template parameterizations : Ideally, all our numerical findings can be expressed in terms of a few physical quantities to facilitate their use in phenomenological studies. In our simulations, all quantities evaluated are dimensionless such that β/H ∗ does not affect the numerical results and only appears when we recover the physical quantities, as indicated in Eq. (2.8). This motivated the authors in Refs. [56, 58] to use the variable Q ' [see Eq. (2.16)] to interpret the numerical results. In the present work, we instead characterize the numerical results based on R ∗ and K 2 int [see Eq. (2.19)], as to allow to capture the essential results in a form as simple as possible (with an almost invariant GW efficiency ˜ Ω GW ), while allowing for deviations with respect to the linear growth of the GW amplitude with the source duration, which is expected for stationary sources (see discussion in Secs. 2.4 and 2.5). We also provide in Sec. 2.6 a definition of K 2 int , exp that allows to incorporate a posteriori the effect of the expansion of the Universe [see Eq. (2.23)].\n- · Development of non-linearities : For strong PTs, and some intermediate PTs with confined hybrids, we observe several phenomena that probably stem from non-linear dynamics of the fluid. In the first place, we observe a decay in the kinetic energy of the fluid at later times (after the PT ends) that could indicate that non-linearities might be leading to a cascading from larger to smaller scales in the fluid perturbations, making the viscous dissipation at small scales more effective. Potentially due to this decay, we find that the dependence of the GW amplitude with the source duration starts to deviate from the expected linear growth, transitioning toward its saturation amplitude at the latest times of the simulations. The cascading of kinetic energy from large to small scales could also impact the UV part of the GW spectrum, leading to a modification from the expected k -3 found for sound waves [40, 50-52, 56, 58-60] towards a shallower spectrum, for example like the one that is found in vortical turbulence, k -8 / 3 [53, 55, 62, 63, 65, 71, 95]. We test the numerical robustness of our results along this section and comment on future studies that would be required to confirm some of our findings. In particular, we present in App. C a preliminary study of the potential development of vorticity in our simulations, and show the vorticity found for an example simulation in the lower panel of Fig. 1.  \nTo test the validity of our numerical results, we will pay special attention to the following points, to be addressed throughout this section. In Sec. 4.1, we study the convergence of our results with respect to the grid spacing, δ ˜ x . In Sec. 4.2, we study the time dependence of the fluid kinetic energy fraction K , and fit the numerical results to the decaying power-law presented in Sec. 2.5. In Sec. 4.3, we test the expected scaling of the GW spectrum with K 2 int and R ∗ , evaluating the evolution of the integrated GW amplitude with the source duration, and compute the GW efficiency ˜ Ω GW , according to Eq. (2.19). We also provide in Sec. 4.3  \nan estimate of the expected GW amplitude in a flat Minkowski space-time, based on the numerical results of Secs. 4.1-4.3, and in an expanding background, using the model presented in Sec. 2.6. Finally, in Sec. 4.4, we study the spectral shape of the GW spectrum, paying special attention to the UV regime, where we find deviations with respect to the expected slope in the sound-wave regime. However, we note that to be able to confirm the presence of a forward cascade in our simulations, and therefore the accuracy of the fits provided for the time decay of K and the resulting GW amplitude, a detailed study of the dependence of the kinetic spectra on the numerical parameters would be required, checking whether the cascading within an inertial range of scales is developed and unaffected by the numerical parameters. Furthermore, one should keep in mind that it is not in general expected that all wave numbers evolve with the source duration in the same way, as shown in Ref. [59].", '4.1 Convergence analysis of the kinetic energy and GW amplitude': "In the present study, for each parameter point { α, v w } , simulations at four resolutions N ∈ { 64 , 128 , 256 , 512 } and two box sizes ˜ L/v w ∈ { 20 , 40 } have been performed. Since our Higgsless simulations use relatively sparse grids compared to simulations with scalar fields [40, 50, 52, 77], the resolution can become an issue to reproduce some of the expected self-similar profiles induced by the uncollided nucleated bubbles at the initial stages of the simulations, especially for parameter points with v w ≲ v CJ , where v CJ is the Chapman-Jouguet speed, determining the transition between hybrids and detonations. In these situations, the fluid profiles become very thin hybrid profiles as v w approaches v CJ . Since the Chapman-Jouguet speed is v CJ = { 0 . 63 , 0 . 73 , 0 . 89 } for α = { 0 . 0046 , 0 . 05 , 0 . 5 } , respectively, for our choice of parameters, we have very thin profiles when v w = 0 . 6 and the PT is weak, when v w = 0 . 72 and the PT has intermediate strength, and when v w = 0 . 8 and the PT is strong. Especially in these cases, the resolution in ξ ≡ r/ ( t -t n ), with r being the radial distance to the nucleation location and t n the time of nucleation, might not be enough to resolve the self-similar profiles at the time when the bubble collides, since the resolution in ξ for a fixed N is initially low and then improves as time evolves. We study this rate of convergence in App. B using singlebubble simulations to understand for each PT what is the required time for the self-similar profiles to converge to the expected ones (see also discussion in Ref. [58]). For reference, we show in App. B (see Fig. 13) the self-similar profiles of the fluid perturbations expected to be produced by uncollided expanding bubbles [49], computed using CosmoGW [96].  \nIn the following, we first analyze the convergence of the numerical results for multiplebubble reference runs to provide estimates of the integrated kinetic energy K int ( ˜ t init , ˜ t end ) and GW amplitude I int sim ≡ ∫ I sim dln k , where I sim ( ˜ t init , ˜ t end , ˜ k ) is defined in Eq. (2.10), for initial and final times ˜ t init = 16 and ˜ t end = 32. These integrated quantities will be used in Secs. 4.3 and 4.4 to estimate respectively the GW efficiency ˜ Ω GW and the spectral shape S ( kR ∗ ). We also provide estimates for the kinetic energy fraction K 0 evaluated at the time when the PT completes, ˜ t 0 ≃ 10. We will then attempt to improve this estimate by including the results of single-bubble runs studied in App. B, tracking the degree of convergence of each bubble at the time when they collide, and leading to a new estimate, K 0 , defined in Eq. (B.4) (see App. B for details). We note that we will only use the improved estimates K 0 to make predictions of the resulting GW template presented in Sec. 5, while we study K 0 computed from the reference runs for the remaining of this section.  \nIn addition to the required resolution for thin self-similar profiles before collisions, as the fluid perturbations become non-linear, we expect large numerical resolutions to be required to fully capture the dynamics during and after collisions. In order to take these effects  \ninto account, we study the numerical results as a function of N and attempt to potentially improve our measurements by extrapolating our results to N →∞ , based on the underlying assumption that the extrapolation method obtained for the computed values of N also applies in this limit. On general grounds, it is possible to define a particular number of grid points N ∗ such that for N ≫ N ∗ a simulation has reached a converged solution, i.e., the numerical results are unaffected within some acceptable tolerance and, hence, we can assume that they accurately represent their continuum values. Empirically, we find that for insufficient resolution, the kinetic energy fraction K is in general underestimated when the grid resolution is insufficient as the values of the velocity profiles around the peak are underresolved. 8 This motivates us to use the following function when extrapolating the numerical values of the kinetic energy fraction:  \nK = K ∞ 1 + ( N ∗ /N ) a , (4.1)  \nwhere a , N ∗ , and K ∞ are found by fitting the numerical results as a function of N . In the few cases when this fit is not valid, in particular for the largest values of K (see footnote 8), we take K ∞ to be the value computed in the simulations for N = 512, and the error ε K is then estimated comparing this value to the one obtained for N = 256. We note that the value of a in Eq. (4.1) indicates the degree of convergence of the numerical results, as the relative error in the kinetic energy ε K can be expanded as  \nε K ≡ K ∞ -K K ∞ = ( δ ˜ x/δ ˜ x ∗ ) a + O ( δ ˜ x 2 a ) , (4.2)  \nwhere δ ˜ x ∗ = ˜ L/N ∗ . In some simulations, as we discuss later, we can observe that the numerical results have already converged for N = 512, indicated by a small relative error ε K ≡ | K -K ∞ | /K ∞ (see values in Tab. 2).  \nWe find empirically that the exponent a in (4.1) usually varies between one and two in our numerical simulations, indicating that the dynamics of the system reduces the effective degree of convergence with respect to the one expected from the numerical scheme, which  \nThe result of applying the convergence analysis based on Eq. (4.1) to the reference runs with multiple bubbles is shown in Fig. 2, where we first consider the value of the kinetic energy fraction K 0 (left panels). For the full evolution of the kinetic energy, see Sec. 4.2 and Fig. 4. The middle panels of Fig. 2 correspond to the rms kinetic energy fraction, computed from K int [see Eq. (2.18)] as ˜ T GW K 2 rms ≡ K 2 int , where ˜ T GW ≡ ˜ t end -˜ t init = 16 corresponds to the time interval over which the GW spectrum is computed. We note that we focus on the simulations with ˜ L/v w = 20 as these cover smaller scales than those with ˜ L/v w = 40 for a fixed N , providing a better resolution of the kinetic energy density in the UV regime. We show in Tab. 2 the resulting values of the fit parameters of Eq. (4.1) for the rms kinetic energy fraction K rms ∞ , normalized by the single-bubble kinetic energy fractions K ξ [see Eq. (2.17)], and a K , and the relative errors ε K as defined in Eq. (4.2) (see also footnote 8), for the set of PTs shown in Fig. 2. We expect K 2 rms ˜ T GW to be the relevant kinetic energy fraction entering the computed GW amplitude, I int sim , according to the model proposed in Sec. 2.5.  \nξ  \n/K  \n0  \nK  \nξ  \n/K  \nrms  \nK  \n2  \n.  \n00  \n1  \n.  \n75  \n1  \n.  \n50  \n)  \nβ  \n∗  \nR  \n2  \nξ  \nK  \nGW  \n˜  \nT  \n(  \n/  \nt  \nin  \nsim  \nI  \n1  \n.  \n25  \n1  \n.  \n00  \n0  \n.  \n75  \n×  \n2  \n10  \n0  \n.  \n50  \n0  \n.  \n25  \n0  \n.  \n00  \n1  \n.  \n2  \n1  \n.  \n0  \n0  \n.  \n8  \n0  \n.  \n6  \n0  \n.  \n4  \n0  \n.  \n2  \n0  \n.  \n0  \n1  \n.  \n2  \n1  \n.  \n0  \n0  \n.  \n8  \n0  \n.  \n6  \n0  \n.  \n4  \n0  \n.  \n2  \n0  \n.  \n0  \n0  \n.  \n0  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \n0  \n.  \n0  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \n0  \n.  \n0  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \nx/vw  \nα  \n= 0  \n.  \n0046  \nδ  \n˜  \nδ  \n˜  \nδ  \n˜  \nδ  \n˜  \nδ  \n˜  \n0  \n.  \n0  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \n0  \n.  \n0  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \nx/vw  \nα  \n= 0  \n.  \n5  \nδ  \n˜  \nδ  \n˜  \nδ  \n˜  \n0  \n.  \n0  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \n0  \n.  \n0  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \n0  \n.  \n0  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \nx/vw  \nx/vw  \nx/vw  \nx/vw  \nx/vw  \nx/vw  \n0.36  \n0.44  \n0.52  \n0.6  \n0.68  \n0.76  \nFigure 2 . Plots showing the kinetic energy fraction K and the integrated GW spectrum I int sim as a function of grid spacing δ ˜ x/v w = ( ˜ L/v w ) /N for runs with ˜ L/v w = 20. Upper and middle panels respectively show K at the time when all the simulation domain is in the broken phase around ˜ t 0 ≃ 10, K 0 , and its rms value, K rms ≡ K int / ˜ T 1 / 2 GW , both normalized by the single-bubble kinetic energy fraction K ξ [see Eq. (2.17)]. Lower panel shows I int sim normalized by a reference value ˜ Ω GW ∼ 10 -2 [40, 50, 52], and by ˜ T GW K 2 ξ R ∗ β , based on the expected scaling of Eq. (2.19). Both K rms and I int sim are computed for ˜ t init = 16 and ˜ t end = 32, with ˜ T GW = 16. Left, middle, and right columns are weak, intermediate, and strong PTs respectively. Solid lines show the least-squares fits of the extrapolation scheme given in Eq. (4.1) when the fit is valid, while dotted lines indicate the numerical trend when the fit is not valid (see footnote 8). Stars indicate the extrapolated values at δ ˜ x → 0.  \n<!-- image -->  \n1  \n.  \n2  \n1  \n.  \n0  \n0  \n.  \n8  \n0  \n.  \n6  \n0  \n.  \n4  \n0  \n.  \n2  \n0  \n.  \n0  \n1  \n.  \n4  \n1  \n.  \n2  \nξ  \n/K  \nrms  \nK  \n)  \nβ  \n∗  \nR  \n2  \nξ  \nK  \nGW  \n˜  \nT  \n(  \n/  \nt  \nin  \nsim  \nI  \n1  \n.  \n0  \n0  \n.  \n8  \n0  \n.  \n6  \n×  \n2  \n10  \n0  \n.  \n4  \n0  \n.  \n2  \n0  \n.  \n0  \nξ  \n/K  \n0  \nK  \nξ  \n/K  \nrms  \nK  \n0  \n.  \n06  \n0  \n.  \n05  \n)  \nβ  \n∗  \nR  \n2  \nξ  \nK  \nGW  \n˜  \nT  \n(  \n/  \nt  \nin  \nsim  \nI  \n0  \n.  \n02  \n×  \n2  \n10  \n0  \n.  \n04  \n0  \n.  \n03  \n0  \n.  \n01  \n0  \n.  \n00  \n1  \n.  \n2  \n1  \n.  \n0  \n0  \n.  \n8  \n0  \n.  \n6  \n0  \n.  \n4  \n0  \n.  \n2  \n0  \n.  \n0  \n1  \n.  \n2  \n1  \n.  \n0  \n0  \n.  \n8  \n0  \n.  \n6  \n0  \n.  \n4  \n0  \n.  \n2  \n0  \n.  \n0  \n0.36  \n0.44  \n0.52  \n0.6  \n0.68  \n0.76  \ncorresponds to second order [56, 58]. We expect that further decreasing δ ˜ x would be required to finding an exact quadratic dependence of the error with δ ˜ x . In any case, we note that for most of the PTs (besides highly confined profiles with v w ≲ v CJ ), we already find absolute errors below 10%, as indicated in Tab. 2. For confined profiles, the relative error is large, and we need to take into account that the extrapolated result K ∞ presents a larger degree of uncertainty. In these simulations, we expect the lack of convergence to also become visible in the GW spectra: For example, we observe that the expected UV behavior, S ( k ) ∼ k -3 , found in the sound-shell model [51, 54, 59, 60] and in PT simulations [40, 50, 52] is in these cases obscured by an exponential decay (see the discussion in Sec. 4.4 and the fit used in Ref. [58]).  \nFurthermore, we also apply the extrapolation of Eq. (4.1) to the integrated GW spectrum obtained in the code I int sim . The resulting values as a function of δ ˜ x are shown in the lower panels of Fig. 2, while the numerical values of the extrapolated I int ∞ , the fit parameter a I , and the relative error ε I are given in Tab. 2. We note that for the GW amplitude, the fit of Eq. (4.1) is not valid when the PT is strong for most of the wall velocities, while it remains valid for weak and intermediate PTs. In the cases when the fit is not valid, we take the extrapolated values in Fig. 2 as those obtained for the largest resolution runs with N = 512. However, the relative errors comparing simulations with N = 512 and N = 256 are already very small compared to the errors in K , making the estimate of I int sim less sensitive to numerical inaccuracies than the estimate of K int .  \nFinally, we also display in Fig. 2 the kinetic energy fraction K 0 at the reference time ˜ t 0 ≃ 10. As we will see later, K 0 is essential to determine the resulting GW amplitude, as can be observed from Eqs. (2.19) and (2.20). Therefore, to correctly capture the GW amplitude we need to accurately reproduce K 0 . A first attempt is to directly take the extrapolated values K 0 ∞ using the fit of Eq. (4.1) (see upper panel of Fig. 2). However, this result does not  \nTable 2 . Numerical values of the fit parameters a , K rms ∞ , and I int ∞ of Eq. (4.1) for the rms value of the kinetic energy fraction K rms /K ξ and the integrated GW amplitude I int sim , respectively shown in the middle and lower panels of Fig. 2. We present the relative errors ε computed by comparing the extrapolated values to those obtained in the largest resolution runs N = 512 when the fit works (see footnote 8). Otherwise (indicated with '-' in the values of a ), the error is computed from the two largest resolution runs, N = 256 and N = 512.  \ntake into account that the self-similar profiles have not reached convergence at the time when the bubbles collide, leading to underestimating the ratio K 0 /K ξ . We present a methodology that estimates the required correction for this underresolution in App. B. The corrected values K 0 /K ξ are presented in Fig. 3, together with a modified 'efficiency' κ 0 , such that  \nK 0 ≡ κ 0 α 1 + α , (4.3)  \nin analogy to Eq. (2.17). We find a general trend that K 0 /K ξ ≳ 1 when v w < c s , while K 0 /K ξ ≲ 1 when v w > c s . If we take an average value of K 0 over the PT parameters v w and α , we find  \nK 0 = 0 . 84 +0 . 24 -0 . 29 K ξ , (4.4)  \nwhere the super and subscripts are the maximum and minimum values found over all wall velocities, indicating that the typical use of K ξ for the kinetic energy would overestimate the GW production by a factor ( K ξ / K 0 ) 2 , which can be as large as 0 . 55 -2 ∼ 3 . 3, for example when α = 0 . 5 and v w = 0 . 8 (see Fig. 2). For different PT parameters ( α and v w ), one could take the values of Fig. 3 to predict the correction of the resulting GW amplitude. According to the locally stationary UETC presented in Sec. 2.5, we find that K 2 int ≡ K 2 rms ˜ T GW is the relevant quantity determining the GW amplitude. Therefore, based on the results of Ref. [58], Ref. [92] used the estimated value K rms ≃ 0 . 6 K ξ . However, we note that using the power-law fit presented in Sec. 2.5 and validated in Sec. 4.2, we can relate the final GW amplitude to the corrected values K 0 and the decay rate b (as we will do in Sec. 4.3).", '4.2 Time evolution of the kinetic energy': 'In this section, we evaluate the time evolution of the kinetic energy fraction K for different numerical resolutions N . We show the results in Fig. 4 for the largest resolution runs N = 512 (upper panel), and for a range of N = { 64 , 128 , 256 , 512 } (lower panel). We find that the kinetic energy is underestimated for low resolution, as argued in Sec. 4.1 (see upper and middle panels of Fig. 2, and lower panel of Fig. 4). This is, for most PT parameters, the case at early times, before the bubbles collide and while the self-similar profiles develop in each nucleated bubble (see App. B). At later times (usually ˜ t ≳ ˜ t 0 ≈ 10), for weak and most of intermediate PTs we find a decay of the kinetic energy with time that becomes less pronounced as we increase the resolution N , directly related to the underresolution of K at earlier times. Since the kinetic energy is typically damped by numerical viscosity, 9 it is in general expected that the decay is less pronounced when the grid spacing is reduced, as can be observed for the weak and most of intermediate PTs.  \nBesides this general trend, the opposite is found for strong PTs, and intermediate ones with thin hybrid profiles ( v w ≲ v CJ ), such that in the decaying phase of the kinetic energy, the decay becomes steeper with smaller grid spacing. To study the dependence of the decay rate with resolution, we show in the lower panel of Fig. 4 the evolution of K ( ˜ t ) normalized by the corresponding values of K 0 , found using the fit presented below. In these cases, the enhancement of the decay with resolution might indicate that as the fluid shells carry larger kinetic energies at the time of collisions, non-linearities are enhanced and might eventually overcome the effect of numerical dissipation. From energy conservation, we then expect that as non-linearities develop, kinetic energy transfers from larger to smaller scales where it can be converted to thermal energy at the scale determined by numerical viscosity.  \n<!-- image -->  \nFigure 3 . Upper panel: Kinetic energy efficiency κ 0 ≡ K 0 (1 + α ) /α obtained from correcting the numerical values obtained for N = 512 (see App. B) for weak (blue), intermediate (red), and strong (orange) PTs, compared to κ ξ (black) for self-similar solutions [see Eq. (2.17)]. Vertical line corresponds to c s , and v CJ is indicated by the dotted gray line. Lower panel: We plot the values of the single-bubble-corrected kinetic energy fraction at the time when the PT completes, K 0 , normalized to the values K ξ for self-similar profiles [see Eq. (2.17)], for weak (left panel), intermediate (middle panel), and strong (right panel) PTs, as a function of v w . Lines in increasing opacity correspond to increasing numerical resolution N ∈ { 64 , 128 , 256 , 512 } . The vertical solid gray line indicates the sound speed, c s , while the dashed lines indicate the Chapman-Jouguet velocity, v CJ . Error bars show the standard deviation from 10 different bubble nucleation histories.  \n<!-- image -->  \nIn addition to the time decay of K , we also find oscillations in time that can be associated to the sound-wave regime, where an oscillatory conversion between kinetic and thermal energies is expected, and confirmed by the fact that we conserve T 00 to machine precision (see results in App. C of Ref. [58]). We then fit the numerical results at times ˜ t > ˜ t 0 when the PT is complete, using the following power-law decay with time, effectively getting rid of the oscillations over time,  \nK ( ˜ t > ˜ t 0 ) = K 0 ( ˜ t ˜ t 0 ) -b , (4.5)  \nFigure 4 . Evolution of the measured kinetic energy fraction K ( ˜ t ) normalized to the single-bubble values K ξ [see Eq. (2.17)] for weak (left columns), intermediate (middle columns), and strong (right columns) PTs, for N = 512 (solid lines) and ˜ L/v w = 20, and the same wall velocities as those in Fig. 2. Dashed lines indicate the fits to the power-law decay of Eq. (4.5) at times ˜ t > ˜ t 0 . Values corresponding to K 0 /K ξ are marked with circles. In the lower panel, the kinetic energy fraction is shown for different numerical discretizations N = { 64 , 128 , 256 , 512 } (solid lines with increasing opacity), normalized to the corresponding values of the fit K 0 at each resolution N . The results for each v w are shifted by a constant to distinguish between wall velocities. This presentation in the lower panel is chosen to emphasize the dependence of the time decay on resolution.  \n<!-- image -->  \nwhere b indicates the power-law decay rate of K . This power-law decay prescription accurately fits the numerical data (see Fig. 4), and we have checked that it remains accurate up to ˜ t end = 64 for an example strong PT with α = 0 . 5 and v w = 0 . 8. We define the half-life of the kinetic energy as the time when K ( ˜ t 0 + ˜ t 1 / 2 ) = 1 2 K 0 , i.e.,  \n˜ t 1 / 2 = ( 2 1 b -1 ) ˜ t 0 . (4.6)  \nb  \n0  \n.  \n0  \n0  \n.  \n2  \n0  \n.  \n4  \n0  \n.  \n6  \n0  \n.  \n8  \n1  \n.  \n0  \n1  \n.  \n2  \n1  \n.  \n4  \n0  \n.  \n3  \n0  \n.  \n4  \n0  \n.  \n5  \n0  \n.  \n6  \n0  \n.  \n7  \n0  \n.  \n8  \nα  \n= 0.0046  \nα  \n= 0.05  \nα  \n= 0.5  \nDecay index  \nb  \nvw  \nFigure 5 . Decay index b (left panel) and half-life ˜ t 1 / 2 (right panel) as a function of v w for N = { 256 , 512 } in increasing opacity for weak (blue lines), intermediate (red lines), and strong (orange lines) PTs. Dashed black lines with colored stars in the right panel correspond to the eddy turnover time ˜ t eddy = ( βR ∗ ) / √ K ξ that we compare with ˜ t 1 / 2 as we expect both time scales to be inversely proportional to K ξ . Error bars in the left panel show the standard deviation from 10 different bubble nucleation histories for N = 512.  \n<!-- image -->  \nWe display in Fig. 5 the fit of the decay index b (left) as well as the half-life ˜ t 1 / 2 (right) as a function of v w for weak, intermediate, and strong PTs. 10 In the right panel of Fig. 5, we also plot the eddy turnover time based on the kinetic energy ratio expected for uncollided bubbles, ˜ t eddy = ( βR ∗ ) / √ K ξ , which corresponds to the time scale of fluctuations in the plasma and it is expected to determine the time decay into turbulent motion. We compare ˜ t 1 / 2 to ˜ t eddy in Fig. 5. The eddy turnover time is ˜ t eddy ≃ 5 for strong PTs, ˜ t eddy ≃ 10-30 for intermediate PTs, and ˜ t eddy ∼ O (100) for weak PTs. Therefore, we expect that the time scale for non-linearities to develop is reached within the simulations for strong and some intermediate PTs, while for other intermediate PTs, the eddy turnover time occurs towards the end of our simulations. To evaluate the development of vortical motion in our simulations, we briefly discuss in App. C the presence of vorticity in our simulations and present some preliminary results.  \nFor weak transitions, the rate of kinetic energy damping is greatly reduced as we increase the resolution, which we interpret as a reduction of the numerical viscosity (see footnote 9).  \nThis observation therefore means that for weak transitions, decay is always dominated by numerical viscosity. Only for the hybrid solution with v w = 0 . 6 ≲ v CJ , when larger fluid velocities can be achieved (see self-similar profiles in Fig. 13), does b (and hence ˜ t 1 / 2 ) appear to stagnate with increasing resolution, pointing towards the onset of resolving the physics responsible for the damping. However, in this extreme case the fluid profile is highly confined and the simulations are far from reaching the converged profiles (see App. B), so it is not completely clear whether the obtained decay rate b is physical.  \nThe results are more interesting in the case of intermediate transitions. For both small and large v w , corresponding to subsonic deflagrations and detonations respectively, b decreases with increasing resolution. However, for a large range of intermediate velocities v w ∈ { 0 . 52 , 0 . 6 , 0 . 68 } ≲ v CJ , the trend is reversed for the highest resolutions. We interpret this point of reversal as a transition from a decay of the kinetic energy dominated by numerical viscosity to a decay determined by the development of non-linearities. We note that in this case, some of the confined hybrids are still underresolved but this is no longer the case for the subsonic deflagration with v w = 0 . 52, indicating that the decay rate seems to be physical (see Fig. 4 and App. B). Furthermore, a similar decay of the kinetic energy was already found for intermediate PT simulations of the scalar-fluid system [52, 77].  \nFor strong transitions, we are universally in the regime where increasing the numerical resolution N leads to a larger decay rate, indicating that the physical non-linear decay dominates over the numerical viscosity. As discussed above, this is expected to be the case, as the expected time scale for non-linearities to develop, i.e., the eddy turnover time, is around ˜ t eddy ≃ 5, occurring during the duration of our numerical simulations.', '4.3 Time evolution of the integrated GW spectrum and GW efficiency': "We show in Fig. 6 (upper panel) the time evolution of the integrated GW amplitude I int sim ≡ ∫ I sim dln k , where I sim ( ˜ t init , ˜ t, k ) is evaluated at ˜ t init = 16, and we allow ˜ t to vary from ˜ t init to ˜ t end = 32 [see Eq. (2.10)]. We find that for weak and intermediate PTs, the evolution with the source duration ˜ t -˜ t init is close to linear in most cases (unless v w ≲ v CJ ), as expected from the usual stationary assumption in the sound-wave regime (see Sec. 2.4) and as argued in previous numerical work [40, 50, 52, 56, 58]. In these situations, K does not significantly evolve with time within the simulations (see Fig. 4). However, when the decay of K is significant, we observe deviations with respect to the linear growth, as expected from the generalized locally stationary UETC proposed in Sec. 2.5. To test the validity of Eq. (2.19), found under this assumption, we plot in the lower panels of Fig. 6 the following ratio  \n˜ Ω GW ( ˜ t ) = I int sim ( ˜ t init , ˜ t ) K 2 int ( ˜ t init , ˜ t )( βR ∗ ) , (4.7)  \nwith the objective to estimate the GW efficiency ˜ Ω GW while including the effect of the decay of K . We note that when K does not significantly decay with time, we recover the linear growth, K 2 int → K 2 ˜ T GW .  \nFor α = 0 . 0046 and v w = 0 . 6, which corresponds to a confined hybrid (see Fig. 13), K presents a sharper decay with time than the other PTs (see Fig. 4) that translates into a decay of the GW amplitude with respect to the linear growth. A similar situation occurs for confined hybrids with α = 0 . 05 and intermediate wall velocities v w ∈ { 0 . 6 , 0 . 68 } ≲ v CJ . We see in Fig. 6 that in these cases again the GW amplitude grows slower than linearly with the sourcing time. However, the proposed ratio ˜ Ω GW ( ˜ t ) of Eq. (4.7), defined with K 2 int , is closer  \nFigure 6 . Dependence of the numerical integrated GW amplitude found in the simulations with ˜ L/v w = 40 and N = 512 as a function of the source duration ˜ t -˜ t init for weak (left columns), intermediate (middle columns), and strong (right columns) PTs. The integrated GW amplitude is normalized as in the lower panel of Fig. 2, while ˜ Ω GW is shown in the lower panel, computed using Eq. (4.7). Dashed lines in the upper panel correspond to the linear growth with the source duration expected under the stationary UETC assumption.  \n<!-- image -->  \n-  \n-  \n-  \nto be constant with time (see lower panel of Fig. 6), allowing us to still compute the GW efficiency and spectral shape. We note that the initial time growth of ˜ Ω GW is a numerical artifact from abruptly starting the GW computation at ˜ t init . For strong PTs, the kinetic energy decays significantly for all wall velocities within the time of our simulations. We then see that the growth of the GW amplitude deviates soon from the linear growth and one needs to incorporate the effect of the time evolution of K . Again, we show in the lower panels of Fig. 6 that the ratio ˜ Ω GW ( ˜ t ) is very close to constant, validating the generalization of the linear growth to a growth proportional to K 2 int ( ˜ t ) as found by the model proposed in Sec. 2.5. Therefore, as long as the fit K ( ˜ t ' ) = K 0 ( ˜ t ' / ˜ t 0 ) -b accurately represents the numerical results at times ˜ t ' ∈ [ ˜ t init , ˜ t ] (see Fig. 4), the growth of the GW amplitude is  \nI int sim ( ˜ t init , ˜ t ) = ˜ Ω GW K 2 0 ( βR ∗ ) ˜ t init ( ˜ t 0 / ˜ t init ) 2 b [1 + ( ˜ t -˜ t init ) / ˜ t init ] 1 -2 b -1 1 -2 b , (4.8)  \nas found in Eq. (2.20) for the case ˜ t 0 = ˜ t init .  \nSince we find that ˜ Ω GW ( ˜ t ) is roughly constant in time after incorporating K 2 int in the scaling of the GW amplitude (see Fig. 6), we take this value at the end of the simulations ˜ t end . The resulting GW efficiency ˜ Ω GW is shown in Fig. 7 for different numerical resolutions N and for both box sizes ˜ L/v w = 20 and 40. We show the values of ˜ Ω GW , computed from Eq. (2.19)  \nFigure 7 . Gravitational wave production efficiency ˜ Ω GW for weak (left), intermediate (middle), and strong (right) first-order PTs. Solid (dotted) lines correspond to ˜ L/v w = 20 (40). Black lines with increasing opacity correspond to increasing resolutions N ∈ { 64 , 128 , 256 , 512 } , while colored lines are the values extrapolated to the continuum limit δ ˜ x → 0. Dots and stars mark ˜ Ω GW as presented in Tabs. 2 and 3 of Ref. [54] corresponding to predictions from the sound-shell model (SSM) for exponential nucleation of bubbles [54] and scalar field-hydrodynamical simulations for simultaneous nucleation [52], respectively. Gray dots correspond to SSM values found using the assumption described in Sec. 2.4, following App. B of Ref. [59], and computed using CosmoGW [96]. Error bars indicate the standard deviation from 10 different bubble nucleation histories for ˜ L/v w = 20 (darker) and 40 (lighter).  \n<!-- image -->  \nusing the extrapolated values to δ ˜ x → 0 of I int sim and K int , as described in Sec. 4.1 (see also Fig. 2). We compare the extrapolated efficiencies with those found using the sound-shell model [51, 54], under the assumptions described in Sec. 2.4 (see also App. B of Ref. [59]), and those obtained from numerical simulations of the full coupled scalar field-fluid system [50, 52]. However, we note that the latter are found using simultaneous bubble nucleation, which in general leads to smaller values of ˜ Ω GW compared to exponential nucleation (see Tabs. 2 and 3 in Ref. [54]). We have modified the values of ˜ Ω GW from Refs. [50, 52, 54] to take into account that they consider βR ∗ = (8 π ) 1 / 3 v w , instead of the corrected βR ∗ = (8 π ) 1 / 3 max( v w , c s ) that we use in Eq. (2.19). Furthermore, we note that the integrated ˜ Ω GW might be modified by the structure that develops below the peak, described in Ref. [59]. However, we neglect this effect for two reasons: (1) we expect that the inclusion of the small wave numbers in ˜ Ω GW is negligible when we are under the assumptions described in Sec. 2.4, (2) the dynamical range in the IR available in our and previous simulations is usually not large enough to clearly reconstruct the exact spectral shape described in Ref. [59], see results in Sec. 4.4.  \nOverall, we find reasonably good agreement for weak PTs, while deviations start to be more significant for intermediate and strong PTs. The error bars in Fig. 7 show the standard deviation obtained from 10 different bubble nucleation histories, corresponding to the 'seeds' set of simulations listed in Tab. 1. For weak transitions ( α = 0 . 0046), the extrapolated values obtained from the Higgsless simulations accurately reproduce not just the numerical values but also the trend of ˜ Ω GW with v w observed in both the sound-shell model and the coupled scalar field-hydrodynamical simulations. This is important for two reasons: (1) the  \nagreement between three independent approaches lends support to the conclusion that the general trend may be physical; (2) since weak transitions are expected to be described by linear dynamics, limit in which the sound-shell model applies, we would a priori expect the Higgsless simulations to accurately reproduce the sound-shell model results, validating this model. However, as α becomes larger, non-linearities become more relevant and full 3D simulations are necessary to push beyond the reach of the sound-shell model. Only a few points of reference data for ˜ Ω GW exist for intermediate PTs ( α = 0 . 05) and so far none 11 for strong PTs ( α = 0 . 5). We note that reference data points ˜ Ω GW in Refs. [52, 54] are computed assuming a linear growth with the source duration as in Eq. (2.14). Hence, incorporating K int as in Eq. (2.19) can modify the value of ˜ Ω GW when the source decays. The extrapolation method described in Sec. 4.1 and presented in Fig. 7 as solid lines seem to behave very well, delivering agreement between the numerical results from both simulation domains ˜ L/v w = 20 and 40.  \nFor intermediate PTs, we begin to see deviations from the sound-shell model, in particular for v w = 0 . 68 ≲ v CJ . We observe that the v w -dependence seen for weak transitions has flattened and that the overall efficiency ˜ Ω GW is larger. Our findings are consistent with the two available data points for scalar field-hydrodynamical simulations from Ref. [52], indicating a departure from linearity in the fluid perturbations and, hence from the sound-shell model. We note that discrepancies with the numerical results of Ref. [52] might be due to the different nucleation histories considered (simultaneous in Ref. [52] and exponential in our simulations). Again, extrapolation seems overall good as the extrapolated values agree well for the simulations with ˜ L/v w = 20 and 40. For strong PTs, we observe even larger efficiencies overall. Besides, the impact of the non-linearities seems to wash out again the dependence of ˜ Ω GW on the wall velocity.  \nWe note that for weak PTs, the relative difference between the extrapolated values and those obtained in our largest resolution runs with N = 512 is still large. Hence, the exact values provided in Fig. 7 might still present numerical errors related to those listed in Tab. 2. Indeed, we find these potential errors to be larger for weak PTs (up to 50% for extremely thin profiles, and usually below 10% otherwise), where we can compare our extrapolated results to those found by the sound-shell model, while for intermediate and strong PTs, the errors of our extrapolated values are below 10% for all wall velocities.  \nAs a final note, we point out that the definition of ˜ Ω GW in terms of the integrated kinetic energy K 2 int reduced the dependence on the wall velocity and the strength of the PT significantly - compared to normalizing it to a stationary kinetic energy ratio (e.g., the one found for the self-similar bubbles or K 0 ) multiplied by the source duration. Partially, this is due to the decay found in the kinetic energy that is not caputred by the stationary assumption for the UETC (see discussion in Sec. 2.5 and numerical results in Sec. 4.2). Furthermore, the universality of ˜ Ω GW is also due to the use of Eq. (2.19) instead of Q ' /K 2 considered in previous work [56, 58], as discussed in footnote 4.  \nThe average values of ˜ Ω GW over v w for each strength α from the simulations are the  \nfollowing  \nwhere the super and subscripts refer to the maximum and minimum values found in the extrapolated values from our simulations. 12 We note that these values only take into account variations of ˜ Ω GW at different values of v w , but not from numerical inaccuracies in our numerical results and hence, in the extrapolated values.  \n10 2 ˜ Ω GW =      1 . 04 +0 . 81 -0 . 67 , for α = 0 . 0046 ; 1 . 64 +0 . 29 -0 . 13 , for α = 0 . 05 ; 3 . 11 +0 . 25 -0 . 19 , for α = 0 . 5 , (4.9)  \nFinally, we note that for most of weak and intermediate PTs we still find a growth rate with the source duration close to linear and, hence, we have not reached the free-propagation regime of the GW amplitude. Therefore, we need to make the usual assumption that the linear growth will persist until the development of non-linearities at ˜ t -˜ t init ≡ ˜ τ sw ∼ ˜ t eddy = ( βR ∗ ) / √ K to then saturate at that time. For PTs where the non-linearities timescale has been reached within the duration of the simulation (strong and some intermediate PTs), we find that even though the GW amplitude is growing with the source duration slower than linearly, it is still growing after ˜ t eddy , and until the final time of our simulations. Based on the decay found for K 2 ( ˜ t ) ∼ ˜ t -2 b in Sec. 4.2, the GW amplitude of Eq. (4.8) asymptotically grows proportional to I int sim ∼ ˜ t 1 -2 b for b < 1 2 , I int sim ∼ ln ˜ t for b = 1 2 , and I int sim ∼ ˜ t 0 for b > 1 2 . Then, we need to extrapolate the resulting GW amplitude by extending K 2 int to times beyond the final time of the simulation using Eq. (2.19). We note that unless b > 1 2 , then the GWamplitude keeps growing unbounded as long as the UETC assumed in Sec. 2.4 describes the source dynamics. However, we expect that the UETC deviates from this description as vortical motion and turbulence development dominates in the simulation [53, 64, 71, 90]. This can effectively be modelled by an appropriate choice of the source duration ˜ τ sw at which to stop the GW sourcing, which we leave as a free parameter in our current estimates. We emphasize that these results seem to indicate that, after the fluid perturbations enter the non-linear regime, the GW amplitude still takes some time to saturate to its free-propagation value and, hence, assuming a linear growth that is cut at ˜ τ sw = ˜ t eddy would underestimate the GW amplitude.  \nWe present in Fig. 8 the dependence with the sourcing time ˜ τ sw ≡ ˜ t -˜ t init of the numerical integrated GW amplitude I int sim , compared to the one computed analytically for the power-law fit K ( ˜ t ) = K 0 ( ˜ t/ ˜ t 0 ) -b using the scaling in Eq. (2.19) to extend the results to times after the end of the simulations, and using the decay rate b found in Fig. 5 for each PT. Furthermore, we assume that the GW production starts at the time when the PT is completed, ˜ t 0 ≃ 10, instead of at the starting time of the numerical GW evaluation 13 at ˜ t init = 16, and that the scaling of Eq. (2.19) can also be extended to times ˜ t ∈ [ ˜ t 0 , ˜ t init ]. Then, the integrated K 2 int in this initial time interval based on the power-law fit of K is added to the numerical values in Fig. 8, and we use the extrapolated values K 0 presented  \nFigure 8 . Dependence of the GW amplitude with the source duration ˜ τ sw ≡ ˜ τ -˜ τ ∗ , assuming that the GW production starts at ˜ τ ∗ = ˜ t 0 , normalized by the reference value ˜ Ω GW ≃ 10 -2 [52] and K 2 ξ R ∗ β , for wall velocities v w = 0 . 4, 0 . 56, and 0.8. Red segments indicate simulation results from Figure 6 re-scaled by the extrapolated values of the kinetic energy fraction at ˜ t 0 , K 2 0 , found in Sec. 4.1 (see Fig. 3), and adding the expected contribution to the GW production from the interval of time ˜ t ∈ [ ˜ t 0 , ˜ t init ]. Solid lines indicate the expected amplitudes extrapolated to times after the end of our simulations when the Universe expansion is ignored ( β/H ∗ → ∞ ), while dashed lines correspond to the expected amplitudes after including expansion for β/H ∗ = 1000 and 100 in decreasing opacity. Dots indicate the eddy turnover time ˜ t eddy = βR ∗ / √ K ξ , which determines the expected scale for non-linearities to develop (they do not appear in the plot for weak PTs with v w = 0 . 4 and 0 . 8).  \n<!-- image -->  \nin Sec. 4.1 to estimate the corrections due to the underresolution of the self-similar profiles. To estimate the effect of the Universe expansion on the GW amplitude, we also include the proxy presented in Sec. 2.6 for the values β/H ∗ = 100 and 1000 [see Eq. (2.22)].", '4.4 Gravitational wave spectral shape': 'In this section, we present the numerical results concerning the spectral shape for weak, intermediate, and strong transitions and a range of wall velocities. We present fits to the data and extract spectral features. Results for weak and intermediate transitions were previously obtained in hybrid simulations in Ref. [56] and Higgsless simulations in Ref. [58]. Utilizing the improved Higgsless code, we update the results of Ref. [58] and present new results for strong transitions. In addition to updating the results, we present scaling relations derived from normalizing to R ∗ rather than β , evidently revealing a better scaling behavior of the knee position in the spectrum associated with the typical bubble size.  \nThe findings in Ref. [58] indicate that the GW spectrum Ω GW ( k ) is characterized by a double broken power law: at small k , a Ω GW ( k ) ∝ k 3 scaling was observed, which is also expected from causality. At large k , the spectrum decays as Ω GW ( k ) ∝ k -3 . These scalings are in agreement with those found in the sound-shell model [59, 60]. At intermediate scales, a linear scaling regime Ω GW ( k ) ∝ k was observed. Due to limited resolution, the spectrum appears to exponentially decay beyond a damping scale k e as a result of numerical viscosity. At scales around or beyond the Nyquist wave number, Ω GW ( k ) behaves erratically and is always neglected for the purpose of analysis and parameter extraction.  \nTo capture the behavior of the spectral shape S ( ˜ k ) ≡ I sim ( ˜ k ) / I int sim , we use the following double-broken power law function,  \nS ( k, k 1 , k 2 , k e ) = S 0 ( k k 1 ) n 1 [ 1 + ( k k 1 ) a 1 ] -n 1 + n 2 a 1 [ 1 + ( k k 2 ) a 2 ] -n 2 + n 3 a 2 e -( k/k e ) 2 , (4.10)  \nwhich corresponds to the shape function used in Ref. [92] with an additional exponential damping factor effective above the damping scale k > k e . We expect that the exponential damping found in the simulated spectra are purely due to numerical viscosity so we disregard the parts of the spectra where the exponential damping is relevant. Assuming k 1 < k 2 and k < k e , the fitting parameters correspond to the slopes n 1 , n 2 , and n 3 , such that S ( k ) ∼ k n 1 at small wave numbers k < k 1 , S ( k ) ∼ k n 2 at intermediate k 1 < k < k 2 , and S ( k ) ∼ k n 3 at large k > k 2 . The parameters a 1 and a 2 allow to control the sharpness/smoothness of the spectral shape around the knee and peak at k 1 and k 2 . S 0 is a normalization constant defined by the condition that ∫ S dln k = 1. We note that the choice a 1 = 2, a 2 = 4, n 1 = 3, n 2 = 1, and n 3 = -3 renders Eq. 4.10 equivalent to  \nS f ( k, k 1 , k 2 , k e ) = S 0 × ( k/k 1 ) 3 1 + ( k/k 1 ) 2 [ 1 + ( k/k 2 ) 4 ] × e -( k/k e ) 2 , (4.11)  \nwhich was previously used in Ref. [56]. Equation (4.10), however, allows for a more adaptable recovery of the GW spectrum peak position and slopes by adapting the sharpness/smoothness of the spectral shape around the knee and the peak to the one found in the numerical data.  \nWe expect the characteristic knee and peak of the GW spectra to be determined by the scale of the fluid perturbations R ∗ . Another important length scale is the fluid shell thickness  \nξ shell := ξ front -ξ rear , (4.12)  \nwhere  \nand  \nξ front = { ξ shock , for deflagrations and hybrids , v w , for detonations , (4.13)  \nξ rear = { v w , for deflagrations , c s , for detonations and hybrids . (4.14)  \nThe scale R ∗ ξ shell is expected to determine the peak of the GW spectrum [22, 40, 50, 52, 54, 59, 97].', 'Fitting to the numerical data': 'We fit Eq. (4.10) to our numerically computed GW spectra and thus extract spectral features from our data. We show in Fig. 9 the numerical GW spectra I sim ( ˜ t init , ˜ t end , ˜ k ) found in the simulations with numerical resolution N = 512 and box sizes ˜ L/v w = 20 and 40, for a range of wall velocities, and for weak, intermediate, and strong PTs, together with the analytical fits. We use ˜ t init = 16 and ˜ t end = 32 to evaluate the GW spectra. In the fitting procedure, we impose the constraint that k 1 < k 2 . However, since k e does not represent a physical scale, we do not require that k 2 < k e , but allow k e to take on any value independently. In the cases where k 2 > k e , the spectral peaks are not resolved properly and suffer from numerical viscosity.  \nIn obtaining the fit, we neglect the first bin for simulations with ˜ L/v w = 20 and the first two bins for ˜ L/v w = 40 to avoid the associated significant statistical scatter. We cut the spectra in the UV where the fit including the exponential damping deviates from the broken power law with no exponential damping. While in Fig. 9 we show fits of Eq. (4.10) to the spectra for different v w , in Fig. 10 we show the fitted spectral features k 1 , k 2 , and  \nFigure 9 . Fits of Eq. (4.10) to the numerical results from weak, intermediate, and strong PTs (in each plot, amplitudes increase with larger α ) with N = 512 for a range of v w , and for ˜ L/v w = 20 in brighter colors (white dots for the numerical data), and ˜ L/v w = 40 in darker colors (black dots for the numerical data). Red lines indicate wave numbers below the knee k 1 , green indicates intermediate wave numbers k 1 < k < k 2 , and blue corresponds to wave numbers above the peak k 2 . The dotted orange lines indicate wave numbers k > k e , where exponential damping dominates. The light and dark gray lines indicate the resulting fitted double-broken power laws excluding the exponential damping. Vertical lines indicate the Nyquist wave numbers k Nyq R ∗ = βR ∗ N/ ˜ L .  \n<!-- image -->  \n∗  \n∗  \nFigure 10 . Fitted characteristic wave numbers k 1 (left column), k 2 and k e (middle column), and k peak (right column) for weak (blue), intermediate (red), and strong (orange) PTs, using simulations with N = 512 and ˜ L/v w = 20 (40) in solid (dotted) lines. Gray regions indicate the Nyquist frequency k Nyq = N/ ˜ L . In the upper panel, wave numbers are normalized as ˜ k ≡ k/β , as presented in Ref. [56], while in the lower panel, they are normalized as k R ∗ . Thick colored lines of low opacity in middle column indicate k e for ˜ L/v w = 20 (40) in solid (dotted) lines. In the upper panel (middle and right columns), thin black-and-color dashed lines indicate 1 /ξ shell , while in the lower panel, they indicate the fitted value ∼ 0 . 5 × 2 π/ ∆ w [see Eq. (4.15)]. In the right column, the lower opacity regions indicate the peak as obtained using the double broken power law fit of Eq. (4.10), neglecting the exponential numerical damping.  \n<!-- image -->  \nk e as functions of v w for weak ( α = 0 . 0046), intermediate ( α = 0 . 05), and strong ( α = 0 . 5) PTs. We present these characteristic wave numbers in units of 1 /β and 1 /R ∗ to evaluate the resulting dependence on v w and determine the scale characterizing the spectral knee and peak. We note that the maximum value of the spectral shape used in Eq. (4.10) is located at k peak , which does not in general exactly coincide with k 2 (see discussion in Ref. [92]) and their relation depends on the fitting parameters. We show the resulting spectral peaks obtained from the fit in the right columns of Fig. 10.  \nExtraction of the scale of exponential damping k e gives us a handle on the reliability of the measurement of other parameters and the peak; clearly, finding k 2 > k e means we are  \nin a regime where damping already dominates on scales larger than the peak in the spectrum. In this case, even though for weak transitions k 2 is found to track 1 /ξ shell well above k 2 > k e (which means that we are potentially recovering a trend expected from physical considerations), caution should be taken in interpreting k 2 and k peak as true physical parameters. However, for intermediate and strong PTs, we do not find any evidence for k 2 to be determined by ξ shell , as previously pointed out in Ref. [58]. Using our numerical results with ˜ L/v w = 20, which present better resolution in the UV, averaged over v w and 10 nucleation histories, we find the following values for k 2 ,  \nk 2 R ∗ 2 π ≃      (0 . 49 ± 0 . 024) / ∆ w , α = 0 . 0046 , 0 . 93 ± 0 . 13 , α = 0 . 05 , 0 . 45 ± 0 . 042 , α = 0 . 5 , (4.15)  \nk 1 R ∗ 2 π ≃ 0 . 39 ± 0 . 1 . (4.16)  \nwhere the indicated uncertainty corresponds to the standard deviation in the measurements among the reference simulations, and ∆ w = ξ shell / max( v w , c s ) is the normalized sound-shell thickness. Sample variance is generally of the order of the scatter with wall velocity. On the other hand, the numerical values at the knee, expected to be related to the fluid perturbations scale R ∗ , are found to be  \nWe note that both scales k 1 R ∗ and k 2 R ∗ (for intermediate and strong transitions) present very small variability with v w , indicating a rather universal behavior. For weak PTs, k 2 R ∗ ∆ w ∼ ˜ k 2 is also almost independent of v w , as expected from the sound-shell model. The values of k 2 R ∗ and k 1 R ∗ used in Ref. [92] are based on the numerical results of Ref. [58]. For weak PTs, we find k 2 R ∗ ∆ w consistent with the values used in Ref. [92], while we find the value of k 1 R ∗ to be twice the one used in Ref. [92]. We note that the extraction of the knee k 1 in the IR part of our spectra is more sensitive to statistical variance, underresolution, and the end time of our simulations.', 'Time evolution of the spectral shape in the simulations': 'We show in Fig. 11 the GW spectrum I sim ( ˜ t init , ˜ t, ˜ k ) for different values of ˜ t of the simulation. We find in general that the causal tail, proportional to k 3 at small k , is initially present from early times and a more complex structure seems to develop below the peak as time advances, potentially consistent with analytical work [59] and numerical simulations [60]. We find that the growth of the GW amplitude with the source duration is faster than linear at small wave numbers, which could be described by the quadratic growth found in Ref. [59], where it is shown that the transition from a quadratic to a linear growth with ˜ t (when no significant decay of K occurs) happens at later times for smaller k . The resulting spectral shape at the end of the simulation, ˜ t end = 32, is then shown in Fig. 9 and used to provide fits of the spectral shape. As discussed in Sec. 2.4, the modelling presented and validated for the integrated GW spectrum is only expected to hold at wave numbers kR ∗ ≫ ( βR ∗ ) / ( ˜ t -˜ t init ). Then, we note that using the GW spectrum shape as measured at ˜ t end as the one that ultimately enters the proposed model in Eq. (2.19) effectively implies the assumption that all wave numbers evolve with the source duration in the same way as the overall amplitude until ˜ t fin > ˜ t end is reached. However, different time evolutions than those validated for the integrated amplitude at wave numbers that do not significantly contribute to the integrated amplitude and/or at times after the end of the simulation, could potentially affect the resulting spectral shape of the  \nFigure 11 . Time evolution of the GW spectral shape I ( ˜ t init , ˜ t, ˜ k ) evaluated at times ˜ t ∈ [17 , 32] with ˜ t init = 16 for weak (left columns), intermediate (middle columns), and strong (right columns) PTs, for v w = 0 . 32 (0.36 for strong PTs) (upper panel), 0.6 (middle panel), and 0.8 (lower panel). The numerical resolution is N = 512 and the size of the simulation box is ˜ L/v w = 40. The GW spectra are normalized by the reference value ˜ Ω GW ≃ 10 -2 and the expected scaling K 2 ξ R ∗ β .  \n<!-- image -->  \n∗  \n∗  \n∗  \nGWs. This can occur within the sound-shell model in the IR regime, as shown in Ref. [59], where a transition from the linear towards a quadratic growth is expected at small k in the stationary case, as well as when the sound-shell model (or its generalization in Sec. 2.5 to decaying sources) is no longer valid due to, for example, the potential development of nonlinear fluid perturbations and vortical motion. In the latter case, the resulting GW spectrum is expected to have a different time evolution than the one expected for compressional motion [55, 62, 63, 65, 70, 71] and we expect that the GW modes would reach their saturation amplitudes in this regime.  \nFigure 12 . Fitted UV index n 3 ≥ -3 for the spectral shape of weak (blue), intermediate (red), and strong (orange) PTs, found for simulations with N = 512 and ˜ L/v w = 20. We note that n 3 = -3 is fixed for weak PTs. The error bars show the standard deviation from 10 different bubble nucleation histories. Black vertical line indicates the sound speed c s , while colored lines the Chapman-Jouguet v CJ for weak and intermediate PTs ( v CJ ≃ 0 . 89 > 0 . 8 for strong PTs).  \n<!-- image -->', 'GW spectral slopes': 'In general, we find a clear n 1 = 3 slope at the smallest frequencies, consistent with the expected causal tail, S ( k ) ∼ k 3 . At intermediate wave numbers, we fix n 2 = 1, although this range of k is not large enough to have a clear prediction of the exact intermediate slope. However, it is clear that a smoothing with respect to the k 3 occurs in this range that eventually leads to the decrease k n 3 with n 3 < 0 at large wave numbers k > k 2 . In this regime, the sound-shell model predicts a slope n 3 = -3 [40, 50-52, 56, 58-60], and our simulations show a clear n 3 ≈ -3 whenever k 2 ≪ k e . However, we allow n 3 to be a parameter in our fits to allow for deviations, potentially due to the development of non-linearities. Looking at Fig. 9, it is apparent that generally for strong PTs, the simulations offer sufficient dynamical range to sample the UV slope of the GW spectrum. This is particularly interesting, since for strong PTs, we expect a departure from n 3 = -3 if non-linearities lead to a cascade of energy into the UV, thus modifying the slope towards a Kolmogorov turbulence spectrum with n 3 = -8 / 3 [53, 55, 62, 63, 65, 71, 95], or a shallower acoustic turbulence spectrum [61, 98]. Hence, we allow n 3 ≥ -3 when deriving the fit for strong and intermediate PTs, while restricting n 3 = -3 for weak PTs since the dynamical depth is typically insufficient to recover the UV behavior in these cases, due to thinner shells and hence k 2 ≳ k e .  \nIn Fig. 12, we plot the fitted values of n 3 . For intermediate transitions, we observe a marginal increase in n 3 towards -2 . 5 as the wall velocity is increased. Strong transitions exhibit a similar trend, while also preferring an optimal n 3 ≲ -2 . 75 for small v w .', 'Smoothing/sharpening of the knees': 'Introducing two new free parameters a 1 and a 2 obviously improves the fits to the numerical data, compared to using the simpler Eq. (4.11), and explicitly shows large degeneracies among the fitting parameters, thereby thwarting meaningful interpretation and extraction of the relevant spectral features.  \nSince the peak of the GW spectrum k peak is of greatest phenomenological interest, we adjust the parameters a 1 and a 2 to constants that universally recover the peak position well for all wall velocities and strengths. Empirically, we find that a slight sharpening of the  \nknee and a slight smoothing of the peak typically improves the peak position recovery and yields good results for the fit overall. Measurements of a 1 benefit from simulations with more data points in the IR, and we use exclusively simulations with ˜ L/v w = 40 for its estimation, whereby a 1 = 3 . 6 (i.e., an increase from a 1 = 2 as used in Ref. [56]) is found suitable. Measurements of a 2 , on the other hand, benefit from resolving the UV, for which we use exclusively simulations with ˜ L/v w = 20, and find that a 2 = 2 . 4 (i.e., a reduction from a 2 = 4 as used in Ref. [56]) is an adequate choice. We use these values for a 1 and a 2 throughout this study, but point out that in principle, the spectral fit could be improved by varying these parameters at the cost of a larger scatter in the parameter extraction (due to degeneracies).  \nThe extraction of the parameters k 1 , k 2 , k e , and k peak are shown in Fig. 10. In the different rows, these parameters are expressed in terms of the physical length scales discussed in the last section. Some example spectra are shown in Fig. 9. The slope of the UV tail is plotted in Fig. 12.', '5 Summary and conclusions': 'We have conducted numerical simulations of cosmological first-order phase transitions (PTs) using the Higgsless approach [58] to compute the fluid perturbations in the primordial plasma induced by a PT and the resulting GW spectra, for a range of PT parameters: α = 0 . 0046 (weak), 0 . 05 (intermediate), and 0 . 5 (strong); and a broad range of wall velocities v w ∈ (0 . 32 , 0 . 8). These results extend the previous numerical results of Ref. [58] to strong PTs, and include a larger number of numerical simulations for weak and intermediate PTs.  \nWe present for the first time results of the GW amplitude and spectral shape sourced by fluid perturbations from strong PTs with α = 0 . 5. We have slightly updated the numerical code, although with no significant impact on the numerical results. We have compared our results to those expected considering a stationary unequal time correlator (UETC), an assumption usually made in analytical computations of compressional motion (e.g., sound waves in the limit of linear perturbations) and commonly used to extrapolate the results from numerical simulations, based on the hypothesis that the GWs are produced by a stationary superposition of sound waves. We find strong numerical evidence for the decay of the kinetic energy fraction K with time for intermediate PTs with highly confined profiles and for strong PTs, and a clear deviation with respect to the linear growth of the GW amplitude with the source duration found in previous numerical simulations and analytical studies, assumed in the GW templates used in the literature. We associate this deviation to the decay of the kinetic energy fraction K and extend the stationary UETC modelling to a locally stationary UETC that allows us to introduce the effect of the numerically found decay rate of K with time.  \nThe numerical results presented in this work have allowed us to generalize the usual stationary UETC assumption to an assumption of a locally-stationary UETC, and to test the validity of this assumption in predicting the integrated GW spectral amplitude. Furthermore, the proposed model has allowed us to numerically find the relevant scales that enter in the resulting GW amplitude and spectral shape; see Eq. (2.19). We also have shown that the GW production might not abruptly stop at the time when non-linearities develop but it might keep increasing for a duration that is uncertain at the moment. For this reason, we present our results as a function of the GW source duration. It is of paramount importance to determine its exact value and, hence, the resulting saturated GW amplitude, to make accurate predictions of the GW spectra expected from first-order phase transitions. The  \nmodelling of the remaining stage will also require numerical simulations (as well as the stage over which we find decay of the kinetic energy in our simulations), as it is deep in the non-linear regime and potentially dominated by turbulence.  \nIn the following, we summarize our numerical results by providing a template that can be used by the community to estimate the GW amplitude from first-order phase transitions, validated for the duration of our simulations and extrapolated to later times, taking into account that some of the values presented might be sensitive to numerical uncertainty. Based on the model presented in Sec. 2.5 and validated with our numerical simulations in Sec. 4.3, we find the following parameterization of the GW spectrum when the Universe expansion can be ignored  \nΩ GW ( k ) = 3 T GW ˜ Ω GW ( H ∗ /β ) 2 K 2 int R ∗ β S ( kR ∗ ) , (5.1)  \nwhere S ( k ) denotes the shape function of the spectrum that is normalized to ∫ d ln k S ( k ) = 1, and K 2 int is the integrated kinetic energy fraction K 2 over ˜ t ≡ tβ , such that it reduces to K 2 τ sw β when K is constant, being τ sw the GW source duration. Therefore, Eq. (5.1) is a generalization of the parameterization used in the stationary UETC assumption previously tested with numerical simulations [40, 50, 52] and usually assumed for sound-wave sourcing of GWs [22, 51, 54, 59, 91, 92, 99] that predicts a linear growth with the GW source duration when K does not decay with time.  \nThe most robust results (i.e., an almost independent value of ˜ Ω GW with the PT parameters) are obtained when the typical bubble separation R ∗ , which determines the length scale of fluid perturbations, is given by the front of the expanding bubbles [22]  \nβR ∗ = (8 π ) 1 / 3 max(v w , c s ) , (5.2)  \nwhere 1 /β parameterizes the duration of the PT, v w is the wall velocity, and c s the speed of sound. This way, the residual dependence on the wall velocity in ˜ Ω GW is quite limited and we estimate from our numerical simulations values for the GW efficiency ˜ Ω GW ∼ O (10 -2 ) for a range of PTs [see Fig. 7 and Eq. (4.9)],  \n10 2 ˜ Ω GW =      1 . 04 +0 . 81 -0 . 67 , for α = 0 . 0046 ; 1 . 64 +0 . 29 -0 . 13 , for α = 0 . 05 ; 3 . 11 +0 . 25 -0 . 19 , for α = 0 . 5 , (5.3)  \nconsistent with previous numerical simulations [40, 50, 52] for weak and intermediate PTs, and with the sound-shell model [51, 54] for weak PTs. For intermediate and strong PTs, we find much less dependence with v w than for weak PTs, clearly showing a departure with respect to the predictions of the sound-shell model (see Fig. 7).  \nWe also provide an estimate of the relevant kinetic energy fraction K 0 at the end of the PT using our numerical results [see Fig. 3 and Eq. (4.4)], given in units of the single-bubble K ξ [see Eq. (2.17)], which, averaged over wall velocities, becomes  \nK 0 = 0 . 84 +0 . 24 -0 . 29 K ξ . (5.4)  \nAs a function of v w , we generally find that K ξ might slightly underestimate K 0 for the smallest v w , while it tends to overestimate it for larger v w (see Fig. 3). This might be a consequence of the energy transfer between thermal and kinetic energies during the phase of collisions and the expected development of the sound-wave regime [51, 54, 59].  \nWe have studied the decay of the kinetic energy fraction K with time ˜ t in Sec. 4.2, and provide a power-law fit K ( ˜ t ) = K 0 ( ˜ t/ ˜ t 0 ) -b , with b > 0 indicating the decay rate, that accurately reproduces the numerical results (see Figs. 4 and 5). For small values of b (hence, negligible decay), one can directly use  \nK 2 int ( b = 0) →K 2 0 ˜ τ sw →K 3 / 2 0 βR ∗ (5.5)  \nin Eq. (5.1), assuming that the duration of the GW sourcing is given by the eddy turnover time ˜ τ sw ∼ ˜ t eddy = βR ∗ / √ K , when non-linearities are expected to develop. In general, we find b ≪ 1 when the eddy turnover time is larger than our final simulation time ˜ t eddy ≫ ˜ t end = 32 (for weak PTs and some intermediate ones).  \nFor non-negligible values of b , we find that the decay of K occurs within the duration of our simulations, potentially indicating that we are already modelling the GW production in the non-linear regime. We indeed find that this might be the case as the eddy turnover time is included in the duration of our simulations for some intermediate PTs and for strong ones, where we find larger values of b . For these PTs, we find that the integrated K 2 int becomes  \nK 2 int ( b, τ sw ) →K 2 0 β t ∗ (1 + τ sw /t ∗ ) 1 -2 b -1 1 -2 b , (5.6)  \nwhen one uses the power-law fit for K ( ˜ t ) and assumes that the GW production roughly starts at the time ˜ t ∗ ≃ ˜ t 0 ≃ 10 (note that the actual value of ˜ t 0 only appears as a consequence of our particular fit). It is unclear what should be the final time of GW sourcing in these cases, as the simulations seem to already be modelling the non-linear regime, so we leave ˜ τ sw as a free parameter. We note that this is an indication that the GW spectrum might still grow once that non-linearities develop in the fluid, such that the use of Eq. (5.5) would in general underestimate the GW production. We compare in Fig. 8 the numerical dependence of the GW amplitude with the source duration ˜ τ sw found in the simulations to the one obtained using Eq. (5.6), extending the analytical fit beyond the time when the simulations end.  \nAs a final remark on the integrated GW amplitude, we note that so far Universe expansion has been ignored, which is not justified for long source durations. Taking into account that the fluid equations are conformal invariant after the PT if the fluid is radiationdominated, we can apply the results from our fluid simulations in Minkowski space-time to an expanding Universe, as long as the PT duration is short ( β/H ∗ ≫ 1) even if the GW source duration is not short (see discussion in Sec. 2.6). Then, as a proxy to estimate the effect of the Universe expansion, we can use the following value for K 2 int [see Eq. (2.23)]  \nK 2 int →K 2 0 Υ b ( τ sw ) ( β/H ∗ ) , (5.7)  \nwhich generalizes the suppression factor Υ = H ∗ τ sw / (1 + H ∗ τ sw ) when the source does not decay [89, 91] to any decay rate b using Eq. (2.24) for the presented power-law decay fit of K ( ˜ t ). We also compare in Fig. 8 the expected evolution of the GW amplitude with the source duration according to Eq. (5.7) for β/H ∗ = 100 and 1000. We note that when one associates τ sw to the eddy turnover time τ eddy = R ∗ / √ K , they should correspond to conformal time intervals, instead of cosmic time, due to the conformal invariance of the fluid equations.  \nRegarding the spectral shape S ( kR ∗ ) in Eq. (5.1), we find that the following template fits accurately our numerical results (see Sec. 4.4)  \nS ( k, k 1 , k 2 ) = S 0 × ( k k 1 ) n 1 [ 1 + ( k k 1 ) a 1 ] -n 1 + n 2 a 1 [ 1 + ( k k 2 ) a 2 ] -n 2 + n 3 a 2 , (5.8)  \nwith n 1 ≃ 3, n 2 ≃ 1, a 1 ≃ 3 . 6, and a 2 ≃ 2 . 4. We note that to compare with our numerical results we have included an exponential damping e -( k/k e ) 2 [see Eq. (4.10)], effective at k > k e , but we omit it here as we expect it to correspond to numerical viscosity and not have physical relevance. The slope of the UV tail is n 3 ≃ -3 for weak PTs, and intermediate ones with small wall velocities v w ≲ c s . The slope becomes slightly shallower (up to -2 . 5) for intermediate PTs with supersonic v w and strong PTs (see Fig. 12). This effect should not play a major role in phenomenological studies but a more detailed description is given in Sec. 4.4. Furthermore, this shallower GW spectral slope, together with the decay of the kinetic energy, seems to indicate the development of non-linearities. To confirm this statement would require a detailed study of the kinetic spectrum properties. For now, we present a preliminary study of the vorticity production in our simulations in App. C.  \nThe most relevant feature of the spectrum is the position of the peak, determined by k 2 . Here, we find a distinction between weak and intermediate/strong PTs (as already previously seen in Ref. [58]). For weak PTs, the peak follows the thickness of the fluid shells, ξ shell , as given in Eq. (4.12), while for intermediate and strong PTs the dependence on the wall velocity is much weaker. As shown in Eq. (4.15) and Fig. 10, we find the following results for k 2 , averaged over all wave numbers and 10 different nucleation histories: k 2 R ∗ ≃ π/ ∆ w for weak PTs, k 2 R ∗ ≃ 2 π for intermediate PTs, and k 2 R ∗ ≃ π for strong ones, where ∆ w = ξ shell / max( v w , c s ) is the normalized sound-shell thickness, which is only found to determine k 2 in weak PTs. Finally, the position of the knee that relates to the typical size of the bubbles does not even depend on the strength of the PT and, quite generally, we find k 1 R ∗ ≃ 0 . 4 × 2 π (see Eq. (4.16) and Fig. 10).', 'Acknowledgments': "This research was supported in part through the Maxwell computational resources operated at Deutsches Elektronen-Synchrotron DESY, Hamburg, Germany. TK and IS acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2121 'Quantum Universe' - 390833306. HR is supported by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy EXC 2094 'ORIGINS'. (No.390783311). ARP acknowledges support by the Swiss National Science Foundation (SNSF Ambizione grant 182044). ARP and CC acknowledge the working space provided during the program on the 'Generation, evolution, and observations of cosmological magnetic fields' at the Bernoulli Center in Lausanne. RJ is supported by JSPS KAKENHI Grant Numbers 23K17687, 23K19048, and 24K07013.", 'A Gravitational wave production from early Universe sources': "The solution to the GW equation for the tensor-mode perturbations h ij , defined such that the line element is  \nwhile the source is active 14 ( t < t fin ) is [17]  \nh ij ( t < t fin , k ) = 6 H 2 ∗ k ∫ t t ∗ Π ij ( t ' , k ) sin k ( t -t ' ) d t ' , (A.2)  \nd s 2 = a 2 [ -d τ 2 +( δ ij + h ij ) d x i d x j ] , (A.1)  \nwhere t ∗ is the initial time at which the tensor of anisotropic stresses, Π ij = Λ ijlm T lm / ¯ ρ , starts to source GWs, being ¯ ρ = 3 H 2 M 2 Pl the critical energy density and Λ ijlm = P il P jm -1 2 P ij P lm the traceless and transverse projector, with P ij = δ ij -ˆ k i ˆ k j . We have assumed that the expansion of the Universe is negligible during the sourcing process, i.e., t fin -t ∗ ≪ H -1 ∗ .  \nAt later times, after the sourcing has ended at t fin , the solution is  \nh ij ( t ≥ t fin , k ) = 6 H 2 ∗ k ∫ t fin t ∗ Π ij ( t ' , k ) sin k ( t -t ' ) d t ' . (A.3)  \nThen, the time derivatives of the strains h ij are  \n∂ t h ij ( t ≥ t fin , k ) = 6 H 2 ∗ ∫ t fin t ∗ Π ij ( t ' , k ) cos k ( t -t ' ) d t ' , (A.4)  \nwhich can be used to find the fractional energy density at present time t 0 [17],  \nΩ GW ( t 0 ) = ρ GW ρ 0 tot = ( a ∗ /a 0 ) 4 12 H 2 0 ⟨ ∂ t h ij ( t 0 , x ) ∂ t h ij ( t 0 , x ) ⟩ . (A.5)  \nWe consider the GW spectrum Ω GW ( t 0 , k ) ≡ dΩ GW ( t 0 ) / dln k , which describes the twopoint correlation function of the statistically homogeneous and isotropic strain derivatives, following the notation of Ref. [66]  \n( a ∗ /a 0 ) 4 12 H 2 0 ⟨ ∂ t h ij ( t 0 , k ) ∂ t h ∗ ij ( t 0 , k ' ) ⟩ = (2 π ) 6 δ 3 ( k -k ' ) Ω GW ( t 0 , k ) 4 πk 3 , (A.6)  \nsuch that Ω GW ( t 0 ) = ∫ Ω GW ( t 0 , k ) d ln k . Substituting Eq. (A.4) into Eq. (A.6), we find  \nΩ GW ( t 0 , k ) = 3 k T GW H 2 ∗ ∫ t fin t ∗ ∫ t fin t ∗ E Π ( t 1 , t 2 , k ) cos k ( t 0 -t 1 ) cos k ( t 0 -t 2 ) d t 1 d t 2 , (A.7)  \nwhere T GW ≡ ( a ∗ /a 0 ) 4 ( H ∗ /H 0 ) 2 is the transfer function and E Π ( t 1 , t 2 , k ) is the unequal-time correlator (UETC) of the anisotropic stresses,  \n⟨ Π ij ( t 1 , k ) Π ∗ ij ( t 2 , k ' ) ⟩ = (2 π ) 6 δ 3 ( k -k ' ) E Π ( t 1 , t 2 , k ) 4 πk 2 . (A.8)  \nAt present time, for modes kt 0 ≫ 1, we can average the product of Green's functions in Eq. (A.7) over oscillations to find  \nΩ GW ( k ) = 3 k 2 T GW H 2 ∗ ∫ t fin t ∗ ∫ t fin t ∗ E Π ( t 1 , t 2 , k ) cos k ( t 1 -t 2 ) d t 1 d t 2 . (A.9)  \nTherefore, once we know the UETC of the source of GWs, E Π ( t 1 , t 2 , k ), in this case produced by the fluid perturbations Π ij = wγ 2 Λ ijlm v l v m / ¯ ρ , we can directly compute the GW spectrum. However, this is in general unknown before the simulation and, hence, we need to estimate it using the numerical results. For this purpose, we can compute the UETC by approximating the ensemble average of the anisotropic stresses in Eq. (A.8) with the average over spherical shells of radius k in Fourier space,  \nE Π ( t 1 , t 2 , k ) = k 2 2 π 2 V ∫ Ω k dΩ k 4 π Π ij ( t 1 , k ) Π ∗ ij ( t 2 , k ) . (A.10)  \nThen, substituting Eq. (A.10) into Eq. (A.9) and taking into account that, since cos k ( t 1 -t 2 ) = cos kt 1 cos kt 2 +sin kt 1 sin kt 2 , the double integral over t 1 and t 2 can be expressed as the following product  \nΩ GW ( k ) = 3 k 3 4 π 2 V T GW H 2 ∗ ∫ Ω k dΩ k 4 π ∫ t fin t ∗ Π ij ( t 1 , k ) e ikt 1 d t 1 × ∫ t fin t ∗ Π ∗ ij ( t 2 , k ) e -ikt 2 d t 2 , (A.11)  \nwhere we have used the fact that the resulting Ω GW ( k ) is real. Finally, defining the following integral over the stress-energy tensor T ij = wγ 2 v i v j ,  \nT ij ( q, k ) = ∫ t fin t ∗ T ij ( t, k ) e iqt d t , (A.12)  \nand using the property Π ij Π ∗ ij = Λ ijlm T ij T ∗ lm / ¯ ρ 2 , we find the following expression  \nΩ GW ( k ) = 3 k 3 4 π 2 V ¯ ρ 2 T GW H 2 ∗ ∫ Ω k dΩ k 4 π Λ ijlm ( ˆ k ) [ T ij ( q, k ) T ∗ lm ( q, k ) ] q = k , (A.13)  \npreviously used in Refs. [56, 58] and referred to as Weinberg's formula, due to its similarity to the expression obtained for the power emitted by isolated deterministic binaries [100]. However, we note that the applicability of Eq. (A.13) is limited to the knowledge of the stochastic variables T ij , which is generated by the simulations from a single realization, and results from the average over shells under the assumption of statistically homogeneity and isotropy. This expression allows us to compute the GW spectrum at present time after averaging over oscillations in time. However, we compute numerically this expression until the end of the simulation at t end , which is, in general, smaller than t fin . Therefore, the numerical result will represent the physical amplitude at a mode k if it has already reached its freepropagation regime by this time, i.e., if its amplitude has already saturated. Otherwise, the numerical result presents the GW amplitude that would be obtained misleadingly assuming that the source is abruptly switched off at t end .", 'B Corrections to the kinetic energy for multiple bubbles': "In this section, we study in more detail the convergence of the numerical simulations of multiple bubbles by comparing their convergence to that of single-bubble simulations, for which we expect the self-similar profiles described in Ref. [49] to develop. For reference, we show in Fig. 13 the self-similar velocity and enthalpy profiles for the PTs considered in the present work (see Tab. 1). We have previously defined the kinetic energy fraction K such that ¯ ρK ( ˜ t ) ≡ ⟨ ρ kin ( x , ˜ t ) ⟩ , where ⟨ ρ kin ⟩ corresponds to the kinetic energy density averaged over the simulation volume V . However, we note that K ξ for single bubbles, defined in Eq. (2.17), is taken as the average of the kinetic energy density fraction induced by a single bubble over the broken-phase volume. Then, defining the ratio of the volume in the broken phase (bp) to the total volume,  \nV ( ˜ t ) = V bp V , (B.1)  \nwe can define the analog of K ξ for multiple bubbles as the ratio K ( ˜ t ) / V ( ˜ t ). Before fluid shells collide, and in the limit of infinite resolution, this ratio should be identical to K ξ after  \nvw  \nFigure 13 . One-dimensional profiles of the fluid velocity (left columns) and enthalpy (right columns) perturbations for a single bubble as a function of ξ ≡ r/ ( t -t n ) for weak (upper panel), intermediate (middle panel), and strong (lower panel) PTs, and for a range of wall velocities. The profiles are computed using CosmoGW [96]. Vertical dashed lines indicate the speed of sound c s and the Chapman-Jouget velocity v CJ [49]. These profiles correspond to the set of simulations presented in this work (see Tab. 1).  \n<!-- image -->  \na very short transient period, over which the fluid profiles develop. Deviations from K ξ before collisions thus correspond to an artifact due to numerical inaccuracy. We plot the ratio K ( ˜ t ) / [ V ( ˜ t ) K ξ ] as solid lines in Fig. 14 for all four resolutions N ∈ { 64 , 128 , 256 , 512 } .  \nIn the multiple-bubble runs, we can similarly define ¯ ρK i ( ˜ t ) = ⟨ ρ kin ,i ( x , ˜ t ) ⟩ as the kinetic energy fraction of each of the single bubbles i before their corresponding first collision, and  \nFigure 14 . Time evolution of the kinetic energy fraction in the broken-phase volume K ( ˜ t ) / V ( ˜ t ) for multiple-bubble simulations (solid lines), normalized by the single-bubble K ξ , for different resolutions N = { 64 , 128 , 256 , 512 } in increased opacity and box size ˜ L/v w = 20. Results are shown for weak (left columns), intermediate (middle columns), and strong (right columns) PTs, and for a range of wall velocities v w = { 0 . 36 , 0 . 48 , 0 . 6 , 0 . 8 } . Dashed lines correspond to the ratio K Σ ( ˜ t ) / V Σ ( ˜ t ) computed from the single-bubble simulations, such that the departures between the solid and dashed lines indicate the time when fluid-shell collisions take place in the multiple-bubble simulations. Black dots are the values of K 0 obtained from the fit K ( ˜ t ) = K 0 ( ˜ t/ ˜ t 0 ) -b studied in Sec. 4.2 for different N . Red and green dots correspond to the estimated values K 0 [see Eq. (B.4)] and K ∞ 0 (obtained from the convergence analysis of Sec. 4.1). Orange stars correspond to the factor S at the collision time, ˜ t coll , used to compute K 0 from K 0 [see Eq. (B.4)].  \n<!-- image -->  \ndefine the ratio K i ( ˜ t ) / V i ( ˜ t ), where V i corresponds to the fractional broken-phase volume occupied by each bubble i .  \nTo monitor the time-dependence of K i , we simulate single bubbles nucleated at the center of the simulation box (see 'single-bubble' runs in Tab. 1). As the convergence of the single-bubble profiles depends on the resolution in ξ ≡ r/ ( t -t i ), where t i is the nucleation  \ntime of the bubble i and r the radial distance to the nucleation center, we empirically find that doubling the resolution from N to 2 N is equivalent to evaluating the profile at time 2( ˜ t -˜ t i ) (see also discussion in Ref. [58]). Hence, the kinetic energy of single-bubble simulations (with t i = 0) obeys K 2 N i ( ˜ t ) / V i ( ˜ t ) = K N i (2 ˜ t ) / V i ( ˜ t ) to an excellent degree and it suffices to run single-bubble simulations for the largest resolution N = 512. These simulations are run approximately until the front of the fluid profile collides with its own mirror image at the edge of the simulation box, which occurs around ˜ t sb end = ˜ Lv w / 2 / max( c s , v w ). 15  \nThen, in the full simulations and before fluid shells collide, the state of the simulation is exactly the superposition of single bubbles nucleated at times ˜ t i < ˜ t in the bubble nucleation history. We thus construct the sum  \nK Σ ( ˜ t ) ≡ ∑ i : { ˜ t i < ˜ t } K i ( ˜ t -˜ t i ) , (B.2)  \nwhich corresponds to the expected kinetic energy fraction for multiple-bubble simulations in the hypothetical case that no single bubble would collide, following the bubble nucleation history up to time ˜ t . Then, before the first fluid-shell collision occurs, we have that K Σ ( ˜ t ) = K ( ˜ t ), while K ( ˜ t ) starts to deviate from K Σ ( ˜ t ) after the first collision at ˜ t coll .  \nSimilarly, we can construct the fractional broken-phase volume occupied by the superposition of single bubbles as  \nV Σ ( ˜ t ) ≡ ∑ i : { ˜ t i < ˜ t } V i ( ˜ t -˜ t i ) , (B.3)  \nwhich can become larger than one, as it ignores interactions between bubbles. However, the ratio K Σ / V Σ is bounded by K ξ . 16 We plot the time evolution of the ratio K Σ / V Σ as dotted lines in Figure 14 using the nucleation history of the reference multiple-bubble simulations with ˜ L/v w = 20.  \nThe ratio K Σ / V Σ indicates the global degree of convergence of the full multiple-bubble simulations in the hypothetical case that all bubbles keep evolving without interacting with other bubbles. Therefore, the ratio K/ V computed in the multiple-bubble simulations is initially identical to K Σ / V Σ at times ˜ t < ˜ t coll . However, as collisions take place, we clearly see in Fig. 14 that both fractions deviate from each other, as a consequence of mainly four phenomena: (1) the self-similar profiles stop converging towards the expected ones when collisions take place, and since the kinetic energy of the uncollided bubbles is in general underestimated, the saturated value of K ( ˜ t ) / V ( ˜ t ) found quickly after a short transient period dominated by collisions, will be underestimated; (2) oscillatory conversion between thermal and kinetic energy; (3) upon collisions, the fluid self- and inter-shell interactions may be non-linear and result in the kinetic energy decay studied in Sec. 4.2, which is again affected by the previous two effects; and (4) numerical viscosity also results in damping of the kinetic energy. The first and last phenomena are purely numerical, while the remaining two are physical effects that might be affected by the numerical accuracy.  \n16 This is the case for all the considered PTs, with the exception of strong PTs with v w = 0 . 36 and v w = 0 . 4, where values K Σ / V Σ ≳ K ξ are found before collisions due to numerical inaccuracy (see Fig. 14). However, as time evolves the fraction asymptotically tends to K ξ .  \nIn Fig. 14, we mark with orange stars the time of first collision ˜ t coll , where we assume that the maximum degree of convergence of the self-similar profiles is reached, as at later times collisions might affect the development of the fluid-shell profiles. Thus, we can attempt to compensate for the underestimation of the kinetic energy fraction due to insufficient resolution at ˜ t > ˜ t coll multiplying K ( ˜ t ) by the factor S = V ( ˜ t coll ) K ξ /K ( ˜ t coll ), effectively correcting to the expected value K ξ at the time when collisions affect the value of K/ V . In particular, the kinetic energy fraction at the time when the PT ends, ˜ t 0 , can be corrected to the following value  \nK 0 = S K 0 = V ( ˜ t coll ) K ξ K ( ˜ t coll ) K 0 . (B.4)  \nIn Fig. 3, we plot K 0 obtained for numerical resolutions N ∈ { 64 , 128 , 256 , 512 } . Compared to the extrapolated value of K 0 in Fig. 2, we find faster convergence in K 0 when we compare the results for the two largest resolutions N = 256 and N = 512. As we are taking into account the known degree of convergence of the self-similar profiles when computing K 0 , we propose that it is a better estimate of the actual value at ˜ t 0 than K ∞ 0 . Furthermore, the resulting values are closer to K ξ and therefore are more conservative estimates, reducing potential deviations with respect to K ξ that might be a numerical artifact. However, one needs to keep in mind that the underresolution at the time of collisions might strongly affect the posterior evolution of the kinetic energy when non-linearities dominate the dynamics.", 'C Preliminary results for vorticity': 'In this section, we present some preliminary measurements of the fluid vorticity in our numerical simulations.', 'C.1 Vorticity on the lattice': 'The vorticity is computed as  \n∇× v = ( ∂v z ∂y -∂v y ∂z ) ˆ x + ( ∂v x ∂z -∂v z ∂x ) ˆ y + ( ∂v y ∂x -∂v x ∂y ) ˆ z . (C.1)  \nOn the lattice, we approximate the derivatives using first-order central differences,  \n∂v i ∂ ˜ x j (˜ x ) ≃ v i [˜ x + δ ˜ x ˆ x j ] -v i [˜ x -δ ˜ x ˆ x j ] 2 δ ˜ x , (C.2)  \nwhere δ ˜ x is the uniform grid spacing.  \nWith this choice of derivative operator, the magnitude of the curl | ˜ ∇× v | is computed at every grid point. 2D simulation slices at different times are shown in the lower panels of Fig. 1. Note that the definition of the numerical derivative operator of Eq. (C.2) inevitably introduces potentially large vorticity at points where the velocity field varies considerably from lattice site to lattice site. This occurs, e.g., around the bubble shock fronts where discontinuities are present. Ideally, the velocity gradients are aligned with the radial direction, in which case no vertical component is present. However, on the lattice, artifacts may arise from the discretization, causing rather strong vorticity to appear at and just around the bubble walls. This is clearly seen in the lower left frame of Fig. 1. The numerical nature of this vorticity is nevertheless clear from the observation that the vortical structure, as we traverse around bubble wall, is seen to inherit the symmetry of the lattice (see, in particular, the  \nlargest bubble in the center). Furthermore, mostly small but spurious oscillations of the fluid velocity occur at the bubble wall interface. These oscillations additionally give rise to extremely local but very steep velocity gradients, potentially showing up as spurious vorticity with large amplitude, confined to very small scales.  \nIn the lower panel of Figure 1, and in particular when presented with the opportunity to study the full time-evolution of the system frame by frame, it is observed that production of vorticity occurs, upon collisions, at the interface of a sound-shell from one bubble crossing over the bubble wall of another, as seen, e.g., immediately to the right of the top section of the central bubble in the left column. The velocity field in the upper panel is included to make vorticity production easy to correlate with the velocity field. In this sense, the resulting vorticity pattern initially appears to track the sweeping of this sound-shell-bubblewall-crossing interface over time. However, during this process, frame-by-frame inspections indicate that convective non-linear motion is induced in the fluid. This motion implies the presence of slowly evolving structures compared to sound waves propagating at the speed of sound, c s . The former convective structures evolve on time scales proportional to the average convection speeds, which are typically much smaller than c s . It therefore appears that, in addition to a pure longitudinal velocity component, a fluid velocity field characterized by convective motion develops as a result of fluid interactions during, and possibly after, the collision phase. This convective component is marginally evident in the right column of Figure 1, where additional small-scale structures are hinted in the velocity field but are absent in the enthalpy field.  \nWe have described the expected and observed presence of spurious vorticity components associated with the choice of the derivative operator and the lattice structure, and spurious osculations around the bubble wall interface. These are very localized effects and do not contribute meaningfully to large-scale vorticity correlated over macroscopic scales. Therefore, vorticity components that emerge from numerically induced oscillations and limited grid resolution will contribute mostly to the UV part of the velocity spectra. Furthermore, the presence of convective motion implies the presence of a transverse component of the velocity spectra. The development of such a component should thus be visible in spectra of the velocity fields decomposed into longitudinal and transverse contributions, and in particular, physical macroscopic contributions should distinguish themselves from numerical contributions through a separation of scales.', 'C.2 Velocity power spectra': "We generate the Fourier transform v ( k ) of the fluid velocity field and then extract the power spectrum from the two-point correlation in Fourier space [101]  \n⟨ v i ( k ) v ∗ i ( k ' ) ⟩ = (2 π ) 3 δ 3 ( k -k ' ) P v ( k ) , (C.3)  \nwhere the ensemble average is performed over momenta with the same absolute value | k | = | k ' | due to statistical homogeneity and isotropy. Along the same lines we also construct the power spectrum of the longitudinal modes  \n⟨ ˆ k i v i ( k ) ˆ k j v ∗ j ( k ' ) ⟩ = (2 π ) 3 δ 3 ( k -k ' ) P ∥ ( k ) , (C.4)  \nwhere ˆ k = saw( k ) / | saw( k ) | is a unit vector according to the saw description, as discussed in Sec. 3.1 [see Eq. (3.3)]. Likewise, we extract the vortical component of the spectra  \n〈 [ ˆ k × v ( k ) ] i [ ˆ k ' × v ∗ ( k ' ) ] i 〉 = (2 π ) 3 δ 3 ( k -k ' ) P ⊥ ( k ) , (C.5)  \nFigure 15 shows the power in longitudinal and vortical modes as well as the fraction of power in vortical modes for strong PTs. The power in longitudinal modes builds up when the first bubbles nucleate while vorticity requires collisions, as expected. There are some artifacts in the UV even before collisions which result from the discretization of space on a grid. Overall, the power in vorticity can reach values of around P ⊥ /P v ≃ 0 . 3 for the deflagration ( v w = 0 . 44), while this fraction is somewhat smaller, P ⊥ /P v ≃ 0 . 1, for the hybrid ( v w = 0 . 80). Still, a sizable fraction of vorticity is observable in both cases. Moreover, the power spectra in vortical modes appear somewhat less steep than the ones in longitudinal modes. All this is consistent with the hypothesis that the energy loss observed over time for  \n<!-- image -->  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n3  \n5  \n7  \n9  \n11  \n13  \n10  \n10  \n-  \n10  \n-  \n0  \n1  \n2  \n3  \n4  \n)  \n˜  \nk  \n(  \n⊥  \nP  \n)  \n˜  \nk  \n(  \nv  \n/P  \n)  \n˜  \nk  \n(  \n⊥  \nP  \n10  \n-  \n10  \n-  \n10  \n10  \n0  \n0  \nk  \n˜  \nk  \n˜  \nk  \nvw  \n= 0  \n.  \n80  \n˜  \nk  \n˜  \nk  \n˜  \nk  \nFigure 15 . Transverse P ⊥ ( k ) and longitudinal P ∥ ( k ) velocity power spectra (in arbitrary units) of the fluid for a strong PT ( α = 0 . 5) with wall velocities v w = 0 . 44 (left columns) and v w = 0 . 8 (right columns). The panels show the power in longitudinal (top) and vortical (middle) modes as well as the fraction of power in vortical modes (bottom), P ⊥ ( k ) /P v ( k ). The different lines correspond to different times in the simulation ˜ t ∈ (2 , 32). For reference, we remind that the first bubbles nucleate around ˜ t ≃ 0, the first collisions happen around ˜ t ≃ 5, and the simulation is filled with the broken phase at ˜ t 0 ≃ 10.  \nsuch that P v ( k ) = P ⊥ ( k ) + P ∥ ( k ).  \n10  \n10  \n10  \n10  \n10  \n1  \n1  \n10  \n10  \n2  \n2  \n)  \n˜  \nk  \n(  \n‖  \nP  \n)  \n˜  \nk  \n(  \n⊥  \nP  \n)  \n˜  \nk  \n(  \nv  \n/P  \n)  \n˜  \nk  \n(  \n⊥  \nP  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n10  \n-  \n3  \n5  \n7  \n9  \n11  \n13  \n3  \n5  \n7  \n9  \n11  \n13  \n10  \n10  \n-  \n10  \n-  \n0  \n1  \n2  \n3  \n4  \n10  \n10  \n10  \n0  \n0  \n0  \n1  \n1  \n1  \n2  \n/β  \n4  \n/β  \n6  \n/β  \n8  \n/β  \n10  \n/β  \n12  \n/β  \n16  \n/β  \n32  \n/β  \n10  \n10  \n10  \n2  \n2  \n2  \nstrong PTs is due to the decay of fluid kinetic energy into vortical motion and eventually turbulence. This point deserves further attention in the future. We also have obtained velocity power spectra for weak and intermediate PTs, finding that P ⊥ /P v < 10 -3 for those scenarios. This indicates that turbulence becomes progressively more important the stronger the PT, and that for strong PTs ( α = 0 . 5) it already has a dominant role in determining the hydrodynamical evolution after the PT ends. This observation underlines the importance of fully non-linear 3D simulations.", 'References': "- [1] LIGO Scientific, Virgo collaboration, Observation of Gravitational Waves from a Binary Black Hole Merger , Phys. Rev. Lett. 116 (2016) 061102 [ 1602.03837 ].\n- [2] LIGO Scientific, Virgo collaboration, GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence , Phys. Rev. Lett. 116 (2016) 241103 [ 1606.04855 ].\n- [3] LIGO Scientific, Virgo collaboration, GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral , Phys. Rev. Lett. 119 (2017) 161101 [ 1710.05832 ].\n- [4] KAGRA, VIRGO, LIGO Scientific collaboration, GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run , Phys. Rev. X 13 (2023) 041039 [ 2111.03606 ].\n- [5] EPTA, InPTA: collaboration, The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals , Astron. Astrophys. 678 (2023) A50 [ 2306.16214 ].\n- [6] NANOGrav collaboration, The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background , Astrophys. J. Lett. 951 (2023) L8 [ 2306.16213 ].\n- [7] D.J. Reardon et al., Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array , Astrophys. J. Lett. 951 (2023) L6 [ 2306.16215 ].\n- [8] H. Xu et al., Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I , Res. Astron. Astrophys. 23 (2023) 075024 [ 2306.16216 ].\n- [9] EPTA collaboration, The second data release from the European Pulsar Timing Array: IV. Implications for massive black holes, dark matter and the early Universe , 2306.16227 .\n- [10] NANOGrav collaboration, The NANOGrav 15 yr Data Set: Search for Signals from New Physics , Astrophys. J. Lett. 951 (2023) L11 [ 2306.16219 ].\n- [11] LISA collaboration, Laser Interferometer Space Antenna , 1702.00786 .\n- [12] M. Colpi et al., LISA Definition Study Report , 2402.07571 .\n- [13] eLISA collaboration, The Gravitational Universe , 1305.5720 .\n- [14] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender et al., Laser interferometer space antenna , 2017.\n- [15] LISA Cosmology Working Group collaboration, Cosmology with the Laser Interferometer Space Antenna , 2204.05434 .\n- [16] M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics and Cosmology , Oxford University Press (3, 2018).\n- [17] C. Caprini and D.G. Figueroa, Cosmological Backgrounds of Gravitational Waves , Class. Quant. Grav. 35 (2018) 163001 [ 1801.04268 ].  \n- [18] E. Witten, Cosmic Separation of Phases , Phys. Rev. D30 (1984) 272.\n- [19] K. Kajantie, M. Laine, K. Rummukainen and M.E. Shaposhnikov, Is there a hot electroweak phase transition at m H ≳ m W ? , Phys. Rev. Lett. 77 (1996) 2887 [ hep-ph/9605288 ].\n- [20] M. Gurtler, E.-M. Ilgenfritz and A. Schiller, Where the electroweak phase transition ends , Phys. Rev. D 56 (1997) 3888 [ hep-lat/9704013 ].\n- [21] F. Csikor, Z. Fodor and J. Heitger, Endpoint of the hot electroweak phase transition , Phys. Rev. Lett. 82 (1999) 21 [ hep-ph/9809291 ].\n- [22] C. Caprini et al., Detecting gravitational waves from cosmological phase transitions with LISA: an update , JCAP 03 (2020) 024 [ 1910.13125 ].\n- [23] S.R. Coleman, The Fate of the False Vacuum. 1. Semiclassical Theory , Phys. Rev. D 15 (1977) 2929.\n- [24] A.D. Linde, Fate of the False Vacuum at Finite Temperature: Theory and Applications , Phys. Lett. B 100 (1981) 37.\n- [25] P.J. Steinhardt, Relativistic Detonation Waves and Bubble Growth in False Vacuum Decay , Phys. Rev. D 25 (1982) 2074.\n- [26] V.A. Kuzmin, V.A. Rubakov and M.E. Shaposhnikov, On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe , Phys. Lett. B 155 (1985) 36.\n- [27] A.G. Cohen, D.B. Kaplan and A.E. Nelson, Progress in electroweak baryogenesis , Ann. Rev. Nucl. Part. Sci. 43 (1993) 27 [ hep-ph/9302210 ].\n- [28] V.A. Rubakov and M.E. Shaposhnikov, Electroweak baryon number nonconservation in the early universe and in high-energy collisions , Usp. Fiz. Nauk 166 (1996) 493 [ hep-ph/9603208 ].\n- [29] A. Riotto and M. Trodden, Recent progress in baryogenesis , Ann. Rev. Nucl. Part. Sci. 49 (1999) 35 [ hep-ph/9901362 ].\n- [30] D.E. Morrissey and M.J. Ramsey-Musolf, Electroweak baryogenesis , New J. Phys. 14 (2012) 125003 [ 1206.2942 ].\n- [31] A. Kosowsky, M.S. Turner and R. Watkins, Gravitational radiation from colliding vacuum bubbles , Phys. Rev. D45 (1992) 4514.\n- [32] A. Kosowsky, M.S. Turner and R. Watkins, Gravitational waves from first order cosmological phase transitions , Phys. Rev. Lett. 69 (1992) 2026.\n- [33] A. Kosowsky and M.S. Turner, Gravitational radiation from colliding vacuum bubbles: envelope approximation to many bubble collisions , Phys. Rev. D47 (1993) 4372 [ astro-ph/9211004 ].\n- [34] M. Kamionkowski, A. Kosowsky and M.S. Turner, Gravitational radiation from first order phase transitions , Phys. Rev. D49 (1994) 2837 [ astro-ph/9310044 ].\n- [35] C. Caprini, R. Durrer and G. Servant, Gravitational wave generation from bubble collisions in first-order phase transitions: An analytic approach , Phys. Rev. D 77 (2008) 124015 [ 0711.2593 ].\n- [36] S.J. Huber and T. Konstandin, Gravitational Wave Production by Collisions: More Bubbles , JCAP 0809 (2008) 022 [ 0806.1828 ].\n- [37] T. Konstandin, Gravitational radiation from a bulk flow model , JCAP 1803 (2018) 047 [ 1712.06869 ].\n- [38] R. Jinno and M. Takimoto, Gravitational waves from bubble dynamics: Beyond the Envelope , JCAP 1901 (2019) 060 [ 1707.03111 ].  \n- [39] D. Cutting, M. Hindmarsh and D.J. Weir, Gravitational waves from vacuum first-order phase transitions: from the envelope to the lattice , Phys. Rev. D 97 (2018) 123513 [ 1802.05712 ].\n- [40] M. Hindmarsh, S.J. Huber, K. Rummukainen and D.J. Weir, Gravitational waves from the sound of a first order phase transition , Phys. Rev. Lett. 112 (2014) 041301 [ 1304.2433 ].\n- [41] J.M. Quashnock, A. Loeb and D.N. Spergel, Magnetic Field Generation During the Cosmological QCD Phase Transition , Astrophys. J. Lett. 344 (1989) L49.\n- [42] A. Brandenburg, K. Enqvist and P. Olesen, Large scale magnetic fields from hydromagnetic turbulence in the very early universe , Phys. Rev. D 54 (1996) 1291 [ astro-ph/9602031 ].\n- [43] J. Ahonen and K. Enqvist, Electrical conductivity in the early universe , Phys. Lett. B 382 (1996) 40 [ hep-ph/9602357 ].\n- [44] P.B. Arnold, G.D. Moore and L.G. Yaffe, Transport coefficients in high temperature gauge theories. 1. Leading log results , JHEP 11 (2000) 001 [ hep-ph/0010177 ].\n- [45] B. von Harling and G. Servant, QCD-induced Electroweak Phase Transition , JHEP 01 (2018) 159 [ 1711.11554 ].\n- [46] D. Bodeker and G.D. Moore, Electroweak Bubble Wall Speed Limit , JCAP 05 (2017) 025 [ 1703.08215 ].\n- [47] A. Kobakhidze, C. Lagger, A. Manning and J. Yue, Gravitational waves from a supercooled electroweak phase transition and their detection with pulsar timing arrays , Eur. Phys. J. C 77 (2017) 570 [ 1703.06552 ].\n- [48] J. Ignatius, K. Kajantie, H. Kurki-Suonio and M. Laine, The growth of bubbles in cosmological phase transitions , Phys. Rev. D 49 (1994) 3854 [ astro-ph/9309059 ].\n- [49] J.R. Espinosa, T. Konstandin, J.M. No and G. Servant, Energy Budget of Cosmological First-order Phase Transitions , JCAP 06 (2010) 028 [ 1004.4187 ].\n- [50] M. Hindmarsh, S.J. Huber, K. Rummukainen and D.J. Weir, Numerical simulations of acoustically generated gravitational waves at a first order phase transition , Phys. Rev. D92 (2015) 123009 [ 1504.03291 ].\n- [51] M. Hindmarsh, Sound shell model for acoustic gravitational wave production at a first-order phase transition in the early Universe , Phys. Rev. Lett. 120 (2018) 071301 [ 1608.04735 ].\n- [52] M. Hindmarsh, S.J. Huber, K. Rummukainen and D.J. Weir, Shape of the acoustic gravitational wave power spectrum from a first order phase transition , 1704.05871 .\n- [53] P. Niksa, M. Schlederer and G. Sigl, Gravitational Waves produced by Compressible MHD Turbulence from Cosmological Phase Transitions , Class. Quant. Grav. 35 (2018) 144001 [ 1803.02271 ].\n- [54] M. Hindmarsh and M. Hijazi, Gravitational waves from first order cosmological phase transitions in the Sound Shell Model , JCAP 1912 (2019) 062 [ 1909.10040 ].\n- [55] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky, Numerical simulations of gravitational waves from early-universe turbulence , Phys. Rev. D 102 (2020) 083512 [ 1903.08585 ].\n- [56] R. Jinno, T. Konstandin and H. Rubira, A hybrid simulation of gravitational wave production in first-order phase transitions , JCAP 04 (2021) 014 [ 2010.00971 ].\n- [57] J. Dahl, M. Hindmarsh, K. Rummukainen and D.J. Weir, Decay of acoustic turbulence in two dimensions and implications for cosmological gravitational waves , Phys. Rev. D 106 (2022) 063511 [ 2112.12013 ].\n- [58] R. Jinno, T. Konstandin, H. Rubira and I. Stomberg, Higgsless simulations of cosmological phase transitions and gravitational waves , JCAP 02 (2023) 011 [ 2209.04369 ].  \n- [59] A. Roper Pol, S. Procacci and C. Caprini, Characterization of the gravitational wave spectrum from sound waves within the sound shell model , Phys. Rev. D 109 (2024) 063531 [ 2308.12943 ].\n- [60] R. Sharma, J. Dahl, A. Brandenburg and M. Hindmarsh, Shallow relic gravitational wave spectrum with acoustic peak , JCAP 12 (2023) 042 [ 2308.12916 ].\n- [61] J. Dahl, M. Hindmarsh, K. Rummukainen and D. Weir, Primordial acoustic turbulence: three-dimensional simulations and gravitational wave predictions , 2407.05826 .\n- [62] A. Kosowsky, A. Mack and T. Kahniashvili, Gravitational radiation from cosmological turbulence , Phys. Rev. D 66 (2002) 024030 [ astro-ph/0111483 ].\n- [63] G. Gogoberidze, T. Kahniashvili and A. Kosowsky, The Spectrum of Gravitational Radiation from Primordial Turbulence , Phys. Rev. D 76 (2007) 083002 [ 0705.1733 ].\n- [64] C. Caprini, R. Durrer, T. Konstandin and G. Servant, General Properties of the Gravitational Wave Spectrum from Phase Transitions , Phys. Rev. D 79 (2009) 083519 [ 0901.1661 ].\n- [65] C. Caprini, R. Durrer and G. Servant, The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition , JCAP 12 (2009) 024 [ 0909.0622 ].\n- [66] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky and S. Mandal, The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence , Geophys. Astrophys. Fluid Dynamics 114 (2020) 130 [ 1807.05479 ].\n- [67] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal and A. Roper Pol, Circular polarization of gravitational waves from early-Universe helical turbulence , Phys. Rev. Res. 3 (2021) 013193 [ 2011.05556 ].\n- [68] A. Brandenburg, E. Clarke, Y. He and T. Kahniashvili, Can we observe the QCD phase transition-generated gravitational waves through pulsar timing arrays? , Phys. Rev. D 104 (2021) 043513 [ 2102.12428 ].\n- [69] A. Brandenburg, G. Gogoberidze, T. Kahniashvili, S. Mandal, A. Roper Pol and N. Shenoy, The scalar, vector, and tensor modes in gravitational wave turbulence simulations , Class. Quant. Grav. 38 (2021) 145002 [ 2103.01140 ].\n- [70] A. Roper Pol, S. Mandal, A. Brandenburg and T. Kahniashvili, Polarization of gravitational waves from helical MHD turbulent sources , JCAP 04 (2022) 019 [ 2107.05356 ].\n- [71] A. Roper Pol, C. Caprini, A. Neronov and D. Semikoz, Gravitational wave signal from primordial magnetic fields in the Pulsar Timing Array frequency band , Phys. Rev. D 105 (2022) 123502 [ 2201.05630 ].\n- [72] P. Auclair, C. Caprini, D. Cutting, M. Hindmarsh, K. Rummukainen, D.A. Steer et al., Generation of gravitational waves from freely decaying turbulence , 2205.02588 .\n- [73] R. Sharma and A. Brandenburg, Low frequency tail of gravitational wave spectra from hydromagnetic turbulence , Phys. Rev. D 106 (2022) 103536 [ 2206.00055 ].\n- [74] R. Jinno and M. Takimoto, Gravitational waves from bubble collisions: An analytic derivation , Phys. Rev. D95 (2017) 024009 [ 1605.01403 ].\n- [75] A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A.G. Tevzadze and T. Vachaspati, The dynamo effect in decaying helical turbulence , Phys. Rev. Fluids. 4 (2019) 024608 [ 1710.01628 ].\n- [76] A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A.G. Tevzadze and T. Vachaspati, Evolution of hydromagnetic turbulence from the electroweak phase transition , Phys. Rev. D 96 (2017) 123528 [ 1711.03804 ].  \n- [77] D. Cutting, M. Hindmarsh and D.J. Weir, Vorticity, kinetic energy, and suppressed gravitational wave production in strong first order phase transitions , Phys. Rev. Lett. 125 (2020) 021302 [ 1906.00480 ].\n- [78] F. Giese, T. Konstandin and J. van de Vis, Model-independent energy budget of cosmological first-order phase transitions: A sound argument to go beyond the bag model , JCAP 07 (2020) 057 [ 2004.06995 ].\n- [79] F. Giese, T. Konstandin, K. Schmitz and J. van de Vis, Model-independent energy budget for LISA , JCAP 01 (2021) 072 [ 2010.09744 ].\n- [80] Y. Gouttenoire, R. Jinno and F. Sala, Friction pressure on relativistic bubble walls , JHEP 05 (2022) 004 [ 2112.07686 ].\n- [81] A. Azatov, G. Barni, R. Petrossian-Byrne and M. Vanvlasselaer, Quantisation Across Bubble Walls and Friction , 2310.06972 .\n- [82] J. Ellis, M. Lewicki, J.M. No and V. Vaskonen, Gravitational wave energy budget in strongly supercooled phase transitions , JCAP 1906 (2019) 024 [ 1903.09642 ].\n- [83] R. Jinno, T. Konstandin, H. Rubira and J. van de Vis, Effect of density fluctuations on gravitational wave production in first-order phase transitions , JCAP 12 (2021) 019 [ 2108.11947 ].\n- [84] S. Blasi, R. Jinno, T. Konstandin, H. Rubira and I. Stomberg, Gravitational waves from defect-driven phase transitions: domain walls , JCAP 10 (2023) 051 [ 2302.06952 ].\n- [85] D.J. Weir, Revisiting the envelope approximation: gravitational waves from bubble collisions , Phys. Rev. D 93 (2016) 124037 [ 1604.08429 ].\n- [86] K. Subramanian and J.D. Barrow, Magnetohydrodynamics in the early universe and the damping of noninear Alfven waves , Phys. Rev. D 58 (1998) 083502 [ astro-ph/9712083 ].\n- [87] D. Bodeker and G.D. Moore, Can electroweak bubble walls run away? , JCAP 05 (2009) 009 [ 0903.4099 ].\n- [88] T. Krajewski, M. Lewicki and M. Zych, Bubble-wall velocity in local thermal equilibrium: hydrodynamical simulations vs analytical treatment , JHEP 05 (2024) 011 [ 2402.15408 ].\n- [89] H.-K. Guo, K. Sinha, D. Vagie and G. White, Phase Transitions in an Expanding Universe: Stochastic Gravitational Waves in Standard and Non-Standard Histories , JCAP 01 (2021) 001 [ 2007.08537 ].\n- [90] C. Caprini, O. Pujol'as, H. Quelquejay-Leclere, F. Rompineve and D.A. Steer, Primordial gravitational wave backgrounds from phase transitions with next generation ground based detectors , 2406.02359 .\n- [91] A. Roper Pol, A. Neronov, C. Caprini, T. Boyer and D. Semikoz, LISA and γ -ray telescopes as multi-messenger probes of a first-order cosmological phase transition , 2307.10744 .\n- [92] LISA Cosmology Working Group collaboration, Gravitational waves from first-order phase transitions in LISA: reconstruction pipeline and physics interpretation , 2403.03723 .\n- [93] A. Kurganov and E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations , Journal of Computational Physics 160 (2000) 241 .\n- [94] M. Frigo and S.G. Johnson, The Design and Implementation of FFTW3 , IEEE Proc. 93 (2005) 216.\n- [95] C. Caprini and R. Durrer, Gravitational waves from stochastic relativistic sources: Primordial turbulence and magnetic fields , Phys. Rev. D 74 (2006) 063521 [ astro-ph/0603476 ].\n- [96] A. Roper Pol, GitHub project cosmoGW , .  \n- [97] C. Gowling and M. Hindmarsh, Observational prospects for phase transitions at LISA: Fisher matrix analysis , JCAP 10 (2021) 039 [ 2106.05984 ].\n- [98] B.B. Kadomtsev and V.I. Petviashvili, Acoustic Turbulence , Soviet Physics Doklady 18 (1973) 115.\n- [99] C. Caprini et al., Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions , JCAP 04 (2016) 001 [ 1512.06239 ].\n- [100] S. Weinberg, Gravitation and cosmology , John Wiley and Sons, New York (1972).\n- [101] A.S. Monin and A.M. Yaglom, Statistical fluid mechanics: Mechanics of turbulence , vol. 2, MIT press, Cambridge, MA, USA. (1975)."}