Gravitational waves from decaying sources in strong phase transitions
Chiara Caprini Ryusuke Jinno Thomas Konstandin Alberto Roper Pol Henrique Rubira Isak Stomberg
=1





./fig/



=10000




 




  
   





12
16
13
32
34
14
112
23
43















































































 


































































    KOBE-COSMO-24-03, TUM-HEP-1522/24, DESY-24-131a,b]Chiara Caprini,c]Ryusuke Jinno,d]Thomas Konstandin,a,1]Alberto Roper Pol,Corresponding author: mailto:alberto.roperpol@unige.chalberto.roperpol@unige.che]Henrique Rubira,d,2]Isak StombergCorresponding author: mailto:isak.stomberg@desy.deisak.stomberg@desy.de[a]Département de Physique Théorique, Université de Genève,
        CH-1211 Genève, Switzerland[b]Theoretical Physics Department, CERN, CH-1211 Genève, Switzerland[c]Department of Physics, Kobe University, Kobe 657-8501, Japan[d]Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany[e]Physik Department T31, Technische Universität München
        James-Franck-Straße 1, D-85748 Garching, Germanychiara.caprini@cern.chjinno@phys.sci.kobe-u.ac.jpthomas.konstandin@desy.dealberto.roperpol@unige.chhenrique.rubira@tum.deisak.stomberg@desy.de
        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.
        Gravitational waves from decaying sources in strong phase transitions
            [
            September 9, 2024
        =====================================================================
        
        
        
        § 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 <cit.>.
        The forthcoming observing run by the LIGO-Virgo-KAGRA (LVK) collaboration are expected to accumulate more events <cit.>.
        Efforts among Pulsar Timing Array (PTA) collaborations have furthermore unveiled convincing evidence of a stochastic gravitational wave (SGWB) background at nano-Hertz frequencies <cit.>.
        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. <cit.>).
        These breakthroughs in GW detection give us ears to astrophysical events and cosmological history inaccessible through other means of observation.
        
        Looking into the 2030s, the launch of the Laser Interferometer Space Antenna (LISA) mission <cit.>, designed to probe
        GWs in the unexplored milli-Hertz frequency band,
        is setting the stage for a potential
        overhaul of modern cosmology <cit.>.
        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.
        
        Cosmological 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. <cit.> and references therein).
        Of particular interest is the phenomenon of first-order phase transitions (PTs) <cit.> that could have occurred in the early Universe.
        While the Standard Model (SM) predicts a crossover at the electroweak scale <cit.>, many theories beyond the Standard Model accommodate a first-order electroweak PT (see Ref. <cit.> 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 <cit.>.
        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 <cit.>.
        While bubble collisions themselves are an important source of GWs <cit.>, it has been shown in Ref. <cit.> 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,[These results have been found in the absence of magnetic fields
        and for small fluid perturbations. A
        primordial magnetic field generated or present during the
        PT, and/or the production of non-linear fluid perturbations, can efficiently induce vortical motion <cit.> due to the
        high conductivity and Reynolds number of the primordial plasma <cit.>.]
        which is expected to be
        the case unless the PT is dominated by vacuum, e.g., for supercooled PTs <cit.>.
        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 <cit.>.
        
        The production of GWs from fluid perturbations can be
        decomposed in two contributions: sound waves (or acoustic/compressional turbulence) <cit.> and vortical turbulence <cit.>.
        While analytical modeling of GW production
        from fluid perturbations is important <cit.>, 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 <cit.>. 
        The non-linear evolution is inevitable due to  the large Reynolds number of the fluid in the early Universe <cit.>.
        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.
        
        Currently, large-scale simulations have been performed by the Helsinki-Sussex group, which numerically solve a coupled scalar field-fluid system <cit.>.
        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 electroweak 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 <cit.>.
        In this scheme, the microphysics of the wall is introduced as a non-dynamical (although space and time dependent)
        energy-injecting boundary condition within the bag equation of state <cit.>,
        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) <cit.>.
        See also Refs. <cit.> 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.
        
        In this paper, we explore the previously uncharted realm of GWs
        sourced by fluid perturbations induced in
        strong first-order PTs.[We clarify that hereby, by strong transitions, we refer to α = 0.5 (see alpha for definition), still far from the supercooled regime where the scalar field potential energy dominates the energy content of the Universe <cit.>.]
        We also update results for weak and intermediate transitions and compare with other results in the literature <cit.>.
        We perform approximately 1000 simulations involving a parameter scan over wall 
        velocities ∈ [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.
        
        The outline of the paper is as follows.
        In sec:setup, we review the GW production in a PT.
        We first describe the bubble nucleation history and the conservation
        laws of the fluid perturbations in
        nucleation_histhydro_eqs.
        Then, GW_production presents the
        description used to compute the GW spectrum
        at present time from the numerical results,
        more detail is given in sec:GW_prod.
        We review in GW_sw the GW production under the stationary
        assumption, found, for example,
        in the sound shell model (which assumes linear fluid perturbations) <cit.>, and used to interpret
        numerical results <cit.>.
        We then
        present
        a novel model in sw_extended that extends the unequal time correlator (UETC)
        to a locally stationary UETC, allowing us to introduce
        the effect of the source decay.
        In sec:expansion, we extend this model 
        to estimate the effects of the expansion
        of the Universe.
        In Sec. <ref>,
        we discuss updates to the previous Higgsless simulations <cit.>, and summarize
        the physical and numerical setup of the simulations.
        In Sec. <ref>, we discuss the numerical results.
        We first present a convergence analysis for the kinetic
        energy and the integrated GW amplitude in sec:convergence, which is combined with
        a convergence analysis that includes
        the potential effects of underresolving the fluid perturbations,
        presented in sec:kinetic_ed.
        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 decay_K2 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 sec:vort, 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 GW_spec_time, together with estimates
        of the GW efficiency Ω̃_ GW.
        We investigate the GW spectral shape in sec:shape and
        fit the numerical results to
        a double broken power law, providing
        estimates for the positions of the relevant
        spectral scales.
        We summarize our results in sec:summary, and present a template for the GW spectrum.
        
        
        
        
        § 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 nucleation_hist, we review the nucleation history of
        broken-phase bubbles that is used
        in the Higgsless approach, and in hydro_eqs, we review the
        the relativistic hydrodynamic equations.
        In GW_production, we 
        describe how the GW generation is tackled within the Higgsless simulations, and discuss its applicability (see also App <ref> for details).
        In GW_sw, 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 <cit.>. In sw_extended, 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:expansion.
        We validate this novel proposed model in sec:results with the results from
        our numerical simulations, described in sec:numerical_setup.
        
        
        
        
         §.§ Bubble nucleation histories
        
        
        
        In the Higgsless approach <cit.>, a fundamental assumption is that the
        broken-phase bubbles reach a terminal wall velocity  due to the friction exerted by the fluid particles
        <cit.>,
        such that  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
        <cit.>.
        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
        
        β = d/dt(S/T)|_t = t_n = - H  T d/dT(S/T)|_T = T_n  ,
        
        where H is the Hubble scale.
        
        A detailed description of how such bubble nucleation histories are constructed is found in Ref. <cit.>,
        and examples of how modified bubble nucleation histories can be constructed are found in Refs. <cit.>.
        We note that in previous numerical work
        <cit.>, bubbles
        are nucleated simultaneously.
        This is expected to have an impact on the spectral peak and the amplitude, but not on the spectral
        shape <cit.>.
        
        
        
         §.§ 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
        
            ∂_μ T^μν=0 .
        
        
        These 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 conformally invariant and
        hence reduce 
        to those in Minkowski space-time by a conformal transformation <cit.>.
        
        We take T^μν to be that of a perfect fluid 
        
        
            T^μν=u^μ u^ν w-  η^μν
            p  ,
        
        
        where 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 ^2 = 1/3,
        
        
            p=1/3 a T^4-ϵ  ,    w=T d p/d T=4/3 a T^4,
        
        
        with T being the temperature.
        The bag constantϵ<cit.>, defined
        as the difference in vacuum energy density between the symmetric and broken phases, is thus promoted to a time- and space-dependent quantity 
        
        
            ϵ(t, )={[                 0  inside bubbles ,;                 ϵ outside bubbles , ].
         
        
        whose time evolution is uniquely determined for each bubble nucleation history by the terminal wall velocity .
        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 <cit.>,
        and also that strong PTs can lead to runaway behavior <cit.>. 
        
        The relevant quantity for boundary conditions of the fluid at the Higgs interface is 
        the difference in the trace of the energy-momentum tensor, θ = w-4p,
        normalized to the enthalpy <cit.>,
        α = Δθ/3 w = 4ϵ/3 w ,  
        
        where 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 <cit.>.
        The conservation laws for a relativistic perfect fluid are
        
        
            ∂_t T^00+∇_i T^i0    =0  ,
            ∂_t T^j0 +∇_i T^i j(T^μ0, ϵ)     =0   .
        
        
        Note that T^i j(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 <cit.>.
        
        
        
        
         §.§ Gravitational wave production
        
        
        
        The GW spectrum as a present-time observable
        is computed using the following relation
        
            (k) = 1/ρ_totd ρ_GW/d ln k = 
                3   T_ GW   (H_∗/β)^2    I (k) ,
        
        where  ρ_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<cit.>, 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.
        
        The 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 twop_Pi in sec:GW_prod), multiplied by the Green's function of the GW equation.
        We have introduced the prefactor (H_∗/β)^2 in OmGW_Ik to 
        express I in normalized time and wave number units
        t̃≡ t β and k̃≡ k/β
        (see App. <ref>, in particular OmGW_aver, for its derivation from the solution to the GW equation)
        
            I (t̃_∗, t̃_ fin, k̃) = k/2 ∫_t̃_∗^t̃_ fin∫_t̃_∗^t̃_ fin
                E_Π (t̃_1, t̃_2, k̃) cosk̃(t̃_1 - t̃_2)
                t̃_1 t̃_2 .
        
        Here t̃_∗ and t̃_ fin denote the initial and final times of action of the GW source. 
        In Eq. (<ref>) 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 <cit.>.
        We discuss in sec:expansion how to extend our results, obtained with simulations performed
        in Minkowski space-time, to an expanding Universe.
        
        The 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 OmGW and Refs. <cit.>),
        given in the normalized units of the simulation t̃, k̃,
        and Ṽ≡ V β^3: 
        
            I_ sim (t̃_∗, t̃_ fin, k̃) = k̃^3/4 π^2 Ṽ∫_Ω_k̃Ω_k̃/4 πΛ_ijlm
                ()
                [T̃_ij (t̃_∗,
                t̃_ fin, q̃, )  T̃_lm^∗ (t̃_∗,
                t̃_ fin, q̃, ) ]_q̃ = k̃ ,
        
        being Λ_ijlm
        the transverse and traceless operator
        
            Λ_ijlm () = P_il P_jm -  P_ij P_lm ,    with    P_ij = δ_ij - k̂_i k̂_j .
        
        The function T̃_ij (t̃_∗,
            t̃_ fin,q̃, )
        is computed from the normalized 
        stress-energy tensor T̃_ij (t, ), sourcing the GWs, as [see Tij_qk]
        
            T̃_i j(t̃_∗, t̃_ fin, q̃, )
            = ∫_t̃_∗^t̃_ fint̃   e^i q̃t̃∫^3   e^-i · T̃_i j(t̃, )  ,
        
        where T̃_ij (t, ) = w γ^2 v_i v_j/ρ̅, and ρ̅= w̅ + ϵ =   (1 + α)  w̅ is the average total energy density [see wpalpha].
        
        From Weinberg
        we can directly obtain the GW spectrum at present time [see OmGW_Ik]
        using the normalized source T̃_i j(t̃_ init, t̃_ end, q̃, ) 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 have not reached this regime, as in previous numerical work <cit.>, one needs to take into account this limitation
        when interpreting the numerical results, as
        we discuss in sec:results.
        For this purpose, we will use the stationary UETC usually assumed for
        sound waves (see GW_sw)
        and its extension to a locally stationary UETC that can incorporate
        the decay of the source with time (see sw_extended).
        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.
        
        
        
         §.§ 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.[Although
        we refer to sound waves along the paper, which applies to the linearized regime
        of fluid perturbations, we generally consider in our simulations the full
        hydrodynamical system from initial compressional motion that can produce
        large fluid perturbations, especially when the PT
        is not weak.]
        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)<cit.>,
        where K = ρ_ kin/ρ̅ is a time-independent kinetic 
        energy fraction with ρ_ kin = ⟨w γ^2 v^2|.
        Under this assumption, Ik becomes
        
            I (t̃_∗, t̃_ fin, k̃) = 
                k^3   K^2
                ∫_t̃_∗^t̃_ fint̃∫_t̃_∗ - t̃^t̃_ fin - t̃cos (k t_-)   f(t̃_-, k̃) t̃_- .
        
        In the sound-shell model of Refs. <cit.>,
        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
        <cit.>.
        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/3max (, ) being the
        mean separation of the bubbles at the end of the PT <cit.>.
        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. <cit.>, 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
        k R_∗≫√(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. <cit.>), such that
        
            I (t̃_∗, t̃_ fin, k̃) = Ω̃_ GW 
                K^2  β R_∗ τ̃_ sw   S(kR_∗) ,
        
        where Ω̃_ GW corresponds to the GW production
        efficiency and S is a normalized spectral shape,
        such that ∫ln k   S(k) = 1.
        Using OmGW_Ik, the final GW spectrum
        can be written as <cit.>
            Ω_ GW (k) = 3   T_ GW Ω̃_ GW  K^2
                 H_∗ R_∗ H_∗τ_ sw S(k R_∗) .
        OmGW_sshell shows that one could divide the function
        I_ sim^ int, obtained in the simulations [see Weinberg]
        and integrated over ln k,
        by K^2  β R_∗ T̃_ GW to estimate Ω̃_ GW, where T̃_ GW = t̃_ end - t̃_ init is the time interval of the simulation in which
        Weinberg is evaluated, as done 
        in previous numerical
        studies <cit.>.
        Alternatively,
        as we work with non-dimensional
        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
        <cit.>,
        
            Q'(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̃) ,
        
        where the prefactor (ρ̅/w̅)^2 in Qprime_lin
        takes into account that the authors in Refs. <cit.> 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:Fourier transform].
        However,
        we confirm
        with simulations in sec:results
        that using R_∗ in OmGW_sshell to describe the GW amplitude, instead of using Q'/K^2<cit.>, allows us to find a value of Ω̃_ GW that is almost independent of .[References <cit.>
        found a strong dependence of the parameter Q'/K^2 ∼ I_ sim/(K^2 T̃_ GW) with  as the authors did not incorporate the β R_∗ term in the parameterization.
        This led the authors to interpreting that I_ sim∼ K_ξ^2 ξ_ shell T̃_ GW. ]
        
        The kinetic energy fraction K is
        usually taken to be the
        one corresponding to a single bubble, which can be expressed
        as
        
            K_ξ≡κ_ξ α/1 + α ,
        
        where α characterizes the strength of the PT
        [see alpha],
         and κ_ξ≡ρ_ kin/ϵ is the single-bubble
         efficiency factor <cit.>.
        We compare in sec:results the kinetic energy fraction found in the
        simulations with the single-bubble result to take into
        account the effect of collisions and non-linear 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 <cit.>.
        Similarly, Refs. <cit.> have also reported maximum
        values of the kinetic energy fraction in their simulations
        K_ max different than K_ξ.
        
        We find in sec:results that OmGW_sshell
        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 sw_extended, 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 GW_spec_time.
        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.
        
        
        
         §.§ Gravitational waves from decaying sources
        
        
        As we show in sec:results,
        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 OmGW_sshell, impeding its use
        to estimate the GW efficiency.
        In such cases, we propose a  generalization of OmGW_stat, 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_+ =  (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 OmGW_sshell can be substituted by K_ int^2,
        
            K_ int^2 (t̃_∗, t̃_ fin) ≡∫_t̃_∗^t̃_ fin K^2 (t̃) t̃ ,
        
        yielding
        
            I (t̃_∗, t̃_ fin, k̃) = Ω̃_ GW  K^2_ int (t̃_∗, t̃_ fin)    (β R_∗)    S(k R_∗) .
        
        We note that K_ int^2 reduces to K^2   (t̃_ fin - t̃_∗) = K^2  τ̃_ sw when K^2 is constant and we recover OmGW_sshell found
        in the stationary assumption.
        A similar UETC has been recently considered in Ref. <cit.>,
        E_Π(t_1, t_2, k) = 2   k^2  √(P_v(t_1, k)   P_v(t_2, k)) cos (k  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 K2rms under the assumptions discussed
        in GW_sw, such that t_- terms can be
        usually neglected in the integral over t_+
        due to the
        assumed small compact support of t_-<cit.>
        (see discussion in Sec. 5 of Ref. <cit.>).
        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 GW_spec_time, we are validating
        the overall amplitude but not necessarily the GW spectral shape.
        
        After validating this model with the results of numerical simulations
        in GW_spec_time, we can
        estimate the GW efficiency Ω̃_ GW even when the kinetic energy
        is decaying with time.
        In particular,
        we show in decay_K2 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.[Since the GW equation is invariant under
        time translations when expansion of the Universe is ignored,
        we can freely choose the origin of time coordinates t̃_ ref in our simulations.
        Then, choosing t̃_ ref = 0,
        the time intervals become Δt̃ = t̃ - t̃_ ref→t̃. ]
        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 K2rms that the dependence of K_ int^2
        with the source duration, τ̃_ sw = t̃_ fin - t̃_∗, is
        
            K^2_ int =
                 K_0^2  t̃_0  (1 + τ̃_ sw/t̃_0)^1 - 2b - 1/1 - 2b .
        
        This expression reduces to
        K_ int^2 → K_0^2  τ̃_ 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:
        
            lim_τ̃_ sw≫t̃_0 K_ int^2 =     K_0^2  t̃_0/1 - 2 b (τ̃_ sw/t̃_0)^1 - 2 b,         when    b <  , 
            lim_τ̃_ sw≫t̃_0 K_ int^2 =     K_0^2  t̃_0/2 b - 1 ,
                       when    b >  .
        
        Hence, K_ int^2 grows unbounded proportional to τ̃_ sw^1 - 2b 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_0^2  t̃_0/(2b - 1).
        
        
        
         §.§ 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 OmGW_stat2 and Refs. <cit.>).
        Similarly, with the proposed locally stationary UETC,
        the GW amplitude would increase unbounded proportional to
        τ_ sw^1 - 2b when the kinetic energy decays with a decay rate b < 1/2 [see growthrate_K], and it would
        only saturate when b > 1/2.
        
        In the stationary UETC case,
        to take into account the expansion of the Universe,
        the linear increase τ_ sw of OmGW_stat2
        can be substituted by the suppression factor Υ (τ_ sw) = τ_ sw/(1 + H_* τ_ sw)<cit.>,
        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 <cit.>.
        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_∗ a_∗ = H_∗.
        
        In the generalized case when we take into account the decay of K(t) [see OmGW_general],
        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
        OmGW_stat to apply in an expanding Universe <cit.>,
        an effective integrated K^2 that can be used in OmGW_general to
        estimate the effect of expansion is the following
        
            K^2_ int, exp≡ (β/H_∗)^2 ∫_τ̃_∗^τ̃_ finK^2 (τ̃)/τ̃^2τ̃= ∫_0^τ̃_ swK^2 (τ̃_∗ + δτ̃)/(1 + δτ̃/τ̃_∗)^2 (δτ̃)
                 ,
        
        where δτ̃≡τ̃- τ̃_∗ 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 <ref>), Δτ̃= δτ̃+ Δτ̃_∗.
        If we again assume that Δτ̃_∗ = Δτ̃_0,
        the resulting integral for 2 b ≠ 1 can be expressed as
        
            K_ int, exp^2 = K_0^2  Δτ̃_0^2b∫_0^τ̃_ sw(δτ̃+ Δτ̃_0)^-2b/[1 + (H_∗/β)  δτ̃]^2 (δτ̃) = 
                K_0^2  Υ_b (τ_ sw)   (β/H_∗) ,
        
        where we have defined a suppression factor
        Υ_b (τ_ sw) = Δ F_b 
        /(1 - 2b) 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)<cit.>.
        The function
        F_b is the following
        
            F_b (H_∗τ) = (Δτ_∗ + τ/Δτ_0)^1 -2bH_∗Δτ_0/(1 - H_∗Δτ_∗)^2 _2 F_1 [2, 1 - 2b, 2 - 2b, - H_∗ (Δτ_∗ + τ)/1 - H_∗Δτ_∗] ,
        
        where _2 F_1 is the hypergeometric function.
        We highlight that the emergence of a hypergeometric function has no deep physical meaning, since
        hyper arises from introducing the chosen fit K(τ)
        in Eq. (<ref>).
        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 Kexp_fit.
        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 Kexp_fit 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_0^2  Υ_b (H_∗τ_ sw)
          (β/H_∗),
        where the β/H_∗ cancels with the normalization introduced in OmGW_Ik
        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
        
            Ω_ GW (k) = 3   T_ GW Ω̃_ GW  K_0^2
                 H_∗ R_∗ Υ_b(τ_ sw)   S(k R_∗) ,
        
        where Υ_b is obtained from the integrated K_ int^2 and in particular reduces to the expression determined
        by hyper when the fit K^2 = K_0^2   (Δτ̃/Δτ̃_0)^-2b holds at all times of GW production.
        
        As discussed above, in decay_K2 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 OmGW_general applies
        within the duration of our simulations
        in GW_spec_time,
        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.
        
        
        
        § NUMERICAL SETUP
        
        
        In this section, we focus on describing the numerical setup of the Higgsless simulations: in subsec:Updates to the simulation, we comment on the updates in the numerical scheme with respect to Ref. <cit.>, and in sec:params, we describe the simulation suite considered for this work.
        
        
        
        
         §.§ Updates to the numerical setup
        
        
        
        In this section, we highlight three updates to the
        Higgsless simulations with respect to Ref. <cit.>
        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 <cit.> used for the Higgsless simulations, we refer to Refs. <cit.>.
        
        Commencing with  (1), in practice, the
        integral in eq:Fourier transform
        must be computed numerically on the grid of space and time. For the space grid, this is accomplished through a fast Fourier transform routine <cit.>. 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:Fourier transform was approximated as
        
        
            T̃_ij (t̃_ init, t̃, q̃, ) = ∑_t̃^' =
            t̃_init ^t̃δt̃  e^i q̃t̃^'T̃_ij(t̃^', )   ,
        
        
        i.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. <cit.> 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. <cit.>).
        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.
        
        Continuing 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 ith 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
        
        k̃_i = 3 - a/δx̃sin(a π l_i/N)  ,
        
        
        where δ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
        
        
            k̃_i = 
            2 π l_i / (N δx̃ ) ,     l_i < N / 2 , 
             
            0 ,     N / 2 , 
             
            2 π (l_i - N) / (N δx̃) ,     l_i>N / 2 .
        
        
        As 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.
        
        The 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. <cit.>.
        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==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. 
        
        In all other regards, the current version of the Higgsless implementation is identical to the first version in Ref. <cit.>.
        
        
        
        
         §.§ Simulations and parameter choices
        
        
        We list the parameters considered in this study in Tab. <ref>.
        We expand upon Ref. <cit.> 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 - )^-13/10<cit.>, 
        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.
        
        The bubble nucleation histories result in an (asymptotic) number of bubbles
        of the order of 
        N_b ≃L̃^3/(8 π ^3),
        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 nucleation_hist,
        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. <cit.> for details).
        In our simulations,
        L̃/
        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̃/ = 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̃/ = 20 cover a larger dynamical range in the UV regime.
        For comparison, the number of bubbles for L̃/ = 40 (N_b ≃ 2500) in our work
        and previous Higgsless simulations <cit.> is in general larger
        than most[We note that Ref. <cit.> uses N_b = 32558 for
        a weak PT with =0.44, while it takes either N_b = 988, 125, or 37 for the rest of the PT parameter space.
        In Ref. <cit.>, 5376 bubbles are used
        for some weak PTs, while 11 and 84 are used for other
        weak PTs, and for intermediate ones.
        Reference <cit.> considers 8 bubbles
        for all simulations.] of the previous numerical simulations of the fluid-scalar system
        <cit.>, 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  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̃/ = 20.
        
        For each of the 76 parameter points {v_ w, α, L̃/},
        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:convergence).
        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 {, α, L̃/}, we have also run single-bubble simulations to track the convergence of the self-similar fluid profiles, leading to 76 “single-bubble” simulations.
        These results are presented in sec:kinetic_ed and will be used
        to improve the extrapolated predictions of the “reference” multiple-bubble simulations in sec:convergence.
        We note that single-bubble simulations are only run until t̃_ end = L̃/[2max(, )], being roughly the time when the fluid shell reaches the edge of the simulation domain.
        
        In 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 ∈{0.32/0.36,  0.6,  0.8}, where the lower = 0.32 is used for weak and intermediate transitions, and =0.36 for strong ones.
        These velocities correspond to deflagrations, hybrids, and detonations, respectively, except for strong transitions for which also =0.8 corresponds to a hybrid. This implies a total of 3 × 3 × 2 × 4 × 9 = 648seed 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:results.
        
        All 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. <cit.>).
        In this regime, we then compute I_ sim (t̃_ init, t̃_ fin, k̃) that allows us to robustly
        test the scaling of OmGW_general and compute
        the GW efficiency Ω̃_ GW and the spectral shape
        S(k R_∗).
        The time t̃_ init = 16 is shortly after
        the time when the broken phase fills up
        the whole volume of the simulation,
        t̃_0 ≃ 10, for the reference nucleation history
        with L̃/ = 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 decay_K2.
        
        In total, we have performed 1028 simulations, which we summarize in Tab. <ref>, with an estimated time of ∼ 10^6 CPU hours.
        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.
        
        
        
        
        § 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:
        
        
          *  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.[Reference <cit.> also provides estimates of the kinetic energy and the integrated
        GW spectra for α = 0.5, but does not
        present results about the spectral shape.]
        Stronger simulations are more challenging when it comes to numerical 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. <ref>
        a template based on the expected GW spectrum from 
        compressional fluid perturbations
        extended to decaying sources (see GW_swsw_extended) to incorporate
        information from our simulations that can be used for phenomenological studies.
        We show in fig:bubbles an example of a simulation
        for a strong PT that corresponds to
        a deflagration with = 0.36.
        
        
          *  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 OmGW_Ik.
        This motivated the authors in Refs. <cit.> to use the variable Q'
        [see Qprime_lin] to interpret the numerical results.
        In the present work, we instead characterize the numerical results
        based on R_∗ and K_ int^2 [see OmGW_general],
        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 GW_swsw_extended).
        We also provide in sec:expansion a definition of K_ int, exp^2
        that allows to incorporate  a posteriori
        the effect of the expansion of the Universe [see Kexp_fit].
        
        
          *  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
        <cit.> towards a shallower spectrum, for example like the one that is found
        in vortical turbulence, k^-8/3<cit.>.
        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 sec:vort 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:bubbles.
        
        
        To test the validity of our numerical results,
        we will pay special attention to the following points, to be addressed throughout this section.
        In sec:convergence, we study the convergence 
        of our results with respect to the grid spacing, δx̃.
        In decay_K2, we study the time dependence of the fluid kinetic energy fraction K, and fit
        the numerical results to the decaying power-law presented in sw_extended.
        In GW_spec_time, we test the expected
        scaling of the GW spectrum
        with K_ int^2 and R_∗,
        evaluating the evolution of the integrated
        GW amplitude with the source duration,
        and compute the GW efficiency
        Ω̃_ GW, according to OmGW_general.
        We also provide in GW_spec_time an estimate of the
        expected GW amplitude in a flat Minkowski space-time, based on the numerical results of
        sec:convergenceGW_spec_time,
        and  in an expanding background, using the model presented
        in sec:expansion.
        Finally, in sec:shape, 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. <cit.>.
        
        
        
        
         §.§ 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̃/∈{20, 40} have been performed.
        Since our Higgsless simulations use relatively sparse grids compared to simulations with scalar 
        fields <cit.>, 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_ 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  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 =0.6 and the PT is weak, when = 0.72 and the PT has intermediate
        strength, and when = 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 sec:kinetic_ed using single-bubble 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. <cit.>).
        For reference, we show in sec:kinetic_ed (see fig:1d_profiles) the self-similar profiles
        of the fluid perturbations expected to be
        produced by uncollided expanding bubbles <cit.>, computed using  CosmoGW<cit.>.
        
        In the following, we first analyze the convergence of the numerical
        results for multiple-bubble reference runs to provide estimates of the integrated kinetic
        energy K_ int (t̃_ init, t̃_ end) and GW amplitude I_ sim^ int≡∫ I_ simln k, where I_ sim (t̃_ init, t̃_ end, k̃)
        is defined in Weinberg, for initial and final times t̃_ init = 16 and t̃_ end = 32.
        These integrated quantities will be used in GW_spec_timesec:shape 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 sec:kinetic_ed, 
        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:calK_0 definition (see sec:kinetic_ed 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:summary, while we study
        K_0 computed from the reference runs for the remaining of
        this section.
        
        In 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 into 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.[As we show in fig:extrapolation,
        we find this to be the case for all simulations when studying K_0 at
        the reference time t̃_0.
        For the integrated K over time, described by K_ rms,
        this is also the case for all simulations but for
        a few cases with large  when α = 0.5.
        In these exceptional cases, as the self-similar profiles become better resolved, hence increasing
        the numerical K before collisions,
        the stronger decay induced by the increase of
        the non-linearities can eventually overcome
        the initial increase of K,
        leading to a decrease of the integrated K_ rms.
        We study the time evolution of K for different N in decay_K2.
        ]
        This motivates us to use the following function when extrapolating the 
        numerical values of the kinetic energy fraction:
        
        
            K = K_∞/1+(N_*/N)^a ,
         
        where 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 <ref>), 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:E_kin_convergence
        indicates
        the degree of convergence of the numerical results, as the relative error in the kinetic energy ε_K can be expanded as
        
            ε_K ≡K_∞ - K/K_∞ = (δx̃/δx̃_∗)^a +  O (δx̃^2a) ,
        
        where δ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. <ref>).
        
        
        
        The result of applying the convergence analysis based on eq:E_kin_convergence to the reference runs
        with multiple bubbles
        is shown in fig:extrapolation,
        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. <ref> and Fig. <ref>.
        The middle panels of fig:extrapolation correspond to the rms
        kinetic energy fraction, computed from K_ int
        [see K2rms] as T̃_ GW  K_ rms^2 ≡ K_ int^2,
        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̃/ = 20 as these
        cover smaller scales than those with L̃/ = 40
        for a fixed
        N, providing a better resolution of the kinetic 
        energy density in the UV regime.
        We show in Tab. <ref> the resulting
        values of the fit
        parameters of eq:E_kin_convergence for the rms kinetic energy  fraction
        K_∞^ rms, normalized by the single-bubble
        kinetic energy fractions
        K_ξ [see Kxi], and a_K,
        and the relative errors ε_K as defined in rel_errorK (see also footnote <ref>),
        for the set of PTs shown in fig:extrapolation.
        We expect K_ rms^2  T̃_ GW to be the relevant
        kinetic energy fraction entering the computed GW amplitude, I_ sim^ int, according to the model proposed in sw_extended.
        
        We find empirically that the exponent a in (<ref>) 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 corresponds to second order <cit.>.
        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_ CJ), we already find absolute errors below 10%, as indicated in Tab. <ref>.
        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 <cit.> and in PT simulations
        <cit.> 
        is in these
        cases obscured by an exponential decay (see the discussion in sec:shape and the fit used in Ref. <cit.>).
        
        Furthermore, we also apply the extrapolation of eq:E_kin_convergence
        to
        the integrated GW spectrum obtained in the code I_ sim^ int.
        The resulting values as a function of δx̃ are shown in the lower panels
        of fig:extrapolation, while the numerical values of the extrapolated
        I_∞^ int, the fit parameter a_ I,
        and the relative error ε_ I are given in Tab. <ref>.
        We note that for the GW amplitude, the fit of eq:E_kin_convergence 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:extrapolation 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_ sim^ int less sensitive to numerical inaccuracies
        than the estimate of K_ int. 
        
        Finally, we also display in fig:extrapolation 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 OmGW_generalKint_decay.
        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:E_kin_convergence (see upper panel of fig:extrapolation).
        However, this result does not take 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 sec:kinetic_ed.
        The corrected values K_0/K_ξ are presented in fig:kappa_eff_kappa, together with
        a modified “efficiency”κ_0, such that
        
            K_0 ≡κ_0 α/1 + α ,
        
        in analogy to Kxi.
        We find a general trend that K_0/K_ξ≳ 1 when <,
        while K_0/K_ξ≲ 1 when >.
        If we take an average value of K_0 over the PT parameters
         and α, we find
        
            K_0 = 0.84^+0.24_-0.29  K_ξ ,
        
        where 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 = 0.8 (see fig:extrapolation).
        For different PT parameters (α and ), one could take the values of fig:kappa_eff_kappa
        to predict the correction of the resulting GW amplitude.
        According to the locally stationary UETC presented in sw_extended, we find that K_ int^2 ≡ K_ rms^2 T̃_ GW 
        is
        the relevant quantity determining the GW amplitude.
        Therefore, based on the results of Ref. <cit.>,
        Ref. <cit.> used the estimated value K_ rms≃ 0.6   K_ξ.
        However, we note that using the power-law fit presented in sw_extended and validated in decay_K2,
        we can relate the final GW amplitude
        to the corrected values K_0 and the decay rate b
        (as we will do in GW_spec_time).
        
        
        
         §.§ 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:E_kin_evolution 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:convergence (see upper and middle panels of fig:extrapolation, and lower panel of fig:E_kin_evolution). 
        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 sec:kinetic_ed).
        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,[In the Kurganov-Tadmor scheme used in our simulations <cit.>, the numerical viscosity is expected to scale proportional to (δx̃)^3<cit.>.]
        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. 
        
        
        
        Besides this general trend, the opposite is found
        for strong PTs, and intermediate ones with thin hybrid profiles (≲ 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:E_kin_evolution 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.
        
        In 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. <cit.>).
        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,
        
            K (t̃ > t̃_0) = K_0 (t̃/t̃_0)^-b ,
        
        where b indicates the power-law decay rate of K.
        This power-law decay prescription
        accurately fits the numerical data (see fig:E_kin_evolution),
        and we have checked that it remains accurate
        up to t̃_ end = 64
        for an example strong PT with α = 0.5 and = 0.8.
        We define the half-life of the kinetic energy as the
        time when K (t̃_0 + t̃_1/2) =  K_0, i.e.,
        
            t̃_1/2 =   (2^1/b-1)    t̃_0  .
        
        
        We display in Fig. <ref> 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.[This time, we refrain from extrapolating to infinite resolution due to the complex behavior of the index.
        Hence, we use the values found for the largest resolution runs with N = 512, and with the best resolution in the UV regime, L̃/ = 20.]
        In the right panel of Fig. <ref>, 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:decay_2.
        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. <ref>
        the presence of vorticity in our simulations and present some preliminary results.
        
        
        
        For 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 <ref>).
        This 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:1d_profiles),
        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 sec:kinetic_ed), so it is not completely clear whether the obtained
        decay rate b is physical.
        
        The 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 ∈{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 = 0.52, indicating that the decay rate seems to
        be physical (see fig:E_kin_evolution and sec:kinetic_ed).
        Furthermore, a similar decay of the kinetic energy was
        already found for intermediate PT simulations
        of the scalar-fluid system <cit.>.
        
        For 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.
        
        
        
         §.§ Time evolution of the integrated GW spectrum and GW efficiency
        
        
        We show in fig:GWgrowth (upper panel) the time evolution of the
        integrated GW amplitude I_ sim^ int≡∫ I_ simln 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 Weinberg].
        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_ CJ),
        as expected from the usual stationary
        assumption in the sound-wave regime (see GW_sw)
        and as argued in previous numerical work
        <cit.>.
        In these situations, K does not significantly evolve with
        time within the simulations (see fig:E_kin_evolution).
        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 sw_extended.
        To test the validity of OmGW_general, found under this assumption,
        we plot in the lower panels of fig:GWgrowth the
        following ratio
        
            Ω̃_ GW (t̃) =  I_ sim^ int (t̃_ init, t̃)/K^2_ int (t̃_ init, t̃) (β R_∗) ,
        
        with 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_ int^2 → K^2  T̃_ GW.
        
        
        
        For α = 0.0046 and = 0.6, which corresponds to a confined
        hybrid (see fig:1d_profiles), K presents a sharper decay
        with time than
        the other PTs (see fig:E_kin_evolution) 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
        ∈{0.6, 0.68}≲ v_ CJ.
        We see in fig:GWgrowth that in these cases again
        the GW amplitude grows slower than linearly with the
        sourcing time.
        However, the proposed
        ratio Ω̃_ GW (t̃) of OmGW_numerical,
        defined with K_ int^2, is closer
        to be constant with time (see lower panel of fig:GWgrowth), 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:GWgrowth 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 sw_extended.
        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:E_kin_evolution),
        the growth of the
        GW amplitude is
        
            I_ sim^ int (t̃_ init, t̃) = Ω̃_ GW  K_0^2   (β R_∗)
                 t̃_ init  (t̃_0/t̃_ init)^2b [1 + (t̃ - t̃_ init)/t̃_ init]^1 - 2 b - 1/1 - 2 b ,
        
        as found in Kint_decay for the case
        t̃_0 = t̃_ init.
        
        Since we find that Ω̃_ GW (t̃)
        is roughly constant in time after incorporating K_ int^2
        in the scaling of the GW amplitude (see fig:GWgrowth),
        we take this value at the end of the simulations t̃_ end.
        The resulting GW efficiency
        Ω̃_ GW is shown in fig:GW_efficiency
        for different numerical resolutions N and for both box sizes
        L̃/ = 20 and 40.
        We show the values of Ω̃_ GW, computed from OmGW_general
        using the extrapolated values to δx̃→ 0 of I_ sim^ int and K_ int,
        as described in sec:convergence (see also fig:extrapolation).
        We compare the extrapolated efficiencies with those found using
        the sound-shell model <cit.>, under the assumptions described in
        GW_sw (see also App. B of Ref. <cit.>),
        and those obtained from numerical simulations of the full coupled scalar field-fluid system
        <cit.>.
        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. <cit.>).
        We have modified the values of Ω̃_ GW
        from Refs. <cit.> to take into account that they consider β R_∗ = (8π)^1/3, instead of the corrected β R_∗ = (8π)^1/3max(, ) that we use in OmGW_general.
        Furthermore, we note that the integrated Ω̃_ GW might be modified by the structure that
        develops below the peak, described in Ref. <cit.>.
        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 GW_sw,  (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. <cit.>, see results in sec:shape.
        
        Overall, 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:GW_efficiency show the standard deviation obtained from 10 different bubble nucleation histories, corresponding to the “seeds” set of simulations listed in Tab. <ref>.
        
        
        
         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  observed in both the sound-shell model and the coupled scalar field-hydrodynamical simulations.
         This is important for two reasons:  (1) the agreement 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[Reference <cit.> presents results of Ω_ GW/Ω_ GW, exp =  I_ sim^ int/ I_ exp^ int, where I_ exp^ int
         would correspond to the value found using OmGW_sshell with K = K_ξ and Ω̃_ GW = 10^-2.
         The ratio that ref. <cit.> presents therefore corresponds to a combined
         estimate of Ω̃_ GW K_ rms^2/K_ξ^2 and extraction
         of Ω̃_ GW for comparison is not straightforward.] for strong PTs (α=0.5).
         We note that reference data points
         Ω̃_ GW in Refs. <cit.> are computed
         assuming a linear growth with the source duration as in OmGW_sshell.
         Hence, incorporating K_ int as in OmGW_general can modify the
         value of Ω̃_ GW when the source decays.
         The extrapolation method described in Sec. <ref> and presented in Fig. <ref> as solid lines seem to behave very well, delivering agreement between the numerical results from both 
         simulation domains L̃/ = 20 and 40.
        
        For 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 -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. <cit.>, 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. <cit.> might be due to the different nucleation histories
        considered (simultaneous in Ref. <cit.> and exponential in our simulations). Again, extrapolation seems overall good as the extrapolated values agree well for the simulations
        with L̃/ = 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.
        
        We 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:GW_efficiency might still present
        numerical errors related to those
        listed in Tab. <ref>.
        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.
        
        As 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 sw_extended and numerical results in decay_K2).
        Furthermore, the universality
        of Ω̃_ GW is also due
        to the use of OmGW_general instead of Q'/K^2 considered
        in previous work <cit.>, as discussed in footnote <ref>.
        
        The average values of Ω̃_ GW
        over  for each strength α from the simulations are the
        following
        
        10^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 , 
        
        where the super and subscripts refer to the maximum and
        minimum values found in the extrapolated values from our simulations.[Based on the parameterization of OmGW_stat2
        and the results of Ref. <cit.>, Ref. <cit.> reported a value A_ sw = 3  Ω̃_ GW≃ 0.11,
        slightly larger than the extrapolated values in eq:Omegatilde with our updated
        numerical simulations and results.]
        We note that these values only take into account variations of Ω̃_ GW at different values of , but not from
        numerical inaccuracies in our numerical results and hence, in the extrapolated
        values.
        
        Finally, 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̃^-2b in decay_K2, the
        GW amplitude of Isim asymptotically
        grows proportional to 
        I_ sim^ int∼t̃^ 1 - 2b for b <,
        I_ sim^ int∼lnt̃ for b =, and I_ sim^ int∼t̃^ 0 for b >.
        Then, we need to extrapolate the
        resulting GW amplitude by extending K^2_ int to 
        times beyond the final time of the simulation using OmGW_general.
        We note that unless b >, then the GW
        amplitude keeps growing unbounded as long as the UETC
        assumed in GW_sw describes the source dynamics.
        However, we expect that the UETC deviates from this description as vortical motion and turbulence development dominates in the simulation
        <cit.>.
        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.
        
        We present in fig:GWmodel the dependence with the sourcing time τ̃_ sw≡t̃ - t̃_ init of the numerical
        integrated GW amplitude
        I_ sim^ int, compared
        to the one computed analytically for
        the power-law fit K (t̃) = K_0   (t̃/t̃_0)^-b using the
        scaling in OmGW_general to extend the results to times after
        the end of the simulations, and
        using the decay rate b found in fig:decay_2 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[We note by the time t̃_0, significant decay of K has already occurred for strong PTs, as can be seen in fig:E_kin_evolution.
        As GWs can start to be produced from the first time of bubble collisions,
        this will lead to an underestimation of the GW amplitude from our extrapolated results, as discussed in sec:convergence, based
        on a low estimate of K_0/K_ξ.
        However, we avoid extrapolating our results to earlier times than t̃_0 as the decay fit for K^2 does not apply, and because in this regime,
        the results are expected
        to strongly depend on the nucleation history and our
        assumptions on the GW production are not expected to be valid.]
        at t̃_ init = 16, and that the scaling of OmGW_general
        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:GWmodel,
        and we use the extrapolated values K_0 presented in sec:convergence 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:expansion for the values β/H_∗ = 100 and 1000 [see Kexp_general].
        
        
        
         §.§ 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. <cit.> and Higgsless simulations in Ref. <cit.>.
        Utilizing the improved Higgsless code, we update the results of Ref. <cit.> 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.
        
        The findings in Ref. <cit.> indicate that the GW spectrum (k) is characterized by a double broken power law:  
        at small k, a (k)∝ k^3 scaling was observed, which is also expected from causality.
        At large k, the spectrum decays as (k)∝ k^-3. 
        These scalings are in agreement with those found in
        the sound-shell model <cit.>.
        At intermediate scales, a linear scaling regime (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, (k) behaves erratically and is always neglected for the purpose of analysis and parameter extraction. 
        
        To capture the behavior of the spectral
        shape S(k̃) ≡ I_ sim (k̃)/ I_ sim^ int, we use the following double-broken
        power law function,
        
        
            S(k, k_1, k_2, k_)=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,
        
        
        which corresponds to the shape function used
        in Ref. <cit.> 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 ln 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. <ref> equivalent to
        
        
            S_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_)^2,
        
        
        which was previously used in Ref. <cit.>. eq:shape function, 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.
        
        We 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
        
        
            ξ_shell :=ξ_front -ξ_rear  ,
        
        where   
        ξ_front =
        ξ_shock ,     for deflagrations and hybrids , 
         ,
            for detonations , 
        and
        
        ξ_rear = 
         ,
            for deflagrations , 
         ,
            for detonations and hybrids .
        
        The scale R_∗ ξ_ shell is
        expected to determine the peak of the GW spectrum <cit.>.
        
        
        
        
          §.§.§ Fitting to the numerical data
        
        
        We fit Eq. (<ref>) to our numerically computed 
        GW
        spectra and thus extract spectral features from our data.
        We show in fig:fits the numerical GW spectra I_ sim (t̃_ init, t̃_ end, k̃) found
        in the simulations with numerical resolution N = 512 and box sizes L̃/ = 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_ does not represent a physical scale, we do not require that k_2<k_, but allow k_ to take on any value independently.
        In the cases where k_2>k_, the spectral peaks are not resolved properly and suffer from numerical viscosity.
        
        
        
        In obtaining the fit, we neglect the first
        bin for simulations with L̃/ = 20 and the first two
        bins for L̃/ = 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:fits we show fits of eq:shape function to the spectra for different v_w,
        in fig:money plot we show the fitted spectral features 
        k_1, k_2, and k_ 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  and
        determine the scale characterizing the spectral knee and peak.
        We note that the maximum value of the spectral shape used in
        eq:shape function is located at k_ peak, which does not in general
        exactly coincide with k_2 (see discussion in Ref. <cit.>) and their relation depends on the
        fitting parameters.
        We show the resulting spectral peaks obtained from the fit in the
        right columns of fig:money plot.
        
        
        
        Extraction of the scale of exponential damping k_ gives us a handle on the reliability of the measurement of other parameters and the peak; clearly, finding k_2>k_ means we are in 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_ (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. <cit.>.
        Using our numerical results with L̃/ = 20, which present
        better resolution in the UV, averaged over  and 10 nucleation histories, we find the following values for k_2,
        
        k_2    R_∗/2 π≃
        (0.49 ± 0.024) / Δ_w  ,     α = 0.0046 , 
        
        0.93±0.13  ,     α = 0.05 , 
        
        0.45±0.042 ,     α = 0.5 , 
        
        where the indicated uncertainty corresponds to the standard deviation in the measurements among the reference simulations, and
        Δ_w = ξ_ shell/max(, ) 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
        k_1    R_∗/2 π≃ 0.39 ± 0.1  .
        
        We note that both scales k_1 R_∗ and k_2 R_∗ (for intermediate
        and strong transitions)
        present
        very small variability with , indicating a rather universal
        behavior.
        For weak PTs, k_2   R_∗ Δ_w ∼k̃_2 is
        also almost independent of , as expected from the
        sound-shell model.
        The values of k_2 R_∗ and k_1 R_∗ used in Ref. <cit.> are based on the numerical
        results of Ref. <cit.>.
        For weak PTs, we find k_2 R_∗Δ_w consistent with
        the values used in Ref. <cit.>, while
        we find the value of k_1 R_∗ to be
        twice the one used in Ref. <cit.>.
        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:GWgrowth_spec 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 <cit.> and numerical
        simulations <cit.>.
        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. <cit.>, 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:fits and used to provide
        fits of the spectral shape.
        As discussed in GW_sw, the modelling presented and validated
        for the integrated GW spectrum is only expected to hold at wave
        numbers k R_∗≫ (β 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 OmGW_general 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 GWs.
        This can occur within the sound-shell model in the IR regime, as shown in Ref. <cit.>, 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 sw_extended to decaying sources) is no longer valid due to,
        for example,
        the potential development of
        non-linear 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
        <cit.> and we 
        expect that the GW modes would reach their saturation amplitudes
        in this regime.
        
        
        
          §.§.§ 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<cit.>, and our simulations show a clear n_3 ≈ -3
        whenever k_2 ≪ k_.
        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:fits, 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<cit.>, or a shallower acoustic turbulence spectrum <cit.>.
        
        
        
         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_. 
        
        In fig:n_3, 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 . 
        
        
        
        
          §.§.§ 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:old_shape_function,
        and explicitly shows
        large degeneracies among the fitting parameters, thereby thwarting meaningful interpretation and extraction of the relevant spectral features.
        
        Since 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 knee 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̃/ = 40
        for its estimation, whereby a_1=3.6 (i.e., an increase from a_1=2 as used in Ref. <cit.>) is found suitable. Measurements of a_2, on the other hand, benefit from resolving the UV, for which we use exclusively simulations with L̃/ = 20, and find that a_2 = 2.4 (i.e., a reduction from a_2=4 as used in Ref. <cit.>) 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).
        
        The extraction of the parameters k_1, k_2, k_e, and k_ peak are shown in Fig. <ref>.
        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. <ref>.
        The slope of the UV tail is plotted in Fig. <ref>.
        
        
        
        
        § SUMMARY AND CONCLUSIONS 
        
        
        
        We have conducted numerical simulations of cosmological
        first-order phase transitions (PTs) using the Higgsless approach
        <cit.> 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 ∈ (0.32, 0.8).
        These results extend the previous numerical results of Ref. <cit.> to strong PTs, and include a larger number
        of numerical simulations for weak and intermediate PTs.
        
        We 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.
        
        The 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 OmGW_general.
        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 modelling 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.
        
        In 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 sw_extended and validated with our numerical simulations in GW_spec_time, we find the following parameterization
        of the GW spectrum when the Universe expansion can be ignored
        
         (k)  =   3   T_ GW Ω̃_ GW  (H_∗/β)^2   K_ int^2  R_∗β S(k R_∗)   ,
            
        where S(k) denotes the shape function of the spectrum that is normalized to ∫ dln k   S(k) =1, and K_ int^2
        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:base_normalization is a generalization of the 
        parameterization used in the stationary UETC assumption previously
        tested with numerical simulations <cit.> and
        usually assumed for sound-wave sourcing of GWs <cit.> that predicts
        a linear growth with the GW source duration when K does not decay
        with time.
        
        The 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 <cit.>β R_* = (8π)^1/3 max (, )  ,
        
        
        where 1/β parameterizes the duration of the PT,  is the 
        wall velocity, and  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:GW_efficiency and eq:Omegatilde],
         
        10^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 , 
        
        consistent with
        previous numerical simulations <cit.> for weak and intermediate PTs, and with the sound-shell model <cit.>
        for weak PTs.
        For intermediate and strong PTs, we find much
        less dependence with  than for weak PTs,
        clearly showing a departure with respect to
        the predictions of the sound-shell model (see fig:GW_efficiency).
        
        We also provide an estimate of the relevant kinetic energy
        fraction K_0 at the end of the PT using our numerical results
        [see fig:kappa_eff_kappa and K0_fit], given in units of the single-bubble K_ξ [see Kxi], which, averaged over wall velocities, becomes
        
            K_0 = 0.84^+0.24_-0.29  K_ξ .
        
        As a function of , we generally find
        that K_ξ might slightly underestimate K_0 for the smallest
        , while it tends to overestimate it for larger 
        (see fig:kappa_eff_kappa).
        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 <cit.>.
        
        We have studied the decay of the kinetic energy fraction K with
        time t̃ in decay_K2, 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 fig:E_kin_evolutionfig:decay_2).
        For small values of b (hence, negligible decay), one can directly use
        
            K_ int^2 (b = 0) → K_0^2  τ̃_ sw→ K_0^3/2 β R_∗
        
        in eq:base_normalization,
        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).
        
        For 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_ int^2 becomes
        
            K_ int^2 (b, τ_ sw) → K_0^2  β  t_∗ (1 + τ_ sw/t_∗)^1 - 2b - 1/1 - 2b ,
        
        when 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 Kint_const would
        in general underestimate the GW production.
        We compare in fig:GWmodel the numerical
        dependence of the GW amplitude with the source duration τ̃_ sw found in the simulations to the one obtained
        using Kint_decaying, extending
        the analytical fit beyond the time when the simulations end.
        
        As 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 radiation-dominated, 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:expansion).
        Then, as a proxy to estimate the effect of the Universe expansion, we
        can use the following value for K_ int^2 [see Kexp_fit]
        
            K_ int^2 → K_0^2  Υ_b (τ_ sw)   (β/H_∗) ,
        
        which generalizes the suppression factor Υ = H_∗τ_ sw/(1 + H_∗τ_ sw) when the source does not decay <cit.> to any decay rate b using hyper
        for the presented power-law decay fit of K(t̃).
        We also compare in fig:GWmodel the expected evolution of the GW amplitude
        with the source duration according to Kint_exp_decay 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.
        
        Regarding the spectral shape S(kR_∗) in eq:base_normalization, 
        we find that the following template fits accurately our numerical results (see sec:shape)
        
        
        S(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  ,
        
        
        with 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_)^2 [see eq:shape function], effective at k > k_, 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 ≲.
        The slope
        becomes slightly shallower (up to -2.5) for intermediate PTs with supersonic  and strong 
        PTs (see fig:n_3).
        This effect should not play a major role in phenomenological studies but a
        more detailed description is given in Sec. <ref>.
        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 sec:vort.
        
        The 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. <cit.>).
        For weak PTs, the peak follows the 
        thickness of the fluid shells, ξ_ shell, as given in eq:xi_shell,
        while for intermediate and strong PTs the dependence on the wall velocity is 
        much weaker.
        As shown in k2R_fit and fig:money plot, 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(, ) 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 k1_fit and fig:money plot).
        
        
        
        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 https://data.snf.ch/grants/grant/208807182044).
        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.
        
        
        
        
        § 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
        
            s^2 = a^2 [- τ^2 + (δ_ij + h_ij)  x^i  x^j] ,
        
        while the source is active[The final time of GW sourcing  can be taken as the time at which the source stops operating such that all GW modes have
        reached a saturated amplitude in their free-propagation regime.
        However, we note that different modes k can saturate (hence reaching a free-propagation regime) at earlier times than .] (t <) is
        <cit.>
            h_ij (t < , ) = 6 H_*^2/k∫_t_*^tΠ_ij (t', ) sin k (t - t')  t',
        
        where 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_ Pl^2 the critical energy density
        and Λ_ijlm = P_il P_jm -  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.
        
        At later times,
        after the sourcing has ended at , the solution is
        
            h_ij (t ≥, ) = 6 H_*^2/k∫_t_*^Π_ij (t', ) sin k (t - t')  t' .
        
        Then, the time derivatives of the strains h_ij are
        
            ∂_t h_ij (t ≥, ) = 6 H_*^2 ∫_t_*^Π_ij (t', ) cos k (t - t')  t' ,
        
        which can be used to find the
        fractional energy density at present time t_0<cit.>,
        
            (t_0) = ρ_ GW/ρ_ tot^0  = (a_*/a_0)^4/12 H_0^2⟨∂_t h_ij (t_0, )  ∂_t h_ij (t_0, )| .
        
        We consider the GW spectrum (t_0, k) ≡Ω_ GW (t_0)/ln k, which describes
        the two-point correlation
        function of the statistically homogeneous
        and isotropic strain derivatives,
        following the notation of Ref. <cit.>
            (a_*/a_0)^4/12 H_0^2⟨∂_t h_ij (t_0, )  ∂_t h_ij^* (t_0, ')| = (2 π)^6  δ^3 ( - ')   (t_0, k)/4 π k^3 ,
        
        such that (t_0) = ∫(t_0, k) ln k.
        Substituting der_strains into two_point_GW, we find
        
            (t_0, k) = 3 k  T_ GW  H_∗^2 ∫_t_*^∫_t_*^
                E_Π (t_1, t_2, k) cos k(t_0 - t_1)  cos k (t_0 - t_2)
                 t_1  t_2 ,
        
        where 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,
        
            ⟨Π_ij (t_1, )  Π_ij^* (t_2, ')| = (2 π)^6  δ^3 ( - ')  E_Π(t_1, t_2, k)/4 π k^2 .
        
        At present time, for modes k t_0 ≫ 1, we can average the product of Green's functions in OmGW_today_k over oscillations to find
        
            (k) =     3 k/2 T_ GW  H_∗^2 ∫_t_*^∫_t_*^
                E_Π (t_1, t_2, k) cos k(t_1 - t_2) 
                 t_1  t_2 .
        
        Therefore, 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 twop_Pi with the average
        over spherical shells of radius k in Fourier space,
        
            E_Π (t_1, t_2, k) = k^2/2 π^2 V∫_Ω_kΩ_k/4π Π_ij (t_1, )  Π_ij^* (t_2, )  .
        
        Then, substituting EPi_Pi into OmGW_aver
        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
        
            (k) = 3 k^3/4 π^2 V  T_ GW  H_∗^2 ∫_Ω_kΩ_k/4π     ∫_t_∗^Π_ij (t_1, )   e^ikt_1 t_1  
            ×    ∫_t_∗^Π_ij^∗ (t_2, )   e^-ikt_2 t_2 ,
        
        where 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,
        
            T_ij (q, ) = ∫_t_*^ T_ij (t, )   e^iqt t  ,
        
        and using the property Π_ijΠ_ij^∗ = Λ_ijlm T_ij T_lm^∗/ρ̅^2, we find the following expression
        
            (k) =     3 k^3/4 π^2 V ρ̅^2  T_ GW  H_∗^2 ∫_Ω_kΩ_k/4π Λ_ijlm () [T_ij (q, )   T_lm^* (q, )]_q=k ,
        
        previously used in Refs. <cit.> and referred
        to as Weinberg's formula, due to its similarity to the expression obtained
        for the power emitted by isolated deterministic binaries <cit.>.
        However, we note that the applicability of OmGW 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 free-propagation 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.
        
        
        
        
        § 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. <cit.> to develop.
        For reference, we show in fig:1d_profiles the self-similar velocity and enthalpy
        profiles for the PTs considered in the present work (see Tab. <ref>).
        We have previously defined the kinetic energy fraction K such
        that ρ̅  K(t̃) ≡⟨ρ_ kin(, 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 Kxi,
        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,
        
            𝒱(t̃)=V_bp/V ,
        
        we 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 a 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:KoK for all four resolutions N∈{64, 128, 256, 512}.
        
        
        
        In the multiple-bubble runs, we can similarly define ρ̅  K_i (t̃) = ⟨ρ_ kin, i (, t̃)| as the kinetic energy fraction
        of each of the single bubbles i before their corresponding first collision,
        and define the ratio K_i(t̃)/ V_i (t̃),
        where V_i corresponds to the fractional
        broken-phase volume occupied by each
        bubble i.
        
        To 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. <ref>).
        As the convergence of the single-bubble profiles depends on the resolution
        in ξ≡ r/(t - t_i), where t_i is the nucleation time of the bubble i and r the radial distance to the nucleation center,
        we empirically find that doubling the resolution from N to 2N is
        equivalent
        to evaluating the profile at time 2(t̃ - t̃_i) (see also discussion in Ref. <cit.>).
        Hence,
        the kinetic energy of single-bubble simulations (with t_i = 0) obeys K_i^2N(t̃)/𝒱_i(t̃) = K_i^N(2t̃)/𝒱_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̃_end^ sb = L̃v_w/2/max(c_s, v_w).[Note that this t̃_end^ sb (where sb stands for “single-bubble” runs) is always smaller than t̃_ end = 32, the final time of the multiple-bubble simulations. Thus, in producing fig:KoK, we extend the fit of the observed convergence for times greater than t̃_end, enforcing
        that in the limit of infinite time, it converges to the value of K_ξ.
        This extrapolation always represents K_i(t̃) accurately from
        the measured values (below 1% error).
        In any case, since we never use values of K_i(t̃) at times larger than  t̃_coll < t̃_ end^ sb, it does not affect the analysis and it is only
        used to indicate
        the expected convergence of the self-similar profiles in  fig:KoK.]
        
        Then, 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
        
            K_Σ(t̃) ≡∑_i : {t̃_i < t̃} K_i(t̃-t̃_i) ,
        
        which 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.
        
        Similarly, we can construct the fractional broken-phase volume
        occupied by the superposition of single bubbles as 
        
            𝒱_Σ(t̃) ≡∑_i : {t̃_i < t̃}𝒱_i(t̃-t̃_i) ,
        
        which can become larger than one, as it ignores interactions between bubbles.
        However, the ratio K_Σ/ V_Σ is bounded by K_ξ.[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_Σ/𝒱_Σ≳ K_ξ are found before collisions due to numerical inaccuracy
        (see fig:KoK).
        However, as time evolves the fraction asymptotically tends to K_ξ.]
        We plot the time evolution of the ratio K_Σ/ V_Σ
        as dotted lines in Figure <ref> using the nucleation history of the
        reference multiple-bubble simulations with L̃/ = 20.
        
        The 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:KoK 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̃)/𝒱(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 decay_K2, 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.
        
        In fig:KoK, 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
        𝒦_0 = 𝒮  K_0 =  V(t̃_ coll)  K_ξ/K(t̃_ coll)  K_0  .
        
        
        In fig:kappa_eff_kappa, we plot 𝒦_0 obtained
        for numerical resolutions N∈{64, 128, 256, 512}.
        Compared to the extrapolated value
        of K_0 in fig:extrapolation,
        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.
        
        
        
        
        § PRELIMINARY RESULTS FOR VORTICITY
        
        In this section, we present some preliminary measurements of
        the fluid vorticity in our numerical simulations.
        
        
        
         §.§ Vorticity on the lattice
        
        The vorticity is computed as 
        
            ∇×𝐯=(∂  v_z/∂ y-∂ v_y/∂ z) x̂+(∂ v_x/∂ z-∂ v_z/∂ x) ŷ+(∂ v_y/∂ x-∂ v_x/∂ y) ẑ .
        
        On the lattice, we approximate the derivatives
        using first-order central differences,
        
            ∂ v_i/∂x̃_j () ≃v_i[+ δx̃ x̂_j]
            - v_i [- δx̃ x̂_j ]
            /2  δx̃ ,
        
        where δx̃ is the uniform grid spacing.
        
        With this choice of derivative operator, the magnitude of the curl |∇̃×𝐯| is computed at every grid point. 2D simulation slices at different times are shown in  the lower panels of fig:bubbles.
        Note that the definition of the numerical derivative operator of
        eq:central derivative 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:bubbles. 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 largest 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. 
        
        In the lower panel of Figure <ref>, 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-bubble-wall-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, .
        The former convective structures evolve on time scales proportional to the average convection speeds, which are typically much smaller than . 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 <ref>, where additional small-scale structures are hinted in the velocity field but are absent in the enthalpy field.
        
        
        We 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.  
        
        
        
         §.§ Velocity power spectra
        
        
        We generate the Fourier transform v() of the fluid velocity field and then 
        extract the power spectrum from the two-point
        correlation in Fourier space <cit.>⟨v_i ()   v_i^∗ (')| = (2 π)^3  δ^3( - ')   P_v (k) ,
        
        
        where the ensemble average is performed over
        momenta with the same absolute value || = |'| due
        to statistical homogeneity and isotropy. Along the same lines 
        we also construct the power spectrum of the longitudinal modes 
        
        ⟨k̂_i   v_i ()  k̂_j   v_j^∗ (')|
        = (2 π)^3  δ^3( - ')   P_∥ (k) ,
        
        
        where =  saw() / | saw()| is a unit vector according to the  saw description, as discussed in subsec:Updates to the simulation [see eq:momenta mapping].
        Likewise, we extract the
        vortical component of the spectra
        
        
            ⟨[× () ]_i[' ×^∗ (')]_i⟩
                = (2 π)^3 δ^3( - ')   P_⊥ (k) ,
        
        
        such that P_v (k) = P_⊥(k) + P_∥(k).
        
        
        
        
        Figure <ref> 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 strong 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.
        
        
        JHEP
