Shallow relic gravitational wave spectrum with acoustic peak

We study the gravitational wave (GW) spectrum produced by acoustic waves in the early universe, such as would be produced by a first order phase transition, focusing on the low-frequency side of the peak. We confirm with numerical simulations the Sound Shell model prediction of a steep rise with wave number k of k 9 to a peak whose magnitude grows at a rate (H/k p)H, where H is the Hubble rate and k p the peak wave number, set by the peak wave number of the fluid velocity power spectrum. We also show that hitherto neglected terms give a shallower part with amplitude (H/k p)2 in the range H ≲ k ≲ k p, which in the limit of small H/k rises as k. This linear rise has been seen in other modelling and also in direct numerical simulations. The relative amplitude between the linearly rising part and the peak therefore depends on the peak wave number of the velocity spectrum and the lifetime of the source, which in an expanding background is bounded above by the Hubble time H -1. For slow phase transitions, which have the lowest peak wave number and the loudest signals, the acoustic GW peak appears as a localized enhancement of the spectrum, with a rise to the peak less steep than k 9. The shape of the peak, absent in vortical turbulence, may help to lift degeneracies in phase transition parameter estimation at future GW observatories.


Introduction
Relic gravitational waves (GWs) from the early universe can reveal valuable information about the underlying generation mechanism and thus the physical conditions at that time [1][2][3].A particularly interesting epoch is the electroweak (EW) era, which may have involved a firstorder phase transition [4][5][6][7].A first-order phase transition is characterized by the formation, expansion, and subsequent merging of bubbles containing the low-temperature phase, leading to a transition of the entire universe to a new phase, and the release of latent heat.During this process, the kinetic energy transferred to the plasma can be a considerable fraction of the total available energy and would therefore be an important source of GWs, peaked at a frequency set by the inverse of the mean bubble spacing.For GWs from the electroweak phase transition, the relevant frequencies lie in the mHz range, which is accessible to the Laser Interferometer Space Antenna (LISA) [8].Studies of GWs from that epoch have therefore attracted considerable attention.
The generation of GWs during a phase transition is often divided into three stages: (i) the bubble collision phase, (ii) the acoustic phase, and (iii) the turbulence phase.The contribution from the bubble collision phase is typically small compared to that from the other two stages, except in the case of a vacuum phase transition, where bubble collision becomes the only source of GW production [9][10][11][12][13].There has been extensive research on GW production by fluid flows in the early Universe through both numerical simulations [14][15][16][17][18][19][20][21][22] and semi-analytical or analytical models [23][24][25][26][27][28][29][30][31][32][33].For phase transitions where the kinetic energy fraction is small, or for those proceeding by detonations, simulations indicate that the kinetic energy in vortical modes is subdominant compared to the acoustic or compressional modes [19].This paper concerns the acoustic contribution to the GW spectrum, and focuses specifically on its shape to the low-frequency side of the peak about which there is some uncertainty in the literature.
Relic GWs can be characterized by their normalized energy spectrum, Ω GW (k), where k is the wave number, and Ω GW d ln k is the fraction of the critical energy density of the universe in gravitational waves.The GW spectrum depends on the spectrum of the hydrodynamic stress, which depends on the spectrum of the velocity field, how the stress is correlated in time, and on how quickly the stress appears [25].
Velocity spectra can be characterized by the wave number of the peak k p (also referred to as the energy-carrying wave number), a declining inertial range spectrum for k > k p , and a rising subinertial range spectrum for k < k p .When the subinertial range spectrum is blue (steeper than that of white noise), the spectrum of the stress is always white -regardless of how blue the turbulence spectrum is [34].
Subject to certain assumptions about the time-dependence of the stress amplitude and its correlations, various arguments can be applied to conclude that the GW spectrum is peaked at wave numbers around k p , and that the behavior of the GW spectrum on the low-frequency side of the peak should be linear [18,20,25,30,35].The linear behavior applies down to a frequency set by the inverse source duration or the Hubble rate, whichever is the shorter.Below this frequency the GW spectrum is white noise.
The most developed semi-analytical model for acoustic production of GWs is the Sound Shell model [26,36].Studying fast transitions, where the peak wave number is much less than the inverse Hubble length, it was found that sound (or acoustic) waves produce a very steep k 9 GW spectrum at low frequencies.Other models, based on modelling the flow by expanding shells of shear stress, have indicated a k 1 spectrum at wave numbers below the peak.Direct numerical simulations of the sound wave phase of the phase transition conducted in ref. [17,37] show a peak, but no steep rise.It remains unclear whether this discrepancy arises due to the inability of the simulations to capture the infrared part of the spectrum accurately, due to limited volume and duration, or for some other reason.
The purpose of the present work is therefore to reconsider the model of ref. [36] with a range of different approaches to clarify the origin of the apparent conflict.Recently, the authors of ref. [38] studied similar questions and find a shallow part in the GW spectrum.While their findings are consistent with ours, we provide a more detailed investigation of how this shallow part appears and also complement our findings by presenting direct numerical simulations of irrotational flows.Furthermore, we compare the GW spectra generated by acoustic and vortical flows in expanding backgrounds with similar stress spectra, confirming that the peak is still visible in the acoustic case, and distinguishes the two types of flow.It is worth noting that the acoustic spectra we study do not capture the onset of turbulence, which may lead to new power laws emerging near the peak [21].
In section 2, we discuss the different approaches used to calculate the GW spectra originating from sound waves.In section 3, we present our findings and compare the results obtained using different approaches discussed in section 2. We discuss our findings in section 4. In Appendix A and Appendix B, we give details on the origin of the linearly rising GW spectrum, while in Appendix C we report on a check on the effect of the growth rate of the stress.In Appendix D, we discuss an initial condition with non-uniform energy density.
Throughout this paper, we adopt units, where the speed of light is unity.
2 Approaches to computing the gravitational wave spectrum

Gravitational wave spectrum
Gravitational waves are represented by the transverse and traceless components of metric perturbations.In our analysis, we consider the background to be a homogeneous, isotropic, and spatially flat expanding universe.The metric describing such a universe with tensor perturbations can be expressed as where a(η) represents the scale factor and η denotes the conformal time.The evolution of h ij is determined by the Einstein equation.By employing this equation, we obtain the following linearized equation of motion for the Fourier space representation of h ij , Here, tildes symbolize quantities in Fourier space; this convention is consistently used throughout this paper.In the above equation, a 2 ΠT T ij represents the transverse traceless part of the energy-momentum tensor Πij of a GW source.In terms of the stress tensor normalized by the energy density at the initial epoch (ρ * ) and Hij = a hij , the above equation reduces to, where Tij = a 4 Πij /(a 4 * ρ * ) is the normalized stress and a * represents the value of the scale factor at the initial epoch.In the radiation-dominated era, using a = a * (η/η * ), and after replacing k → k/(a * H * ) and η → η/η * (η * and H * denote the conformal time and Hubble parameter at the initial epoch, respectively), the above equation reduces to Here, G(η) = 6a * /η.In the limit of wave number large compared with the Hubble rate (k → ∞), one can make the approximation G(η) = 6a * /η * , equivalent to a non-expanding background.With our choice of units and scale factor, G(η) = 6 in a non-expanding background (see ref. [39] for details).Usually, one is interested in estimating the GW spectrum at the present day, Ω GW (k).It represents the GW energy density per logarithmic wave number interval normalized by the present-day critical energy density of the universe.For the case when the source is active during the interval η i to η, where η denotes the time until the turbulence is active, Ω GW (k) is given by [40] where U T (k, η 1 , η 2 ) is the normalized fluid shear stress unequal time correlator.Here H 0 and a 0 denote the Hubble parameter and the scale factor at the present epoch, respectively and (2.6)

The non-expanding universe case
When the effective duration of the GW source is shorter than the Hubble expansion time, it is possible to neglect the expansion of the universe in estimating the GW energy spectrum.In such cases, Ω GW is given by (2.7) The above expression is similar to the expression obtained in ref. [36] after substituting their equation (3.11) into their equation (3.6).The only difference is that in our case we have normalized the stress tensor with the energy density of the initial epoch and the GW energy density with the present-day critical energy density.
The main reason for considering the non-expanding case is that we want to assess the validity of the approximations used in ref. [36], where a non-expanding universe was assumed.This allows us to provide a more detailed understanding of the origin of particular features in the GW spectra.

Calculation based on the sound shell model
The evaluation of the GW energy spectrum requires determining U T (k, η 1 , η 2 ) resulting from colliding sound waves generated by a phase transition in the early universe.The stress tensor for the fluid is given by where w = ρ + p is the enthalpy density consisting of the energy density ρ and the pressure of the fluid p.Here, u i and γ = (1 − u 2 ) −1/2 are the components of the fluid 3-velocity and Lorentz factor, respectively.The sound shell model of ref. [36] assumes that the expanding shells of pressurized plasma surrounding bubbles of the new phase continue to propagate after the phase transition is completed, each of them acting as an initial condition for a sound wave.The velocity field generated by a phase transition involving a large number of expanding bubbles is then at each point a collection of sound waves resulting from many of these sound shells.The velocity can then be treated as a Gaussian random field.Therefore, in calculating U T (k, η 1 , η 2 ) for a non-relativistic fluid using the standard method of expanding the resulting connected four-point correlator using the Wick expansion, any non-Gaussianities are thus ignored.The four-point correlator then reduces to a sum containing products of two-point velocity correlators.The velocity correlator is assumed to be curl-free, like the flow around a single bubble, and to satisfy the linearized wave equation.In Fourier space, it then takes the form ũi (q, η 1 )ũ where (adopting conventions of non-relativistic fluid dynamics) E K (q) is the kinetic spectrum per linear wave number interval, ω = c s q is the angular frequency, and c s = 1/ √ 3 is the sound speed in the plasma, assumed ultrarelativistic.The kinetic spectrum obeys , where V represents a volume average.We write the kinetic spectrum per logarithmic wave number explicitly as qE K (q).
The result for U T (k, η 1 , η 2 ) is equation (3.34) in ref. [36], which can be written as where w = ρ + p is the mean enthalpy density and µ = (q 2 + k 2 − q2 )/2kq is the cosine of the angle between k and q.In our calculations, we assume a kinetic spectrum of the form with k p being the nominal peak wave number, u rms the root-mean-square velocity, and ẼK (k) the normalized kinetic spectrum.Whenever we plot spectra from a simulation, we compute them for wave numbers sampled uniformly in k.For the above form of the kinetic spectrum, the actual peak of the spectrum is at 2 1/6 k p ≈ 1.12 k p .This form is consistent with irrotational and causal [41] flows with shocks [21], as appropriate for phase transitions.The peak wave number is inversely proportional to the mean bubble spacing.
We assume that the kinetic spectrum appears instantaneously, and remains at a constant amplitude.In Appendix C, we consider a simple model for the growth of kinetic energy during the phase transition.By substituting U T (k, η 1 , η 2 ) from equation (2.10) into equation (2.7), the GW spectrum can be obtained as where and (2.14) Here, Ω r represents the radiation energy density fraction at the present epoch, g * s and g 0s denote the effective degrees of freedom in entropy at the initial and the present epoch and g * and g 0 are the corresponding effective degree of freedom in energy density.In equation (2.12), the time dependence is embedded into a kernel ∆, which now has the form where ω ±± = k ± c s q ± c s q, and the two different previously discussed cases regarding the expansion of the universe have been taken into account by defining (2.17) The evaluation of the two time integrals can be carried out analytically in both cases.For the non-expanding (i.e.Minkowski) background, the result can be written in the form while for the radiation-dominated expanding universe the result is (2.19) where Si and Ci are the sine and cosine integrals.Hence, to obtain the shape of the spectrum, only the q and q integrals need to be evaluated numerically.All other numerical prefactors have been absorbed into the constant Ω 0 .
In the analytical study of the GW spectrum in the non-expanding case in ref. [36], the term with ω −− in Eq. (2.18) was shown to grow linearly with conformal time after the q integral is performed, and was argued to dominate.This is consistent with the growth observed in Minkowski space simulations [15][16][17]22].The linearly growing term was shown to behave as k 2n−1 for a logarithmic kinetic spectrum kE K (k) ∼ k n .With n = 5 for k ≪ k p and n = −1 for k ≫ k p , the linearly growing contribution to the GW spectrum has a steep k 9 spectrum at low wave number and a k −3 spectrum at high wave number.In between there is a peak at a similar wave number to that of the peak in the kinetic spectrum, k p .In view of its origin, we call it an acoustic peak.
Here, we consider all terms in equations (2.18) and (2.19).The integrations were carried out using the Nintegrate routine in Mathematica with the adaptive Monte Carlo method, which has proven efficient for the case at hand.The initial time has been chosen so that η i = 1 and the upper limit in the q integral has been taken to be a large finite value of approximately 10 k p .Spectra produced by the adaptive Monte Carlo method have been compared with those produced by quadrature integration schemes, and the results are found to be in a good agreement with each other with respect to accuracy.Our numerical analysis demonstrates that, below the peak, there is indeed a slope of k 9 , although it does not extend to very low wave numbers.Instead, there is a transition to a linear spectrum at a point that depends on the product k p η in the non-expanding case, and on k p in the expanding case. 1  We will provide a detailed discussion of these findings in section 3.1.

Simulations with a kinematic velocity field
To compute the GW field numerically, we need to evaluate the hydrodynamic stress in regular time intervals, δη.In this section, we construct a three-dimensional, time-dependent, random irrotational velocity field in Fourier space, ũ(k, η), as where φ(k) is the Fourier transformed scalar potential of the velocity, ω(k) = c s k is the dispersion relation for sound waves, k = |k| is the wave number.The formulation in equation (2.20) ensures that at η = 0, the pressure is uniform.We construct φ(k) such that kE K (k) has a k 5 subinertial range for k < k p and a k −1 inertial range for k > k p .The potential φ(k) includes a phase factor exp iϕ(k) with random phases in the range −π < φ(k) < π.The amplitude of φ(k) is chosen such that the resulting kinetic spectrum has a broken power law (for details, see section III B of ref. [42]).
The horizon wave number at the time of generation is k = 1.In our calculations, we sometimes include wave numbers below the horizon wave number in order to see the subsequent evolution from the Ω GW ∝ k 3 spectrum at early times to a linearly increasing one at later times for kη > 1.
At each time step, we Fourier transform ũ(k, η) into real space, u(x, η), and compute Π ij (x, η).From this, we compute the transverse tracefree (TT) projection T TT ij (k, t) in Fourier space and express it in terms of plus and cross polarization modes, denoted with the subscript A ∈ {+, ×}.We then compute the strain field polarizations HA (k, η) by solving [18] where G = 6 in all our calculations without expansion, and G = 6/η in our calculations with expansion.For the simulations, we use the Pencil Code [43], which is a massively parallel public domain code, where the relevant equations have already been implemented [39].The computational domain is a cube of size L 3 , where k 1 = 2π/L is the lowest wave number.
Although ũ(k, η) can be computed with high precision, regardless of the choice of δη, we cannot choose the timestep too large, because otherwise the accuracy of H(k, η) will become poor.The error in the solution for H(k, η) is O(δη 2 ), and the coefficient in the error is smaller than the third-order accurate solution of the hydrodynamic equations, when those are computed; see section 2.5 below.

Direct numerical simulations with the Navier Stokes equation
The purpose of the kinematic velocity fields discussed in section 2.4 is to bridge the gap between actual turbulence, which is nonlinear, and the linear random sound waves considered in ref. [36].To model turbulence more realistically, we also solve the compressible Navier-Stokes equations directly, i.e., where are the components of the strain rate tensor, and ν is the viscosity.We construct our initial condition in Fourier space by multiplying a random vector field with a superposition of a vortical and an irrotational contribution, where 0 ≤ p ≤ 1 quantifies the irrotational fraction.This tensor can also be written as In this work, we consider the extreme cases p = 0 and p = 1 for vortical and irrotational flows, respectively.In all cases, the initial density is usually chosen to be uniform and equal to ρ.In the irrotational case, however, this leads to a collection of standing instead of traveling waves.We therefore also consider non-uniform initial energy density distributions in accordance with the continuity equation where the initial Fourier-transformed logarithmic density is given by ln ρ(k) = k • u/(c 2 s k).We provide a detailed illustration of this in Appendix D, where we compare with the case of a uniform initial energy density.The relevant parameters for our simulations are the root-mean-square (rms) velocity, the Reynolds number, Re = u rms /νk p , and the wave number k p of the spectral peak of the initial velocity spectrum.We recall that k = 1 corresponds to the horizon wave number at the initial time, η i = 1.

Spectral evolution with time
In figure 1, we compare Ω GW (k, η)/kΩ 0 for η = 2, 10, 50, and 250 in the non-expanding universe case.The division by wave number makes a linear growth with k clear.The spectra were obtained by numerically integrating equation (2.12).In the calculations leading to figure 1, we used equation (2.11) for the kinetic spectrum with k p = 30.We see that Ω GW (k, η) consists of a marked peak on top of a flat part with Ω GW (k)/k ∝ k 0 .As demonstrated below using numerical simulation, the peak does not appear for vortical flows.The flat part corresponds to Ω GW (k) ∝ k and is therefore also referred to as the shallow part of the GW spectrum.The acoustic peak increases approximately linearly with time.It does indeed have a k 9 slope below the peak over a short range 10 < k < 30.Above the peak, Ω GW (k, η) falls off approximately like k −3 .This corresponds to Ω GW (k, η)/k ∝ k −4 , as is indicated in figure 1.This normalization shows more clearly the extent of the flat part with The reason why we see the k 9 subrange only for a relatively short k range is the presence of an approximately stationary and linearly rising part with a shallow slope Ω GW (k) ∝ k for 1 < k < 10.This shallow rise was not present in the earlier work of ref. [36].As demonstrated in Appendix A, the shallow part emerges because of additional contributions to the time integral that are not negligible in practice.Indeed, as mentioned in the introduction and discussed in Appendix A, a linearly increasing GW spectrum is expected on quite general grounds.Furthermore, as time goes on, this linear GW spectrum continues to extend to progressively lower k values as the horizon grows, smaller k values become causally connected, and GWs begin to oscillate [31].
From figure 1, it appears that the linear scaling part of the GW spectrum remains unchanged with time.This observation raises the question whether this constancy is a result of the initial conditions we used.However, it is important to note that the linear scaling part develops gradually with the saturation occurring over a duration equal to 2π/k for the wave number k.We have provided a time evolution of the Ω GW /kΩ 0 for wave numbers within the region of linear scaling part in Appendix B. In this analysis, we consider a constant kinetic spectrum, a choice based on the assumption that the kinetic energy turns on within a time duration shorter than 1/k p .With time, the acoustic peak with its Ω GW (k) ∝ k 9 rise continues to grow, while the shallow part of the spectrum remains at constant amplitude.The wave number of intersection between the two power laws therefore moves toward lower wave numbers as η 1/8 .The horizon wave number moves toward lower k like 1/η.Therefore, the shallow part with Ω GW (k) ∝ k broadens with time.However, since the height of the Ω GW (k) ∝ k 9 feature continues to grow, this steep part of the spectrum and the acoustic peak does indeed constitute an important contribution to the spectrum.To determine the ratio between the peak value of the spectrum and the shallow part, it is necessary to estimate the duration of the interval when the source remains active.In this article, we primarily focus on the contribution from the sound modedominated stage of the phase transition.The sound mode contribution remains active until the development of turbulence, and the relevant timescale for this is the eddy turnover time corresponding to the peak of the spectrum.This timescale is equal to the inverse product of the peak wave number of the kinetic spectrum and the root-mean-square (rms) velocity of the fluid, i.e., the nonlinear timescale η NL = 1/(k p u rms ).If this timescale exceeds the Hubble expansion time, the effect of expansion on GW production becomes important.
The Hubble expansion time sets an upper bound on the effective duration of the sound mode contribution.We discuss this further in section 3.4.This scenario may be applicable to weak phase transitions and may be of observational relevance.For instance, when the peak of the kinetic spectrum occurs at a wave number 10 times greater than the Hubble rate and the rms fluid velocity is 0.1, the value of η NL relevant for turbulence development, becomes comparable to the Hubble expansion time.In section 3.4, we provide a comparison considering the expansion for this specific case.
In figure 2, we show the compensated GW spectra k p Ω GW and k 2 p Ω GW .This allows us to see that the height of the peak of the GW spectrum decreases like 1/k p and the shallow part of the GW spectrum decreases like 1/k 2 p .Thus, the height of the peak relative to the shallow part of the spectrum grows linearly with k p .η =1.

Non-expanding case with synthesized random sound waves
We now compare with the results from our model where the stress is constructed from synthetic sound waves; see figure 3. We should emphasize here that the GW energies from the semianalytical method and the three-dimensional simulations agree with each other rather well.We see that the kinetic spectra have the expected E K ∝ k 4 subrange for k < k p , followed by a k −2 subrange for k > k p .(This k −2 subrange is expected for acoustic turbulence, especially in the presence of shocks [44], although in the present case, shocks are only present when we solve the Navier-Stokes equation.) In figure 3, the GW spectra per linear wave number interval, Ω GW /k, are shown at two different times η = 1.0 and η = 100.At the early time η = 1, the spectrum starts having a shallow subrange and a steep part and these features become more pronounced at late times; see figure 3(d).It is evident from this figure that the steep and large wave number part of Ω GW is very well described by the spectra obtained by our semianalytical model described in section 2.3.

Results from direct numerical simulations
In figure 4, we show the kinetic spectrum E K , the stress spectrum Sp(T ) ≡ U T (k, η, η) k 2 /2π 2 , and the compensated GW spectrum Ω GW (k)/kΩ 0 , at different times for flows starting both from vortically and irrotationally initialized velocity fields.The parameters of these runs are summarized in table 1, where E GW = Ω GW (k, η * )/k dk is the total GW energy in units of the critical density of the universe at the initial epoch.The value of the Reynolds number ≈ 40 is small compared with what is expected in reality.This is not a major concern as we are not interested in the development of turbulence, but in the relation between the GW spectrum and the kinetic spectrum.
Looking at figure 4, we see that Ω GW (k)/kΩ 0 has a bump at k ≈ 10...20 for irrotational turbulence, while for vortical turbulence there is no such bump.Also, Sp(T ) never shows a bump.The bump in Ω GW (k) is potentially an important characteristic of acoustic turbulence in the GW spectra.We also see that the inclusion of low wave numbers in the simulations (k < 1 for Runs A1 and V1) results in a clearer representation of the Ω GW (k) ∝ k 3 range, which is not visible for the other runs where k 1 = 1 or larger.
The existence of a shallow part and the transition to Ω GW (k) ∝ k 3 for very small k are independent of whether the turbulence is acoustic or vortical.In all these cases, the height of the bump is about one order of magnitude.The relative height is also independent of the overall amplitude of the turbulence.
At early times, some of the spectra show a wavy structure in E K (k).These are not caused by numerical artifacts, but are a consequence of having initialized a velocity pattern with zero energy density fluctuation, so all modes have the same phase.This causes a modulation of the spectrum of the form cos kξ(t), where ξ(t) = c s t is the distance a sound wave has propagated in the time t since the initial condition was applied.This type of wavy feature or ringing in the spectrum was explored and explained in a different context in more detail in ref. [45].As shown in Appendix D, however, this ringing phenomenon is strongly reduced when the density is initialized such that the continuity equation is obeyed.
Table 1.Summary of the parameters of the direct numerical simulations discussed in the paper.All runs have a numerical resolution of 1024 3 mesh points.

of expansion in the sound shell model
As discussed in section 2.1, when the expansion of the universe is taken into account, the effective duration of the acoustic source is limited and the steep part of the GW spectrum is much smaller.This is shown in figure 5, where we show Ω GW (k)/k for k p = 10 at η = 5, 15, and 30.We see that the spectra do not change much at late times after η = 4, and that the results for η = 8 and 15 are almost indistinguishable and close to those for η = 4. Furthermore, the height of the bump is limited to about one order of magnitude.The presence of a linearly rising Ω GW (k) is not immediately apparent in figure 5.However, it becomes more pronounced in figure 6 as we increase the value of k p to 30.We also estimate the spectral index of Ω GW (k) just below its peak.For this purpose, we fit a single power law within the k values approximately in the range 0.5 k p to 0.9 k p .The index is  determined through least-squares fitting.The results for the time η = 10 are summarized in table 2. At time η = 10, the spectrum is almost in a saturated state.From table 2, we conclude that, as we increase the value of k p , the GW spectrum below its peak tends toward a k 9 scaling.In the right panel of figure 6, we show the time evolution of k p Ω(k p ) for different values of k p .It is evident from the figure that all of these cases have a similar time evolution and show saturation after η ≈ 5, independently of the value of k p .

Conclusions
The present work has confirmed that there is indeed the very steep Ω GW ∝ k 9 rise to the peak in the GW spectrum for acoustic gravitational wave production, as originally proposed in ref. [36].However, it may not be very prominent in practice, because it is superimposed on a stationary linearly rising part for all subhorizon wave numbers above k > 1/η and below the spectral peak at k ∼ k p .The absence of this shallow part in earlier analytical calculations is a consequence of having considered only the leading term in Eq. (2.15), as discussed in more detail in Appendix A. The height of the acoustic peak grows linearly in time, so it becomes distinct only at late times, and the k 9 subrange exists only in the range k p η −1/8 k k p .The k 9 subrange therefore becomes prominent only when the flow is long-lived and k p ≫ 1.
In Appendix C we considered the effect of the growth rate of the hydrodynamic stress on the shape of the GW spectrum, using a simple model where the stress approaches its final value exponentially with time.This allows analytic expressions for the time integrals to be maintained.We find that the amplitude of the shallow part of the spectrum is reduced, and decreases the slope below linear.The steep rise to the peak is unaffected.
Given that the k 9 subrange takes time to emerge, we can ask whether indications of it have still been seen in earlier work on acoustically generated GWs.Acoustic GW production has been considered in several direct numerical simulations [18,22,36,46], but none has reported a k 9 slope.However, a steepening of the slope leading to a peak with time has been observed in ref. [22], up to approximately k 5 .Our analysis indicates that this steepening would continue in a longer simulation with more dynamic range between the peak wave number of the velocity field and the inverse grid size.On the other hand, vortical flows show no such peak.[18,20,46].Both the presence and the shape of the peak therefore probe the nature of the flow (irrotational or vortical) and provide new information about the value of k p , which may help lift potential degeneracies in parameter estimation identified in ref. [47].More comparative studies are needed before the shape of the peak in the GW spectrum can be used as a diagnostic tool in future observational studies with LISA.Considering the limit of large k p /k, the integral is dominated by the upper limit on x, while the y can be replaced by x, and x − y is O(1).Hence, given that ρ(1, x, x) ∼ x −2 , we have This linear scaling with k/k p is shown in figure 2.

B Early time dependence of the shallow part of the GW spectrum
In this appendix, we present the time evolution of Ω GW /(kΩ 0 ) for wave numbers within the linear scaling region of the GW spectrum.Specifically, in figure 8, we show the timedependent behavior of Ω GW /(kΩ 0 ) for values of k equal to 2, 4, 6, and 8, for the case in which k p = 30.In the right-hand panel of this figure, we normalize the abscissa using a wave number-dependent time scale denoted as η k = 2π/k.This normalization suggests an inverse relationship between the saturation time scale and wave number.

C Effect of kinetic energy growth period on the GW spectrum
In section 3.1, we assumed that the kinetic energy remains unchanged throughout the evolution.We have found that such a case produces a GW spectrum characterized by a steep decline just below its peak, smoothly transitioning to linear scaling at lower wave numbers.However, in the context of a phase transition, the kinetic energy increases and reaches a constant over a certain time span.
In this appendix, we consider a simplified model for the time evolution of the kinetic energy and study its effect on the linearly rising part of the GW spectrum at low wave numbers.We assume that the kinetic energy evolves such that the quantity ∆ has the form  In figure 9, we show Ω GW /kΩ 0 for different values of k p , η d , and η.The upper panel shows the spectrum for k p = 10 and the lower panel corresponds to for k p = 30.Each panel is further divided: the left side shows Ω GW /kΩ 0 at η = 10, whereas the right-hand side corresponds to η = 100.In each plot, the red, blue, and black curves correspond to the case with no time evolution (η d → 0), η d = 1/(0.1kp ), and η d = 1/(0.03kp ), respectively.From this figure, it is evident that the time evolution affects the amplitude of the shallow part of the GW spectrum towards the peak, thus slightly reducing the slope.However, this effect reduces as we decrease the value of η d .

D Case of an non-uniform initial energy density
To consider the case of a non-uniform initial energy density such that ln ρ(k) = k • ũ/(c 2 s k), we have conducted new simulations A1 ′ and A3 ′ , analogous to the runs A1 and A3.In figure 10, we provide a comparison between these simulations.The red and black curves represent the uniform case, while the red and green curves correspond to the non-uniform case at simulation times t = 1.5 and t = 41, respectively.From this figure, it is evident that the ringing effect is reduced in the non-uniform case.Additionally, it is worth noting that the resulting GW spectra at t = 41 exhibit negligible differences between the two cases.Consequently, we conclude that initializing the energy density in accordance with the continuity equation is important, although its impact on the resulting GW spectrum is negligible.

2 Ω 2 Figure 1 .
Figure1.GW spectrum without expansion at different times.The amplitude of the spectrum in the range 1 < k < 10 is unchanged.However, it continues increasing at k = 30 with time ∝ η.By assuming that the GW spectra below the peak maintain the k 8 spectrum, the transition from k 8 to the flat regime moves toward the left with a speed proportional to η 1/8 .The transition to k 2 propagates ∝ η in the horizontal direction.

3 Figure 2 .
Figure 2. Scaled GW spectrum for different values of k p at time η = 10.In the left panel, the spectrum is normalized by k p and by k 2 p in the right panel.

Figure 3 .
Figure 3. Kinetic spectra (left) and GW spectra (right), obtained numerically from synthesized sound waves in a non-expanding background, as described in section 2.4.Red curves show the semianalytic spectra from the sound shell model, blue curves represent numerical spectra.Oscillations in the GW spectra at large k are an artifact of large time steps.

Figure 4 . 2 Figure 5 .
Figure 4. Kinetic, stress and mean GW spectra for irrotational runs (solid lines) and vortical runs (dashed lines).The red and black curves show spectra at η = 1.5 and η = 41.0,respectively.Note the presence of the bump in all irrotational runs, absent in all vortical runs.

Figure 6 .
Figure 6.Left: GW spectra for k p = 30 at times η = 2, 4, 6, and 8 for cases with expansion included.Right: time evolution of Ω(k p )/(k p Ω 0 ) for the same case, but for different values of k p .

Figure 8 .
Figure 8.In this figure, we show the evolution of Ω GW /(kΩ 0 ) with time for k = 2, 4, 6, and 8.In the right panel, we normalized the time with η k = 2π/k.

8 Figure 9 .
Figure 9. Compensated GW spectra, Ω GW /kΩ 0 , for non-expanding cases with different values of growth time η d , peak wave number k p , evaluated at conformal time η for η = 10 and 100 in the leftand right-hand panels.In the top (bottom) panels, we have k p = 10 (k p = 30).

Figure 10 .
Figure 10.Kinetic and normalized GW spectra at different times.The red and blue curves represent the spectra at t = 1.5, while the black and green curves correspond to t = 41.0.

Table 2 .
The spectral index of GW spectrum below the peak for the expanding case, computed in the sound shell model.