{'fraction_non_alphanumeric': 0.04153457682099675, 'fraction_numerical': 0.015297847242701269, 'mean_word_length': 4.293833544736072, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generation of gravitational waves from freely decaying turbulence

We study the stochastic gravitational wave background (SGWB) produced by freely decaying vortical turbulence in the early Universe. We thoroughly investigate the time correlation of the velocity field, and hence of the anisotropic stresses producing the gravitational waves. With hydrodynamical simulations, we show that the unequal time correlation function (UETC) of the Fourier components of the velocity field is Gaussian in the time difference, as predicted by the"sweeping"decorrelation model. We introduce a decorrelation model that can be extended to wavelengths around the integral scale of the flow. Supplemented with the evolution laws of the kinetic energy and of the integral scale, this provides a new model UETC of the turbulent velocity field consistent with the simulations. We discuss the UETC as a positive definite kernel, and propose to use the Gibbs kernel for the velocity UETC as a natural way to ensure positive definiteness of the SGWB. The SGWB is given by a 4-dimensional integration of the resulting anisotropic stress UETC with the gravitational wave Green's function. We perform this integration using a Monte Carlo algorithm based on importance sampling, and find that the result matches that of the simulations. Furthermore, the SGWB obtained from the numerical integration and from the simulations show close agreement with a model in which the source is constant in time and abruptly turns off after a few eddy turnover times. Based on this assumption, we provide an approximate analytical form for the SGWB spectrum and its scaling with the initial kinetic energy and integral scale. Finally, we use our model and numerical integration algorithm to show that including an initial growth phase for the turbulent flow heavily influences the spectral shape of the SGWB. This highlights the importance of a complete understanding of the turbulence generation mechanism.


Introduction
Gravitational wave (GW) signals from the early Universe have the potential to provide a new observational perspective on high energy physics phenomena. In this context, first order cosmological phase transitions provide a compelling source of GWs. This was first proposed many years ago (see Refs. [1][2][3][4][5][6]), before it became clear that the electroweak (EW) symmetry breaking proceeds as a crossover in the Standard Model [7]. However, it has since emerged that many scenarios beyond the Standard Model (BSM) lead to first order phase transitions at -and beyond -the EW scale, reopening the case for the study of GW production from first order phase transitions (for a GW-oriented review, see Ref. [8]). This is particularly interesting in the context of the Laser Interferometer Space Antenna (LISA), which is sensitive to a frequency window in the millihertz range [9]. In the context of primordial GWsourcing processes that are localized in time, such as a first order phase transition, the characteristic frequency of the GWs today can be connected to the characteristic time/length scale of the source anisotropic stresses L * , via f 16.5 · 10 −3 mHz 1 L * H * T * 100 GeV g * 100 Here a subscript * denotes the epoch at which the phase transition occurs. This shows that LISA can potentially detect the GW signal from transitions in the window 100 GeV -1 TeV, if the characteristic time/length scale of the anisotropic stresses is of the order L * ∼ 10 −2 to 10 −3 of the Hubble scale H −1 * at the time of the phase transition. These can be considered as typical values for L * H * , given that L * is related to the mean bubble separation R * (see e.g. [8,10] and references therein). LISA will therefore offer a new way to probe BSM physics, complementary to the Large Hadron Collider (see e.g. Ref. [11]).
The first 3-dimensional (3D) simulations of the coupled system of a scalar field and a relativistic fluid have shown that sound waves, produced in the fluid by expanding bubbles, are also a promising source of GWs [23][24][25]. These simulations have been performed using the code SCOTTS (Simulations of COsmic Thermal TransitionS) on which we will also rely in the present work. Analytical insight into the generation of GWs by sound waves has been developed in Refs. [26,27], and a new simulation technique for computing the GW spectra in the sound shell limit has been proposed in Ref. [28].
Refs. [23][24][25] performed simulations of first order phase transitions of weak to intermediate strength, e.g. α 0.1, where α is the ratio of the trace anomaly of the energy momentum tensor and the thermal energy. It was found that the bulk fluid motion consisted primarily of compressional modes, which correspond to sound waves. In this case, sound waves are expected to be the dominant GW source. However, shocks will develop in the fluid and are expected to convert the acoustic phase into a turbulent one [29,30]. The characteristic time of shock formation can be estimated as τ sh ∼ R * /v , where v denotes the root-mean-squared (rms) velocity of the acoustic motion. For stronger transitions, in which the initial value of v is typically larger, the onset of turbulence is expected to occur earlier and give a more substantial contribution to the total GW signal. Furthermore, recent simulations of transitions with α ∼ 0.5 found that, if the bubble wall velocity is subsonic, flows with substantial vorticity are already generated during the transition itself when fluid shells and bubble walls collide [31].
In the present work, we study the GW signal generated by a hypothetical turbulent phase in the aftermath of a first order phase transition. The first analyses of the GW signal from turbulence have relied on analytical modelling of the turbulent flow, and semi-analytical estimates of the GW signal [6,[32][33][34][35][36][37]. In particular, Ref. [38] evaluates the GW signal from all components (compressional, vortical, magnetic field) of both standard, and helical, freely-decaying magneto-hydrodynamic (MHD) turbulence. Various spectral shapes and scalings of the GW power spectrum per logarithmic wave number interval have been found, depending on the relative amplitude of the compressional and vortical components, and on the presence of helicity. In Ref. [38], it was also argued that the correct auto-correlation time to be used in the equations describing GW production by turbulence is the Eulerian eddy turnover time (see Ref. [39]). Far in the inertial range, the Eulerian eddy turnover time becomes τ E (k) ∼ (kv sw ) −1 , where v sw v/3 is a locally uniform velocity field sweeping the vortices, in accordance with the Kraichnan random sweeping approximation [40], and v is the rms velocity of the turbulence. Earlier work used the Lagrangian eddy turnover time τ (k) ∼ (kv ) −1 , where v denotes the characteristic velocity on the scale ∼ 1/k [33,34,36]. Our simulations confirm that the Eulerian eddy turnover time is the correct option.
Simulations of both non-helical and helical MHD turbulence and the resulting GW generation have been carried out in Refs. [41][42][43][44]. In these works there is an initial phase in which the MHD turbulence develops, starting from a nearly monochromatic electromotive force, or a kinetic forcing. They find that the spectral shape of the GW signal varies, depending on whether a fully-developed turbulent magnetic and/or velocity field is present in the initial conditions, or whether it is sourced by the forcing.
In this paper, we develop a semi-analytical model of GW generation by freely-decaying kinetic, vortical, non-helical turbulence, which we support with the results of simulations. We rely on the relativistic hydrodynamic code SCOTTS, developed by the authors of Refs. [23][24][25]31], but we use it to study the evolution of vortical fluid motion, excluding the dynamics of the scalar field undergoing the transition (see Section 5). We have chosen to over-simplify the turbulence model (no compressional modes, no helicity, no magnetic field, no initial forcing phase) in order to have full control of our analytical understanding of the GW production. Inclusion of these additional features will be tackled in future works.
The layout of the paper is as follows. We begin with a short review of stochastic gravitational wave backgrounds in the context of vortical turbulence in Section 2. We first introduce the equations of motion for the metric perturbations in an FLRW background, and write their solutions in terms of Green's functions and the anisotropic shear stress. The gravitational wave energy density power spectrum in the radiation era is then constructed, assuming a statistically homogeneous and isotropic signal. The solution for the gravitational wave energy density power spectrum is given in terms of an integral over the Green's function and the unequal time correlator (UETC) of the anisotropic stress. The section concludes by reviewing the computation of the anisotropic stress UETC in the context of vortical turbulence, decomposing the anisotropic stress UETC into a kernel which contains the fluid velocity UETC.
We continue by modelling freely decaying turbulence in Section 3. In Section 3.1, we introduce the velocity power spectrum for a divergence free velocity field, along with important quantities such as the integral scale of the flow ξ, and the rms velocity v. In our modelling, we assume that the turbulent flow is described by a von Kármán spectrum, which has a Kolmogorov power law in the inertial range and a causal slope of k 5 at large scales, where k refers to the power spectrum wavenumber. The evolution of the velocity field in decaying turbulence is considered in Section 3.2. Here we review the self-similarity properties of the turbulent power spectrum, starting from dimensional arguments and the assumption that there is a single length scale in the flow. On further assuming that the eddy turnover time, τ ξ , is the only principal timescale in the flow, one arrives at exponential decay laws for the kinetic energy and integral scale, where the decay exponents of both quantities are in principle related. We conclude with Section 3.3, where we introduce a model for the velocity UETC. To model the decorrelation in time, we follow Ref. [38] and adopt the exponential decorrelation provided by the short-time analysis of Ref. [39]. This consistently reproduces the random sweeping model [40] in the inertial range, which is supported by numerical works simulating decaying isotropic turbulence (see e.g. Ref. [45]).
In Section 4 we expand on the discussion of UETCs by considering them in the context of positive kernels and Mercer's condition. Indeed, any two-point correlation function of a random variable satisfies a set of inequalities known as Mercer's condition. Positive kernels are defined as the functions satisfying Mercer's condition, thus they are candidates to model UETCs. Earlier works, such as Ref. [36], encountered the question of how to define a valid correlation function for a non-stationary process. Here, we consistently address this issue in the context of the theory of positive kernels, and propose to model the unequal time correlator (UETC) of the turbulent velocity field as a Gibbs Kernel. This provides a way to symmetrize the unequal time correlator and, most importantly, it guarantees that the anisotropic stress two-point function is a valid correlation function. As a consequence, this leads to a GW energy density power spectrum which is always positive, as it should be.
After discussing the theory, we move on to describe our simulations in Section 5. We introduce the Minkowski space simulation code, SCOTTS, and discuss how it has been modified to study relativistic vortical turbulence. To prepare the initial conditions, we initialize the velocity field in momentum space, where each mode is randomly drawn using Gaussian statistics from a von Kármán power spectrum. We also project the velocity field so that it is purely vortical. After outlining our suite of simulation runs, we specify the methodology for outputting the velocity power spectra and UETC, as well as the procedure for evolving the metric perturbations and obtaining the gravitational wave power spectrum in the simulations.
In the following two sections we present our results. First, we focus on the evolution of the velocity field in the simulations, in Section 6. We study the evolution of the kinetic energy and of the integral scale of the flow in Section 6.1, and link it to the model of Section 3.2. In the simulations, the velocity power spectrum displays self-similar scaling towards the end of the simulation, but the theoretical relation between the decay law exponents and the velocity spectrum power law at large scales is not satisfied. We also test our model for the velocity UETC in Section 6.2. The combination of the Gibbs Kernel model of the UETC with the decorrelation of Ref. [39] performs well for both the inertial range and for large scale modes. We compare our model to that of Ref. [38], and find that our model has a closer agreement with the simulation data.
Our results for the gravitational wave power spectrum are given in Section 7. There are three different cases. We begin by evaluating the gravitational wave power spectrum under the assumption that the turbulence is stationary, in Section 7.1. Next we cover the case of instantaneously generated turbulence in Section 7.2, which contains our main results regarding the SGWB spectrum. Here our initial conditions are a fully developed turbulent power spectrum. We first show the SGWB extracted from the simulations. We then calculate the gravitational wave signal using our model of freely decaying turbulence outlined in Sections 3 and 4. In this case, to calculate the gravitational wave power spectrum, we use Monte Carlo integration with importance sampling to compute the integrals found in Section 2. Using this technique, we confirm that the precise values of the decay law exponents are unimportant for the final spectral shape. We also describe the constant source approximation, in which the source of gravitational waves is taken to be constant in time, before it is turned off abruptly [46]. We obtain in this context a ready-to-use formula for the SGWB spectrum from hydrodynamic turbulence. We compare the gravitational wave power spectrum under all of these methods, finding good agreement in their mutual regions of applicability. The final case we consider in Section 7.3 is the inclusion of a simplified turbulence growth phase. Using our turbulence model and Monte Carlo integration, we calculate the gravitational wave power spectrum in the case of linear growth for both a C 0 and a C 1 growth phase. The inclusion of a growth phase substantially affects the shape of the gravitational wave power spectrum, underlining the importance of studying the creation mechanism of turbulence in the early Universe. Finally, we summarize our results and conclusions in Section 8.
To help the reader, in Table 1 we summarize some timescales and temporal parameters which appear in different sections of this paper.
2 Stochastic background of gravitational waves from vortical turbulence

Generation of gravitational waves
In the cosmological context, GWs are described by transverse and traceless tensor spatial perturbations h ij of the background FLRW metric for a homogeneous and isotropic Universe:  Table 1. Important timescales and temporal parameters used in this article. A * subscript indicates quantities evaluated at τ * . In the above, k refers to the magnitude of the non-zero wavenumbers on the discrete momentum space lattice in the simulations.
Here a(τ ) is the scale factor, and τ the conformal time. It follows from the linearized Einstein equations that h ij (k, τ ) satisfies where H =ȧ/a is the comoving Hubble parameter, · = d/dτ , k denotes a comoving wave-vector, and Π is the transverse traceless part of the anisotropic stress tensor. The anisotropic stress is defined by with T µν the fluid energy momentum tensor, and p its pressure. The transverse traceless component of the anisotropic stress is extracted by means of the projector and We work in the radiation era, and therefore we further define Π (T T ) ij = (4ρ/3)Π ij with ρ the energy density in the early, radiation dominated Universe, and p = ρ/3. Using Friedmann's equation, Eq. (2.2) can be rewritten asḧ Assuming that there were no changes in the relativistic degrees of freedom, H = τ −1 , and changing variable from h ij (k, τ ) to τ h ij (k, τ ), Eq. (2.7) becomes the equation for a forced harmonic oscillator whose Green's function is known (see Ref. [47] for a review). The solution and its time derivative are where τ ini is the conformal time at which the source turns on.

GW energy density power spectrum
We assume that the GWs are generated by a stochastic process (specifically vortical turbulence) and consist of signals coming from patches in the sky that were causally disconnected at the moment of their emission in the very early Universe. Assuming statistical homogeneity and isotropy, the two-point correlation function ofḣ ij can be written as which defines the spectral density Pḣ. The fractional GW energy density power spectrum at conformal time τ in the radiation era is then . (2.11) The relevant statistical object of the GW source is the unequal time correlation function (UETC) of the anisotropic stress, defined by From Eq. (2.9), one then obtains 14) The GW energy density fraction Eq. (2.11) becomes then where we have used the Friedmann equation 3H 2 = 8a 2 πGρ. As can be appreciated from Eq. (2.15), the GW signal is sourced by the UETC of the anisotropic stress energy of the turbulent fluid, PΠ(k, τ, ζ).

The unequal time anisotropic stress correlator
We now review the computation of the UETC of the turbulent fluid anisotropic stress, defined in Eq. (2.12). The starting point is the spatial, off-diagonal part of the energy momentum tensor of a relativistic fluid. In order to simplify the computation, we neglect the spatial dependence of the fluid enthalpy density, and set the Lorentz factor γ = 1. This amounts to assuming non-relativistic turbulence. As we shall see, we find an excellent agreement between the semi-analytical results based on the non-relativistic source model and the relativistic simulations, also in what concerns the GW signal. Therefore, the assumption of non-relativistic turbulence does not significantly affect our final results, within the tested range of initial rms velocity. We then approximate We assume Gaussianity of the velocity field to rewrite the four-point velocity correlator in terms of two-point correlators using Wick's theorem: (2.18) Similarly to the anisotropic stress UETC of Eq. (2.12), we can define the turbulent velocity UETC.
Under the further assumption that the velocity field is divergence-free, statistically isotropic and homogeneous, the velocity UETC is Then on using Λ ij m (k)Λ ijrs (k) = Λ mrs (k), Eq. (2.18) becomes (2.20) where q = k − p (note that the first term in Eq. (2.18) does not contribute). The contraction over indices can be carried out explicitly using Eqs. (2.4) to (2.6), We can now extract the kernel of the anisotropic stress defined in Eq. (2.12) and obtain This result agrees with Ref. [36], but differs from Ref. [38]. The reason is that in Ref. [38], the tensor indices on the two terms in square brackets on the right-hand side of Eq. (2.20) were assumed identical, but this is not the case as can be seen explicitly in Eq. (2.21) and Eq. (2.22). The following Sections 3 and 4 are devoted to the construction of a model for the velocity UETC P v (p, τ, ζ) in Eq. (2.24).

Modelling freely decaying turbulence
In this section, we provide definitions and review properties of freely decaying turbulence. This will be complemented with new theoretical input regarding positive kernels in Section 4, leading to the new velocity UETC model which we propose in this work. Our new model will be validated with relativistic, simulations in Section 6.

Velocity power spectrum
The real-space two-point correlation function of a statistically homogeneous, isotropic and divergencefree velocity field reads [48] b is the rms velocity squared, and the functions Σ t and Σ represent respectively the transverse and longitudinal correlation functions. Since the velocity field is divergencefree, we have that ∂b ij /∂r i = 0 and the correlation functions are related through The spectral density In spherical coordinates and using Eq. (3.2), this becomes We assume that the correlation function of Eq. (3.1) has compact support because it vanishes outside the causal horizon. As a consequence, the leading order term in Eq. (3.6) is r 3 Σ ∞ 0 = 0, and the spectral density on large scales is given by Therefore, causality in the early Universe points to Batchelor turbulence, at least on super-horizon scales [32]. The coefficient of the amplitude of the large-scale part of the spectrum in Eq. (3.7) is the Loitsyansky integral dr r 4 Σ (r). This quantity was conjectured to be constant in time in the absence of long-range correlations (see e.g. Ref. [48]), which would then be consistent with the causality requirement. The power spectrum P v (k) encodes how v 2 , or equivalently the non-relativistic kinetic energy per unit enthalpy, is distributed into the different length-scales: (3.8) In our notation, it is related to the spectral density through Note that the power spectrum as defined here is related to the energy spectrum E(k, τ ) as defined in Ref. [48] through P v (k, τ ) = 2kE(k, τ ). (3.10) To characterize the typical length-scale of the system, we also define the integral scale ξ in terms of the longitudinal correlation function, or equivalently in terms of the power spectrum We use the integral scale to define dimensionless, time-dependent wavenumbers with A a normalization constant given in Eq. (3.14) below.
Concerning the shape of the velocity power spectrum, we assume in our modelling that in the inertial range it is determined by the well-known Kolmogorov K −2/3 law [49], and on large scales by the causal K 5 slope, motivated by Eq. (3.7). These properties are captured by the von Kármán spectrum [36,38,48,50] The coefficients A and B are chosen to ensure that the definitions of the kinetic energy and the integral scale of Eqs. (3.8) and (3.11) are consistent 1 : √ πΓ(17/6) ≈ 4.02 (3.14) Eq. (3.13) assumes that the von Kármán spectral shape is preserved during the time evolution of the turbulent field. This power spectrum will be used for the semi-analytical evaluation of the GW signal (see Section 7.2.2): self-similarity in time is therefore an underlying assumption of the resulting GW spectra. Furthermore, Eq. (3.13) will also be used as the initial condition for our hydrodynamical simulations. The simulations will show, on the other hand, that the von Kármán spectral shape is not fully preserved by time evolution: while for K 1, a power law consistent with K −2/3 returns after O(10) eddy turnover times, a shallower one than K 5 emerges for K 1 (see Section 6.1). Nevertheless, the GW spectra output by the simulations, and those resulting from the semi-analytical integration, are in very good agreement in their common wavenumber K region: we therefore anticipate that differences in the time evolution of the large-scale part of the turbulent spectrum do not significantly affect the final results. This is due to the fact that the production of gravitational waves is fast, as we shall see.
As a consequence of the fast GW sourcing, though, the particular choice of the initial velocity power spectrum determines the final GW spectral shape. Equation (3.13) is the only input we use for exploring the SGWB signal from decaying turbulence. Nevertheless, we are able to draw general lessons from our analysis.

Evolution of the velocity field in decaying turbulence
In this section we provide a very simplified review of the time evolution of the freely-decaying turbulent power spectrum, and of average quantities such as the kinetic energy per unit enthalpy v 2 (τ ) and the integral scale ξ(τ ).
In line with the model of Section 3.1, we assume that there is only one principal length scale in the flow, ξ (see e.g. Ref. [52] for a turbulence description not satisfying this hypothesis). Since the power spectrum is dimensionless, it must generally take the form (consistently with Eq. (3.13)) where the function ϕ(kξ(τ )) must satisfy integrability conditions following Eqs. (3.8) and (3.11), and * indicates quantities evaluated at the start of turbulence decay. At large scales, the absence of relevant dynamical processes implies a simple power law spectrum where, for the moment, the spectral index β is left unspecified. Further assuming that at large scales the power spectrum is constant in time 2 , leads to v 2 (τ )ξ(τ ) 1+β = const .
The individual time evolution of the rms velocity and the integral scale can be derived, if we further assume that there is only one principal timescale in the flow, the eddy turnover time τ ξ = ξ/v. With this assumption, one can write where c v is a constant. Together with Eq.
where the decay exponents p and q take the values 24) and N e is a constant representing the number of initial eddy turn-over times that the kinetic energy takes to decay, and the integral scale to grow. In particular p and q satisfy p = (1 + β)q = 2(1 − q) . (3.25) Note that this model does not predict the value of β, although consistency with the arguments of Section 3.1 would indicate β = 4 as would follow from the constancy of the Loitsyansky integral. In Section 6.1 we will report on our numerical investigation of the evolution of the velocity power spectrum in the framework of the model described above. The validity of the decay laws given in Eqs. (3.21) and (3.22), together with the actual value of the parameter β (and consequently of the decay law exponents through Eq. (3.25)), are currently open problems of the theory of turbulence. Invariance of the Navier-Stokes equations under appropriate rescaling leads to the self-similarity of the velocity power spectrum, and links the decay law exponents p and q to the scaling parameter [59,60]. The scaling parameter can then be related to the slope of the initial velocity power spectrum at large scales, if one assumes that the latter also satisfies the self-similarity conditions (see e.g. Ref. [61]). In this context, an initial power spectrum leading to β = 4, as motivated in Section 3.1, provides p and q values consistent with the classical Kolmogorov theory, namely p = 10/7 and q = 2/7 [48,62] (see also [63] for a numerical result). The relation between the scaling parameter and the slope of the initial power spectrum at large scales has also been derived in the context of renormalization group analyses of the incompressible MHD equations, as holding at the fixed point [64,65]. However, the initial velocity power spectrum does not need to necessarily satisfy either the self-similarity condition, nor does it need to exhibit power law behaviour. Relaxing these assumptions changes the prediction on the subsequent free decay model, disconnecting the value of β from the initial conditions [56,58]. Recently, new models of the turbulent free-decay have emerged, following which the turbulence decay is governed by the integral scale rather than by the asymptotic behaviour at large scales [52,[66][67][68][69]. According to these findings there would be then no universal decay regime, since the details of the spectral shape near the peak are related to the turbulence production mechanism.
Note that in both this Section and Section 3.1 we have neglected the dissipative scales. The turbulent kinetic energy dissipates at the viscous scale λ, at which the Reynolds number is Re(λ) = v λ λ a/ν = 1, where ν denotes the kinematic viscosity. The viscous scale λ increases with time, and the turbulent velocity on a given scale 2π/k goes to zero once the dissipation scale has grown to λ(τ ) = 2π/k. Furthermore, λ(τ ) grows faster than the integral scale ξ(τ ): the turbulent inertial range is therefore fully dissipated when the two scales become equal, so that Re(ξ) = 1. As previously demonstrated [36,37], turbulence is expected to free-decay for several Hubble times, before it is fully dissipated. For example, with the initial set of parameters T * = 100 GeV, v * = 0.03, ξ * H * = 0.01, q = 2/7, p = 10/7, one obtains Re(ξ) = 1 at T fin ∼ 0.3 GeV, corresponding to about 1000 τ ξ * . The kinetic viscosity in the early Universe is indeed extremely small: for the same initial conditions, one finds Re(ξ * ) ∼ 10 11 , and the inertial range spans about 8 orders of magnitude in wavenumber. In principle, we should have inserted a cutoff at k > 2π/λ(τ ) in the turbulent power spectrum given in Eq. (3.13), and we should account, in the GW production, for the viscous dissipation of the kinetic energy at each scale k in addition to the overall free decay described by Eqs. (3.21) and (3.22). However, as we will discuss in Section 7.2, the bulk of the GW signal is generated faster than the turbulent decay: therefore, the dynamics at the viscous scale plays no role in practice in the present analysis, and it is justified to neglect it altogether.

Unequal time correlator of the velocity field
Together with the overall free-decay, we also need to model the time decorrelation of the turbulent velocity field. This is encapsulated by the normalized velocity UETC 26) The first model to have been developed was Kraichnan's random sweeping approximation [40], in which the decorrelation function in the inertial range is a Gaussian in the variable kv|τ − ζ|, for small time differences. The thinking behind the model is that modes are decorrelated by being "swept" by the large-scale flow, whose rms velocity is v. This model has good numerical support in the inertial range for both stationary and decaying turbulence (see for example Refs. [45,70,71]). We review it in Section 3.3.1. However, for GW production calculations, we need to model the decorrelation around the peak of the velocity power spectrum. To this end, we study -for the first time with hydrodynamical simulations -the decorrelation of the turbulent velocity field outside the inertial range, as far as possible in the low wavenumber part of the spectrum. The simulation results are presented in Section 6.2, and show that a good model is obtained by replacing the rms velocity in the Gaussian with a time and scale dependent average decorrelation velocity, based on the sweeping velocity proposed in Ref. [39], which we present in Section 3.3.2.
Furthermore, we also need to take into account the formal requirement that the UETCs (velocity and anisotropic stress) must be positive definite kernels, which is not manifest in the simple Gaussian model. This will be presented in Section 4.

Kraichnan random sweeping approximation
The random sweeping approximation provides an explanation for the origin of the Gaussian functional form for the decorrelation. Kraichnan's sweeping hypothesis is based on the assumption that Fourier modes of the velocity field in the inertial range are advected without deformation by a locally uniform velocity field V, which may be time-dependent [40]. This amounts to neglecting the non-linear terms and writing ∂v ∂τ Since V is locally uniform, the evolution of modes with different wavenumbers k is decoupled. The turbulent velocity field can then be explicitly integrated: from Eq.
where ζ < τ . Under the assumption that the locally uniform velocity field V is statistically independent of the turbulent velocity field v at the initial time, the unequal time correlator of v can be expressed as The average of the exponential can be calculated assuming that its argument is a Gaussian random variable. It then follows that the velocity field decorrelates with a characteristic Gaussian law 3 , This expression can be further simplified assuming that on average the locally uniform velocity has the same amplitude in the three spatial directions [40], leading to where v sw (τ, ζ) is the average sweeping velocity which takes the form [45,72] v 2 sw (τ, ζ) = Kraichnan's model assumes stationary turbulence and sets the sweeping velocity to the rms velocity divided by √ 3, under the hypothesis of statistical isotropy [40,73]. Since this model has good numerical support in the inertial range [45,70,71], The sweeping velocity should reduce to at equal time and in the inertial range, with C 2 v a numerical factor whose value will be inferred from simulations. As presented in Section 6.2, simulations indicate C 2 v 1,mserving as a verification of the random sweeping approximation and of statistical isotropy.

Scale dependent sweeping velocity
The random sweeping approximation reviewed in Section 3.3.1 applies in the inertial range. However, in order to evaluate the GW generation, we need to model the time decorrelation of the velocity 3 We recall that for a Gaussian random variable, field on all scales (see for example Eq. (2.24)). The decorrelation dynamics at scales larger than the integral scale in freely decaying turbulence have received little attention: numerical studies such as those performed in Refs. [45,72], for example, study decorrelation only in the inertial range. The largest scale analysed in Ref. [70] is k = 1, corresponding to the scale of the forcing, namely the peak of the power spectrum: it appears from this analysis that even this scale decorrelates more slowly than those of the inertial range.
In the context of the literature dedicated to GW production by turbulence, various decorrelation models have been proposed, see Refs. [33,36,38]. Ref. [36] assumed that the large scales do not decorrelate. The only time dependence of the velocity spectrum outside the inertial range was therefore due to the free decay, and an exponential decorrelation was inserted for wavenumbers in the inertial range by means of a step function. This introduced a non-physical discontinuity (see Eq. (57) of [36]). Furthermore, as pointed out in Ref. [38], Refs. [33,36] used the Lagrangian eddy turnover time as typical decorrelation time, which is scale-dependent in the inertial range, and therefore in conflict with the numerical results of Refs. [45,70].
Ref. [38] recommended the use of the model of Ref. [39], which improves Kraichnan's random sweeping approximation by taking interactions between modes into account. The result is an Eulerian eddy turnover time which is weakly scale-dependent for wavenumbers up to O 10 2 times that corresponding to the integral scale. We adopt the same form for the Eulerian eddy turnover time to construct our model of the velocity UETC, but correct it to ensure the mathematical consistency of the correlation, as discussed in Section 4.
Following Refs. [39,74], we define an instantaneous scale-dependent decorrelation velocity where the function h(x) is given by 4 Note that the calculations performed in Ref. [39] in principle apply only to the inertial range. We will see that extending it to lower wavenumbers is justified by the numerical data. Eq. (3.36) reduces to v 2 dc (τ ) v(τ ) 2 /3 far in the inertial range, consistent with the result of the random sweeping hypothesis, and gives C 2 v = 1 (see Eq. (3.35)). However, it also has another term, relevant at scales closer to the integral scale ξ, which is proportional to −v 2 k ∼ −P v (k). This negative contribution comes from the analysis at small time intervals and appears to slow down the process of decorrelation. In order to account for decorrelation at all scales, one can therefore attempt to use Eq. (3.36) as an extension of the sweeping velocity [38]. Pushing Eq. (3.36) to large scales, one can appreciate that v dc (k, τ ) provides a continuous interpolation from v 2 /3 in the inertial range to 2v 2 /15 at large scales, as shown in Fig. 1. In the semi-analytical GW evaluation, we will solve the integral of Eq. (3.36) numerically, although a good analytical fit to that we compare to the numerical result in Fig. 1 (note that this fit differs slightly from what was given in Ref. [38]). In order to describe how the relevant scales in the velocity spectrum decorrelate, one might substitute, in the exponential of Eq. (3.33), the decorrelation velocity given in Eq. (3.36). This is the approach adopted in Ref. [38]. However, this breaks the time symmetry of the UETC P v (k, τ, ζ) (see Eq. (2.19)). To restore this symmetry, the decorrelation velocity entering the exponential in Eq. (3.33) must depend on both times. The next section is devoted to ensuring the symmetry and positivity requirements that are necessary to model a valid UETC.

Unequal time correlators as positive kernels
The two-point correlation function of any random variable must be a positive kernel [75], in other words it must satisfy certain properties -most notably Mercer's theorem. This applies to the turbulent velocity field UETC Eq. (2.19), as well as to PΠ(k, τ, ζ) Eq. (2.12), describing the random anisotropic stresses arising from the turbulent field. The model of the velocity field -and consequently of the anisotropic stress -UETCs that we construct in this section ensures that they are positive kernels. According to Mercer's theorem, this implies that the GW energy density power spectrum is indeed also positive, see Eq. (2.15). This issue had already been raised in Ref. [36], but it is analysed here in more depth directly on the velocity field UETC.

The unequal time correlation function as a non-stationary Gaussian kernel
As we have discussed in Section 3, the turbulent rms velocity and integral scale are evolving in time, while the turbulent source decorrelates. GWs are therefore generated by a non-stationary, random process, whose UETC we need to model. To do so, we make use of some general properties of two-point correlators of stochastic processes, which we review in the following.
The two-point correlator of a complex valued stochastic process φ(τ ) is defined by We assume that the real and imaginary part of the stochastic process are uncorrelated, in which case K(τ 1 , τ 2 ) is real. K(τ 1 , τ 2 ) is a positive semi-definite kernel which satisfies several properties: symmetry K(τ 1 , τ 2 ) = K(τ 2 , τ 1 ); the Cauchy-Schwartz inequality; and the Mercer condition (see e.g. Ref. [75]) for any interval I. This follows from the fact that, for any well-behaved function f , the stochastic variable X = dτ f (τ )φ(τ ) must have positive variance |X| 2 ≥ 0. One can show that the linear combination and multiplication of two kernels yield another kernel [76].
The velocity UETC P v (k, τ, ζ) is the two-point correlation function of the stochastic process v i (k, τ ) (see Eq. (2.19)), and is therefore a positive semi-definite kernel. In Section 3.3 we analysed the time evolution and decorrelation properties of the velocity field, based on our understanding of the turbulence physical process. In this section, we explain how the fact that a UETC must be a positive semi-definite kernel leads to our choice of the decorrelation velocity at unequal times v dc (k, τ, ζ), and how to symmetrize P v (k, τ, ζ).
Most importantly for us, stationary kernels also include the Gaussian kernel with α > 0, which is the form provided by the Kraichnan sweeping scenario Eq. (3.33). However, here we would like to model the power spectrum of freely decaying turbulence, which is, by definition, non-stationary. A simple but well-defined departure from stationary kernels is provided by locally stationary kernels, of the form where K 1 is a non-negative function and K 2 a positive stationary kernel [78]. The UETC of a freely-decaying turbulent velocity field, though, cannot assume this form. This is because the sweeping velocity depends explicitly on time, as can be seen from Eq. (3.33). For the problem at hand, we therefore also need to go beyond local stationarity. We do so by introducing process-convolution kernels [79]. This specific class of non-stationary kernels is defined through a two-valued complex function g(τ, u) as follows: It is easy to show that process-convolution kernels satisfy Mercer's condition Eq. (4.2). The processconvolution kernel adapted to our case is the non-stationary Gaussian kernel, defined with the function where J(τ ) is a free function.

Symmetrized unequal time correlator of the velocity field
To define the kernel representing the velocity unequal time correlator, we choose for a given k where v dc (k, τ ) is given in Eq. (3.36). Inserting this definition into Eq. (4.4) and further into Eq. (4.3), one finds a family of non-stationary Gaussian kernels .

(4.6)
Note that these kernels satisfy K k (τ, τ ) = 1, and are therefore of a suitable form to be used as a normalized UETC R(k, τ, ζ).
We can now use this result to construct a UETC P v (k, τ, ζ) which satisfies the Mercer condition. First we note that the UETC evaluated at equal times is the power spectrum, P v (k, τ, τ ) = P v (k, τ ). Hence, multiplying the kernel Eq. (4.6) by P v (k, τ )P v (k, ζ) will give a function of the two times which correctly reduces to the power spectrum at equal times. This function is also a kernel, as is a separable kernel, and the product of two kernels is a kernel. Our model is then where the average decorrelation velocity v dc (k, τ, ζ) is given by the harmonic average of the equal-time decorrelation velocity: .
Note that this is symmetric in the two times τ and ζ.
To summarize, we have shown that the velocity UETC of Eq. (4.7) is a positive kernel, i.e. a well-defined two point correlator. This guarantees that the anisotropic stress UETC PΠ(k, τ, ζ) of Eq. (2.24) is also a positive kernel, being the product of two positive kernels. Finally, the integrand of Eq. (2.15) is also a positive kernel. Indeed, G is a separable kernel since it results from the square of Eq. (2.9). As a consequence, the GW energy density power spectrum defined in Eq. (2.15) is always positive, as it should be.
Having managed to model the decorrelation of the freely-decaying turbulence with Eq. (4.7), we can go beyond the approach of Ref. [36]. In that paper the shear stress decorrelation was modelled directly in the anisotropic stress UETC as a top-hat function of the kind Θ[k|τ − ζ| − x c ]. We remark that such a function is not a positive kernel, since it is the Fourier transform of a sinc function: in order to derive a positive definite GW power spectrum, one must indeed take x c < π [36,37].

Simulations
We carry out a series of direct hydrodynamical simulations in order to study the freely-decaying turbulent flow, in particular the UETCs of the velocity field and its overall time evolution, as well as the SGWB signal produced. For our simulations, we use a modified version of the relativistic hydrodynamics code SCOTTS, previously developed to study the coupled evolution of the scalar field and the relativistic fluid, and the production of GWs during a thermal phase transition [23][24][25]31].
Here we are mainly interested in the dynamics of the fluid after the phase transition has completed. We therefore initialize the scalar field to lie in the broken phase at the start of the simulation. Given the form of the effective potential and equation of state in Ref. [31], this leaves us solely with the evolution of the relativistic fluid.
Ideally, one would want to simulate the coupled field-fluid model [80], in order to observe the generation of turbulence at all stages from the bubble collisions onwards [31]. In this work, however, we do not analyse the complete system up to fully developed turbulence. We rather concentrate on the decay of turbulence, inserting the turbulence parameters (e.g. kinetic energy and integral scale) as input values, disconnected from the rest of the phase transition dynamics. Here, we will initialize only the vortical component of the fluid velocity. This permits a cleaner analysis of the results and comparison to the semi-analytical modelling of Sections 3 and 4. We defer studying the compressional component and simulating the turbulence generation to future work. See Fig. 2 for a visualisation of the velocity field in a slice through one of our simulations at initialisation and after many eddy turnover times have passed.
In the rest of this section, we specify the equations of motion of the fluid that are solved by the code, and then explain how we fix the initial conditions and evaluate the UETCs. Finally, we discuss how we compute GW power spectra from the simulations.

Hydrodynamic evolution
Given that we are discarding the scalar field dynamics, the energy momentum tensor of the system becomes that of a perfect fluid, with the internal energy density in the fluid, p the pressure, U µ = γ(1, v) the fluid four-velocity and with Lorentz factor The equation of state is that of an ultrarelativistic gas with p = /3. It has been shown that the hydrodynamic equations of motion of an ultrarelativistic fluid in an expanding universe expressed in comoving coordinates and conformal time are the same as the hydrodynamic equations in Minkowski space-time, provided the dynamical quantities are replaced by appropriately scaled variables [81]. However, the source term for the GWs does not scale in the same way. Therefore, the Minkowski space-time code SCOTTS can be used to study hydrodynamic turbulence in an expanding background, but the output GW power spectrum can be matched to the expanding background result only for sources which are short-lived compared with the Hubble time [24].
The dynamical quantities that the code evolves are the fluid energy density E = γ with equation of motion [24] and the fluid momentum density Z i = γ 2 ( + p)v i , the components of which evolve according tȯ The evolution algorithm follows the approach taken in Ref. [82] with a leapfrog and operator splitting approach to updating the dynamical quantities. For this work we use a van Leer scheme for the advection update [83,84], whereas earlier uses of the SCOTTS code used an upwind donor cell scheme. For all of our simulations we pick the lattice spacing and timestep such that their ratio is dx/dτ = 2.5 5 . We denote the starting time of our simulations as τ = τ * . Note that while our simulations do not contain a physical viscosity, the introduction of a lattice spacing dx effectively sets a dissipation scale. This is due to the numerical viscosity introduced by the discretization of the fluid equations. Therefore, simulations with larger ξ * /dx have a larger effective Reynolds number.

Initial conditions of the simulation
In our simulations the velocity field v is not necessarily divergence-free. However, as outlined in the rest of this section, we initialize the velocity field to be divergence-free in our simulations. In Appendix F we show that in this case the longitudinal component of the velocity field remains negligible throughout our simulations, and so we treat v as divergence-free in subsequent sections. At the start of our simulations, we initialize the fluid energy density E everywhere. For our simulations we choose to set E dx 4 = 1/16. Next we initialize the velocity field of the simulations  Table 2. List of the simulations used in this work. We list the number of sites in the lattice for each simulation L 3 , the initial rms velocity v * = v 2 (τ * ) (see Eq. in Fourier space. We do this indirectly by initializing the spatial components of the four-velocity U = γv, such that where x and k refer to the discrete real and momentum space coordinates on the lattice, and L 3 is the total number of lattice sites. Each mode U k is randomly distributed and follows Gaussian statistics with mean zero and variance determined by the desired power spectrum serving as initial condition. Initializing U k rather than v k prevents individual modes of v k being drawn from the tails of the Gaussian with unphysical velocities larger than one. The U k field is then projected onto its vortical component with the projector ⊥ ij (k), see Eq. (2.6). Since U is real, we impose that U k = U * −k . Once the U k field is correctly initialized we perform a Fourier transform to find the real space configuration of the U field, and from this we can calculate the velocity field v.
This method allows us to initialize the simulation with arbitrary choices for the U field power spectrum P U (k, τ ) = k 3 P U (k, τ )/π 2 . Here P U (k, τ ) corresponds to the spectral density of the U field, defined via the equal time correlator (see also Eq. (3.3)) where in our initial conditions we have set the longitudinal component of U to zero. We initialize the simulations with a relativistic version of the von Kármán spectrum Eq. (3.13), supplemented with a high-wavenumber cutoff: (5.6) The two constantsk and C determine respectively the peak location and amplitude. The values of C andk are chosen so that the initial rms velocity v * = v 2 (τ * ) (defined in Eq. (3.8)) and integral scale ξ * = ξ(τ * ) (defined in Eq. (3.11)) extracted from the simulation will be close to a desired value. A list of parameters for the simulations we perform is given in Table 2. In Eq. (5.6) we insert an exponential cutoff to ensure that the power does not extend to the lattice scale, with k max dx = π/4 used in our simulations.

Velocity power spectra in the simulations
During the simulation we regularly output the rotational velocity power spectra P v (k, τ ). We validate that the longitudinal piece of the velocity correlator P v (k, τ ) is negligible over the full duration of all our simulations, providing it is initialized to be zero, see Appendix F. To find the velocity power spectrum within our simulations, we first perform a discrete Fourier transform on the real space field to obtain the Fourier modes The rotational power spectrum is then approximated by projecting and binning these modes in momentum space. The power spectrum for the r-th bin is constructed via where k r = 2π(r + 1/2)/Ldx is the central value of the bin in k, and ∆k = 2π/Ldx is the bin width. Note that the power spectra from the simulation are formed from an average over a spherical shell in k of Fourier modes.

Unequal time correlators in the simulations
We are also interested in studying the UETCs of the velocity field in these simulations as these enable us to model the GW source, see Eqs. (2.15) and (2.24). However, evaluating P v (k, τ, ζ) for many values of τ and ζ is very costly within a simulation, in part because each pair of τ and ζ represents snapshots of the field v that must be stored concurrently in memory. Instead, we define a reference time τ ref at which we store in memory the field v(k, τ ref ), and then compute the UETC at regular intervals with the field on the current timestep, i.e.
We evaluate the unequal-time correlator at times τ = τ ref + n∆τ UETC , where ∆τ UETC is the interval between UETC outputs (see Table 2) and n is some positive integer. Using a single value of τ ref per simulation is in principle a limitation of our analysis. Nonetheless, we are confident that the latter still captures the general behaviour of the time decorrelation of the turbulent field, especially given the good agreement we find between the simulations results and the theoretical model, see Section 6.

Gravitational waves
We also output GW power spectra calculated from our simulations. As SCOTTS is a Minkowski space code, we neglect cosmological expansion in the following. To find the metric perturbations h ij , we evolve an auxiliary tensor u ij according to Ref. [85] (see also Eq. (2.2)) The metric perturbations can then be recovered from our auxiliary tensor by performing a projection in momentum space, where Λ ij, m (k) is the projector onto transverse-traceless rank two tensors, Eq. (2.5). We can then define the power spectral density of the metric perturbations, Pḣ(k, τ ) via 12) and from this construct the power spectrum of the metric perturbations, From Pḣ(k, τ ) and a suitable scaling, one can then reconstruct the gravitational wave power spectrum (see Eq. (2.11)) dΩ gw dln k τ = Pḣ(k, τ ) 8πGρ * , (5.14) where ρ * is the conserved average total energy density (note that we insert a * subscript to make the connection with the expanding Universe quantity, see Eq. (5.17)). Dimensionful quantities in our simulations are always output normalized according to a corresponding length scale, which we have chosen to be the initial integral scale ξ * (defined later in Eq. (3.11)). Noting that the gravitational constant G has been scaled out of our equation of motion in Eq. (5.10), the output from our code is In this expression, the power spectra of the metric perturbations in the simulations are binned and constructed analogously to the velocity power spectra, The GW power spectrum is thus recoverable from the simulations by fixing a value for ξ 2 * /G. While the GW evaluation is carried out in flat space-time, we can still infer a value for the Hubble rate via the Friedmann equation This allows us to rewrite the scaling for the gravitational wave power spectrum as (note that ρ * is a conserved quantity throughout the simulation) The SGWB output of the simulations will be compared with the semi-analytic results in Section 7.2.4. Since the simulations neglect the expansion of the Universe, one can only compare the spectral region k > H * . The wavenumbers accessible in the simulations satisfy kξ * > 2πξ * /(L dx) > 0.03 (see Table 2), and we will therefore compare with semi-analytical results calculated setting H * ξ * = 0.001 (see Section 7.2.2).
6 Results: evolution of the velocity field In this section, we verify the model of decaying turbulence presented in Section 3 with the hydrodynamic simulations. We compare the time evolution laws of the turbulent kinetic energy and correlation scale with the simulation results in Section 6.1. Then in Section 6.2, we validate the symmetrized velocity field decorrelation function constructed in Sections 3.3 and 4, and thereby show that it can also describe decorrelation at large scales, outside the inertial range. Note that the turbulence model we adopted has been developed in the context of non-relativistic turbulence. The highest value of the initial rms velocity that we have tested in the simulations is v * 0.3 (see Table 2), finding good agreement with the predictions of the non-relativistic turbulence modelling, as far as the UETC, free decay evolution and SGWB spectra are concerned.

Velocity power spectrum and averaged quantities
Here we use the results from simulation (A) ( Table 2) to investigate the time evolution of the velocity power spectrum, and the values of the parameters β, p and q, introduced in Section 3.2.
In Fig. 3 we show the time evolution of the velocity power spectrum. Fig. 3a shows the velocity power spectrum in its raw form, whilst Fig. 3b shows the power spectrum rescaled with Eq. (3.16), indicating that there is indeed only one principal length scale in the flow, given our initial condition. The self-similarity is established starting from about 40τ ξ * , the time to which the first power spectrum line in Fig. 3b corresponds. Note that the large-scale slope of the late-time power spectrum is   significantly less steep than the k 5 corresponding to the Batchelor power spectrum imposed in the initial condition Eq. (3.13): therefore, the large-scale power spectrum is not constant in time in the simulation, before about 40τ ξ * . While the large-scale slope could be subject to finite volume effects, its constancy in time seems to be well established at late times.
Turning to the form of the velocity power spectrum at wavenumbers higher than the peak, we see from Fig. 3 evidence for a power law in the range 1 kξ * 5. This power law is approximately consistent with a Kolmogorov power spectrum (k −2/3 ) at late times. Although this was also the power spectrum in the initial conditions, one can see that this power law is lost in the early evolution, but later regained. We do not have the inertial range to give a firm value for the late-time power law index. For wavenumbers above kξ * 5, we see a small build up in power followed by a sharp fall-off for kξ * 10. The fall-off is set by the viscous damping scale of the numerical scheme, while the build up seems to be an associated feature of the van Leer advection scheme 6 .
In Fig. 4a we plot the combination v 2 ξ 1+β , predicted to be constant when the large-scale power spectrum is constant as well, see Eq. (3.18). It can be seen that this combination indeed remains approximately constant after about 40τ ξ * , provided β = 3. Consequently, in Fig. 4b we analyse again the self similarity of the power spectrum, this time rescaled according to Eq. (3.19) setting β = 3. Again, we only plot the rescaled spectra for τ −τ * 40 τ ξ * . It appears that the self-similarity condition Eq. (3.19) is satisfied by the simulation with β = 3. This result is consistent with the findings of Ref. [58]. However, the large-scale power spectrum is less steep than k 4 , which is the power law which would be inferred setting β = 3.
Simulation (A), therefore, seems broadly consistent with the self-similar behaviour outlined in Eq. (3.19), at late times. We have seen that β is not set by the slope of the power spectrum at low wavenumbers in the initial condition; rather, it seems to be part of the dynamics of self-similarity. We would need even larger lattices to study in detail the evolution of the power spectrum at low wavenumbers, where there appear to be correlations established at scales greater than the integral scale. Limitations of the size of the computational domain are known to possibly affect the interpretation of the turbulent decay features [66]. In (a) Evolution of v 2 ξ 1+β .  and provide a good fit to the simulation outcome. We shall see later in Section 7.2.2 that the precise values of the decay exponents do not have a significant effect on the final SGWB spectrum. Instead, we will see in Section 7.3 that the initial phase preceding the generation of the turbulence plays a much larger role in determining the spectral shape. We further analyse the behaviour of the power law exponents p and q of the kinetic energy and integral scale in Fig. 6. In Fig. 6a, we adopt the same approach as Ref. [58] and plot the state of the system in (p, q) space for a range of times during the simulation. The values of p and q are found by solving Eqs.     Fig. 6a. The self-similarity lines p = (1 + β)q are shown for various values of β as dashed lines. As in clear from the figure, the point (p, q) is still evolving at the end of the simulation. However, the evolution is moving in the direction of the intersection of the self-similarity line p = (1 + β)q with β = 3, and the scale invariance line, p = 2(1 − q). In Fig. 6b, we display the evolution of p and q against time, and also show the values that can be predicted from requiring scale-invariance and self-similarity for β = 3, namely p = 4/3 and q = 1/3.

Unequal time velocity correlations
In this subsection we report on our results for the decorrelation of the velocity field at different times. We introduce a new model for the unequal time correlator and test it against the numerical simulations.
Our model for the normalized UETC is obtained by combining Eqs. (3.33), (3.36) and (4.8) giving The factor in front of the exponential guarantees that the shear stress UETC is a positive definite kernel and hence that the resulting gravitational wave power spectrum is positive definite, as explained in Section 4.1. This is to be compared with the model of Ref. [38], In Fig. 7, we compare the two models to data from simulation (A), by plotting combinations of data and model functions which should be of Gaussian form. Results from the other simulations are given in Appendix D. The construction of the numerical UETC, P v (k, τ, ζ) is described in Section 5.4; in our simulations we measure the correlation of the velocity field at times ζ = τ ref and τ .  For our model, the data is plotted as the quantity . We choose a range of values for the wavenumber k, both larger and smaller than the inverse integral scale ξ * . The reference time is taken to be approximately one initial eddy turnover time, τ ref τ ξ * = ξ * /v * (see Table 2). The success of the models can be judged by how well all the data from simulations collapses onto a Gaussian in the x-axis argument, indicated by the black lines. From Fig. 7 it appears that both models work well in the inertial range kξ * 1. This demonstrates that the sweeping model, with sweeping velocity v/ √ 3, characterizes very well the decorrelation of the system at high wavenumbers, and that we may take C 2 v = 1 in Eq. (3.35). This also indicates that models using the scale-dependent Lagrangian eddy turn-over time in the inertial range [6,33] give a poor fit in the inertial range.
On the other hand, the decorrelation model based on Eq. (6.2) degrades somewhat for 0.125 ≤ kξ * ≤ 1. Our decorrelation model is a significant improvement on these scales, which we ascribe partly to the fact that the average decorrelation velocity takes into account the slowing of the decorrelation with the decay of the kinetic energy.
Because of limitations in the dynamical range of the simulations, it is not possible to analyse very small values of kξ * . However, we will see that the GW power spectrum decreases at small kξ * , and so it is less important to model accurately the UETC at wavenumbers well below the peak.
Note that according to recent work in Ref. [71], the Gaussian form of the random sweeping hypothesis should change to an exponential form at large time differences. Again, this effect should be most influential at low k, below the peak of the power spectrum. We leave further study of UETC models for future work.

Results: the gravitational wave spectrum
This section is organized as follows. In Section 7.1 we study a stationary source, meaning that we assume that P v (p, τ, ζ) only depends on the time difference τ − ζ. This is not the correct model for freely decaying turbulence, but we consider it nonetheless since it is useful to illustrate some properties of the GW power spectrum of a freely decaying source.
We present the case of freely decaying turbulence in the subsequent Sections 7.2 and 7.3. Section 7.2 contains the main results of our work, and is dedicated to the case of instantaneously generated turbulence, meaning that we insert as initial conditions a fully developed turbulent spectrum. Finally, in Section 7.3, we go beyond the instantaneous generation scenario, and include a growth phase for the turbulence kinetic energy.
Before continuing, we first set the basis for the calculation of the stochastic GW background that will follow. Combining Eqs. (2.15) and (2.24), one obtains (recall that q = k − p) From this expression, it appears that the computation of the GW background involves performing a five-dimensional integration for each mode k. Integrating over the azimuthal angle of p immediately gives a factor of 2π and one is left with a four-dimensional integration. Technical details on how to tackle this four-dimensional integration are presented in Appendix A. In particular, we advocate the change of variables p = k/2 + h (and hence q = k/2 − h): this makes the symmetry around the zero of the cosine of the declination angle α =ĥ ·k explicit. After the complete decay of turbulence, occurring at a time that we denote τ fin , the GW power spectrum in the radiation era evolves as This can be obtained by matching the sourced and free solutions of Eq. (2.2) at time τ fin . Note that the function G of Eq. (2.14) depends explicitly on τ , meaning that the GW spectrum in principle still depends on time after the source has decayed. We are interested in the GW spectrum at a time τ in the radiation era long after the source has decayed, and such that all relevant wavenumbers are inside the horizon, kτ > 1. This way we can neglect the second and third terms in Eq. (2.14), proportional to 1/(kτ ) and 1/(kτ ) 2 , since they become subdominant with respect to the first one. Furthermore, we can expand the first term of Eq. (2.14) to obtain The residual time dependence consists of irrelevant, rapid oscillations in time. We can therefore perform a time average of the free SGWB solution over a long time interval around τ , satisfying ∆τ 1/k for every interesting k, to obtain This is the effectively time-independent GW power spectrum that we aim to calculate. Recall that Ω gw is normalized to the radiation energy density, see Eq. (2.15). Therefore, long after the source has decayed, the averaged energy density of the GW signal would be expected to decrease as a(τ ) −4 .

The stationary assumption
In this section, we simplify the evaluation of the GW spectrum integral of Eq. (7.4) by assuming that the turbulence is stationary. Under this assumption, the rms velocity v and the integral scale ξ remain constant in time. We further set the decorrelation velocity to v dc ≡ v * / √ 3 on all scales. In other words, we neglect the k-dependence in Eq. (3.38). Motivated by what we see in the simulations, we also set C 2 v = 1 in Eq. (3.35). Consequently, the only time dependence left in the UETC in Eq. (4.7) is the time difference in the exponential. In order to evaluate Eq. (7.4) under the hypothesis that the source is stationary, it therefore becomes natural to change the time integration variables to τ mid = (τ + ζ)/2 and τ diff = τ − ζ, where the integration over the latter is symmetric around zero. We give the form that Eq. (7.4) takes under this stationary assumption as Eq. (B.1).
In order to proceed analytically, we rely on the exponential decorrelation, and extend the integration limits of τ diff to the interval [0, ∞[. Furthermore, one can also assume that the time difference τ diff is small compared to the absolute time τ mid . Under these two further assumptions, which are not strictly speaking a direct consequence of stationarity, one can integrate Eq. (B.1) analytically.
A stationary source in principle lasts forever. However, this is not very realistic in the early Universe setting. We will therefore fix the source duration to a finite value that we express in terms of N cut , the number of eddy turnover times for which the turbulence lasts. In other words, we integrate from τ * until τ * + N cut τ ξ * . The result then depends explicitly on N cut .
In Appendix B we provide the details of the integration procedure. The asymptotic behaviour of the resulting SGWB spectrum, and its scaling with the turbulence parameters are given by (note that K = Aξk, and the value of the numerical coefficients c i are given in Appendix B). Eqs. (7.5) and (7.6) show that the scaling of the GW signal with the parameter H * ξ * depends on the source duration. If the sourcing lasts less than one Hubble time, meaning N cut τ ξ * = N cut ξ * /v * < H −1 * , one can neglect the second term in the denominator of Eqs. (7.5) and (7.6), and the SGWB scaling becomes quadratic in H * ξ * . On the other hand, if the source lasts longer than one Hubble time, the scaling is linear in H * ξ * (see for example Refs. [8,24,25,32,36]). If we had set the duration of the source to infinity, as required in principle by the stationary assumption, we would have consistently found a linear dependence on H * ξ * .
Gravitational wave generation from sound waves also presents the same duration-dependent scaling with H * ξ * . This was previously derived within the sound shell model [24,27,86,87]. In the sound shell model one can also assume stationarity and set N cut ∼ 1 as in Ref. [87]. Note that, in the acoustic case, the decorrelation time depends on the sound speed, instead of on the rms fluid velocity. The typical lifetime of the source is connected to the kinetic energy. For the turbulent case, low velocities correspond to longer eddy turnover times τ ξ * = ξ * /v * , and therefore sources that last longer. The kinetic energy is in turn expected to be connected to the PT strength (see for example Ref. [88]). The latter then ultimately influences the scaling of the SGWB amplitude [8].
The turbulent kinetic energy also naturally enters the SGWB amplitude, through v * . In particular, the stationary assumption leads to a stronger scaling with v * than both the case of instantaneous generation (see Section 7.2) and the SGWB signal arising from sound waves [27,86]: at large scales and/or for low initial turbulent velocity, Eqs. (7.5) and (7.6) show that the SGWB amplitude from turbulence within the stationary assumption scales as v 6 * , as opposed to v 4 * for both sound waves and an instantaneous turbulence generation (as we shall see in Section 7.2).
We remark that the value of v * also influences the high-frequency slope of the SGWB, see Eq. (7.6). This is a consequence of the exponential decorrelation in the velocity power spectrum Eq. (4.7). For large v * > 0.5, the second term in Eq. (7.6) dominates: consequently, the SGWB spectrum becomes shallower, and its amplitude scales as v 4/3 * . The slope of the SGWB at low wavenumber, on the other hand, is due to the steep (causal) increase of the velocity power spectrum Eq. (3.13), causing the convolution in Eq. (2.24) to be flat for k < 1/ξ * (uncorrelated in space for scales larger than ξ * ) [32]. The complete spectral shape can be obtained by integrating numerically Eq. (B.4): it smoothly interpolates between the slopes given by Eqs. (7.5) and (7.6), peaking at the wavenumber k 2v * /(Aξ * ). Examples of spectra are given in Appendix B.

Instantaneous turbulence generation
In this section, we consider turbulence generated instantaneously, and freely decaying afterwards. We consider a situation where gravitational wave production begins when the fluid has the velocity spectrum of fully developed turbulence, at time τ * . We use three different methods to compute the SGWB spectrum. First, we measure the SGWB in the simulations. Second, we numerically integrate Eq. (7.4) inserting the model for freely-decaying turbulence derived in the first part of this paper, summarized in Eq. (4.7) (see also section Section 6). Third, we assume that the GW source is constant in time, and analytically integrate Eq. (7.3) under this assumption 7 . We show that there is good agreement between the three methods for the resulting SGWB. In particular, the analytical approximation provides an easy-to-use formula, Eq. (7.16), in which the scaling with the turbulence parameters (initial rms velocity, integral scale, duration) is apparent.

Simulations
In Fig. 8 we present GW power spectra from three simulations with v * ≈ 0.1, corresponding to simulations (A), (D) and (E) in Table 2. Equivalent plots for the other simulations are given in Appendix E. We explore different regions of the power spectrum with each simulation. Simulation (D) has a large value of L dx /ξ * and is therefore able to probe the GW power spectrum around the peak and towards smaller k, whereas simulation (E) is better able to probe the inertial range due to the smaller value of dx/ξ * . Simulation (A) has more lattice sites, allowing for a larger dynamic range. For this simulation we chose L dx /ξ * and dx/ξ * such that they lie in between the values of simulations (D) and (E). Simulation (A) then probes the peak of the spectrum and intermediate wavenumbers, confirming the power laws at low and high wavenumbers in these regions.
The SGWB power spectrum builds up very rapidly -this is further discussed in Section 7.2.3; its shape reflects the initial power spectrum, P U of Eq. (5.6), representing fully developed turbulence.
By averaging over the gravitational wave power spectrum in the latter half of each simulation, we are able to smooth out time-dependent oscillations in the GW power spectrum. This averaged power spectrum is shown by the black dashed line in Fig. 8. A clear k −8/3 power law is found for the inertial range in simulations (A) and (E), and an approximate k 1 power law at small wavenumbers is seen in (A) and (D). Taking a time average over the second half of the simulation is justified since the simulation wavenumbers satisfy k > 2/(τ end − τ * ) for all simulations. This can be verified by comparing the k-range in Fig. 8 (figures for the other simulations are in Appendix E) with the values given in Table 2 (see also Table 1).
When we compute the power spectrum from our simulations, we bin modes which have a wavenumber within a given shell of k, as given by Eq. (5.16). At late times, the frequency of the oscillations of each mode increases and the modes in each bin begin to oscillate out of phase. Consequently, our simulations average over these oscillations and do not capture them, even when we do not perform an average over time. This process can clearly be seen for high wavenumbers in Fig. 8.
To show the variation of the GW power spectrum for different v * and for different choices of ξ * /dx, we plot the averaged GW power spectrum from all simulations in Fig. 9. As can be seen, simulations with different choices of ξ * /dx explore different parts of the spectrum, with small ξ * /dx mapping out low wavenumbers and the peak of the spectrum, and larger ξ * /dx exploring the inertial range. Indeed, for v * 0.1 we find good agreement between simulations (A) and (D) from low wavenumbers up to the vicinity of the peak, and with (E) at high wavenumbers. This gives us confidence that the results obtained for other v * are also valid.   . The sampling of relatively few k-space modes into discrete bins leads to some noise in the very lowest wavenumbers. It is also possible that, for these lowest wavenumbers, finite volume effects lead to a flattening of the power spectrum as this is not seen in our results using numerical integration. We have cut off the spectrum at high wavenumbers as it progressively gets polluted due to numerical precision errors in projectinguij(k) toḣij(k).
In subsequent sections we will show how these simulations can be compared with the results of the four dimensional numerical integration, and with the analytical result obtained under the assumption of a constant source in time. We plot all three approaches in Fig. 13 and find good agreement.

Gravitational wave spectrum via numerical integration of the anisotropic stress UETC
We have developed a numerical method allowing exact integration of the time-averaged SGWB spectrum, Eq. (7.4). This includes integration of the full angular dependence as well, improving on previous attempts [36,38]. We insert Eq.  We also need to insert the time evolution laws for v 2 (τ ) and ξ(τ ). In accordance with the initial conditions adopted in the simulations (see Section 5.2), we assume instantaneous generation of turbulence. Specifically, prior to τ * the velocity field is zero everywhere, then at time τ * , turbulence appears fully developed with initial rms velocity v * and initial correlation scale ξ * . Turbulence then starts decaying on a timescale of order τ ξ * , following Eqs. (3.21) and (3.22). To summarize, for the instantaneous generation scenario: In the following we set N e = 5, the value that best fits the simulations, see Fig. 5. For the decay exponents p, q, we tested both sets of values that one infers from Eqs.  [89,90], an iterative and adaptive Monte Carlo scheme. More details on the implementation are given in Appendix A.
As discussed at the end of Section 3.2, turbulence is fully dissipated after about 1000τ ξ * (the inertial range has disappeared), for reasonable choices of the initial parameters T * , v * , ξ * H * . However, there is no need to time-integrate Eq. (7.4) for so long: the bulk of the GW power spectrum is sourced on a much shorter timescale. This is due to two reasons. First, the generation of GWs is very localized in time: GWs with wavenumber k are produced by the source on a timescale O(1/k), since the Green's function in Eq. (7.1), given in Eq. (2.14), has period 2π/k. Second, the decay of turbulence occurs over a timescale of the order of the initial eddy turnover time τ ξ * , which effectively cuts off the source (more precisely, about five τ ξ * , since N e 5 is the value that best fits the simulations -see Section 3.2). All GW wavenumbers with k 1/(N e τ ξ * ) are therefore sourced before the turbulent kinetic energy has decayed appreciably. This includes the peak region of the SGWB spectrum. On the other hand, GW wavenumbers with k 1/(N e τ ξ * ) are excited when the decay of the turbulent kinetic energy has already occurred: GW production on these wavenumbers is therefore negligible. This is the large scale spectral region, proportional to k 3 . We therefore expect the GW power spectrum to stop growing after a timescale of the order of 1/(N e τ ξ * ). We have numerically checked that the integration in Eq. (7.1) indeed converges to a fixed result after a few tens of τ ξ * .
In the code, we actually integrate for 100 τ ξ * : therefore, in practice we output the SGWB spectra long after the source has decayed, when the spectral amplitude does not grow in time any longer (the integral in Eq. (7.4) has converged). Furthermore, the k-range of interest is such that k > 1/(100 τ ξ * ) for every chosen value of v * and H * ξ * , as can be appreciated from Fig. 10. Therefore, we are allowed to integrate directly the time-averaged GW spectrum Eq. (7.4). We have indeed checked that numerically integrating the full expression Eq. (7.3) provides a GW spectrum which is oscillating around the amplitude of the integral of the time-averaged GW spectrum, Eq. (7.4), shown in Fig. 10. We do not plot these spectra, as they do not convey any additional information. Fig. 10 shows the GW power spectrum computed from the numerical integration code, for different values of the parameters v * and H * ξ * and for both β = 3 and β = 4. It can be appreciated that the peak amplitude scales like (H * ξ * ) 2 . The spectral peak is independent of the value of v * , and corresponds to the peak of the velocity power spectrum Eq. (3.13), occurring at kξ * 2.7/A. This is consistent with the finding of Ref. [77] for a source which is discontinuous in time.
The spectra which have H * ξ * = 10 −3 in Fig. 10 are minimally affected by the expansion of the Universe: the chosen values of v * and H * ξ * satisfy 100 τ ξ * ≤ H −1 * , and so the dynamics take place fully at sub-horizon scales. It is therefore appropriate to compare these spectra to the time-averaged spectra resulting from the simulations, in the region k > H * (see Eq. (5.18) and discussion thereafter). We will do so in sub-Section 7.2.4. On the other hand, for larger values of H * ξ * in Fig. 10, we integrate for longer than one Hubble time.
We find that the actual values of the exponents p and q do not substantially impact the SGWB spectrum, indicating that the precise details of the turbulent decay are unimportant in this respect. This is due to the nature of the free decay of turbulence in combination with the typically short timescale associated with GW production. Fig. 10 shows the SGWB spectrum from the numerical integration for both β = 3 and β = 4. For all values of v * and H * ξ * the spectra only differ in the region k < 1/(N e τ ξ * ), as expected for reasons we have outlined above. For the same reasons, we can safely ignore the role of the viscous scale λ in all our evaluations.

Constant source approximation
In this section, we show that the GW source can be approximated as almost constant in time while the majority of the GW signal is being generated. This has already been pointed out in Ref. [46], in the case of MHD turbulence. Since the Green's function in Eq. (7.1) has period 2π/k, the GW production for each wavenumber k occurs on a timescale of the order of 1/k. Therefore, over the typical decay time for the turbulence source, namely N e τ ξ * -as we assumed for Eq. (7.7) -there is sourcing of gravitational waves for wavenumbers k > 1/(N e τ ξ * ) ≡ v * /(N e ξ * ). This constitutes a large fraction of the GW spectrum, including the region around the peak k ξ * 2.7/A. Consequently, the SGWB signal on the scales of interest -specifically, around the peak -is fully established before the turbulence has decayed appreciably. On scales larger than the peak, on the other hand, the GW production is less efficient due to the decay of the source. On even larger scales, one expects the k 3 increase typical of uncorrelated sources in both time and space [77].
We assume in what follows that the turbulent GW source is instantaneously generated, remains constant for a duration N cut τ ξ * , and sharply turns off at a time τ * +N cut τ ξ * . The number N cut is of the same order as N e , and we set N cut to 7 as it gives the best fit to the result of the numerical integration. We neglect both the overall time dependence of the turbulent fluid velocity power spectrum and its time decorrelation in Eq. (4.7) This means that the anisotropic stress spectral density Eq. (2.24) becomes time-independent, and the time and momentum integrals in Eq. (7.1) decouple. We can then rewrite Eq. (7.2) as dΩ gw d ln k τ * +Ncutτ ξ * = 8 3π 2 k 3 T gw (k, τ * , N cut τ ξ * ) PΠ(k, τ * ) . (7.8) Let us first analyse the part due to the double time integral T gw . We evaluate it at the source turn-off: Restricting our consideration of this integral to wavenumbers satisfying k > 1/(τ * + N cut τ ξ * ), we are again entitled to drop the second and third term from the Green's function Eq. (2.14), to obtain [46]: The upper limits of both time integrals must be set to τ fin = τ * + N cut τ ξ * in order to match the results of both the simulations and the numerical integration (see Section 7.2.4). As pointed out in Ref. [46], the abrupt switching off of the source means that the evolution after τ fin would give an enhancement of the SGWB power at high frequency, which would break this agreement. Since we evaluate the time integrals at τ fin , we cannot perform a time average; this is why we do not use Eq. (7.4) in this Section. In Ref. [46], a simplified expression for Eq. (7.10) has been found, tuned to interpolate the envelope of the k-oscillations. In this work, we want to compare with the results of the simulations (Section 7.2.1) and of the numerical integration (Section 7.2.2), which both evaluate the time-averaged SGWB. We therefore adapt the result of Ref. [46] to account for time averaging, as follows: T gw (k, τ * , N cut τ ξ * ) ≈ ln 2 (1 + N cut τ ξ * /τ * ) for √ 2 k < 1/(N cut τ ξ * ) , ln 2 1 + ( √ 2kτ * ) −1 for √ 2 k > 1/(N cut τ ξ * ) .
The transition occurs at the wavenumber corresponding to the source lifetime, namely N cut τ ξ * (see Ref. [77]).
We now turn to the description of the anisotropic stress spectral density at initial time PΠ(k, τ * ), also entering Eq. (7.8). This is given by the convolution in momentum of Eq. (2.24) (without residual time dependence). It takes the form (we use the notation of Appendix B): where the factors of 3.0 and 3.5 come from an analytical fit, tuned to match the exact numerical integration of Eq. (7.13), and PΠ(0) is the value of the anisotropic stress spectral density at k = 0: In Fig. 11, we compare the anisotropic stress spectral density from simulation (A) with both the numerical result Eq. (7.13), and the analytical approximation Eq. (7.14). As can be appreciated from the figure, at the initial time τ * there is excellent agreement between the simulation result and the anisotropic stress spectral density derived in the context of our turbulent model. For illustration purposes, in Fig. 11 we also display the time evolution of the anisotropic stress spectral density from simulation (A): it is decaying, due in turn to the decay of the turbulent source itself. The spectral shape of the SGWB from a constant source, Eq. (7.8), is given by the combination of the k-dependence coming from both T gw (k, τ * , N cut τ ξ * ) and PΠ(k, τ * ). Recalling that the values of A and B can be found in Eqs. (3.14) and (3.15), we find dΩ gw d ln k = 7 √ π Γ(13/6) The k-dependence of T gw (k, τ * , N cut τ ξ * ) has been described in Ref. [77] in the context of coherent sources, meaning those for which the UETC can be factorized. Eq. (7.10) is effectively a time-domain Fourier cosine transform, therefore the decay properties of T gw (k, τ * , N cut τ ξ * ) at large k depend on its time differentiability -and in the case being studied here, it is discontinuous in time. This gets combined with the behaviour of PΠ(k, τ * ), which is uncorrelated for wavenumbers k < 1/(Aξ * ), while at high k it shares the same power-law decrease as P v (k, τ * ) (hence, for an initially von Kármán spectrum, it decreases as k −11/3 at high k) [20]. Altogether, for instantaneous turbulence generation of a von Kármán spectrum, the expected slopes for low v * are (with notation K = Aξ * k): This behaviour is in accordance with that found in Refs. [46,77]. Note that the PΠ/k 2 slope observed in Ref. [91] at high k can be interpreted in light of what has been discussed above as well as the results of Ref. [77]. In Fig. 12, we plot the constant source approximation, Eq. (7.16), and compare it to the result of the full numerical integration presented in Section 5. As discussed above, we set N cut = 7 to match the SGWB spectra obtained from the numerical integration. The spectral shapes of the SGWBs evaluated under the constant source approximation are in good agreement with those produced with the numerical integration, over the tested range of values for the input parameters v * and H * ξ * . The discrepancy around the transition between regimes at k = 1/( √ 2N cut τ ξ * ), see Eq. (7.16), is most notable for sources that last longer than one Hubble time (meaning N cut τ ξ * H * > 1), where the intermediate power law for wavenumbers smaller than the peak is less pronounced. Note that in Ref. [46] it was possible to adjust the source duration (i.e. the value of N cut ) in the constant approximation, as a function of v * and H * ξ * . This was done by comparison with the SGWBs output by the simulations, which covered also the large scale part of the spectrum. Here, on the other hand, we fix N cut = 7 for all initial v * and H * ξ * .
The UETC of decaying turbulence Eq. (4.7) has both a factorizable (coherent) component and a decorrelating one. The coherent component shapes the SGWB spectrum at large scales and around the peak. On the other hand, the decorrelation influences the spectral shape at small scales, as discussed in Section 7.1. The discrepancy that can be observed in Fig. 12 at high wavenumber for large v * can therefore be interpreted in light of what found under the assumption of a purely stationary source, see the asymptotic behaviour at high k of Eq. (7.6). The constant source approximation, in which the decorrelation is absent, predicts a steeper power law -given in Eq. (7.17) -than that found by the numerical integration, which contains the full time evolution of the source.
All in all, it emerges that both the overall decay of turbulence as well as its decorrelation play a relatively minor role in shaping the SGWB signal: the first influences the spectral shape only at large scales kξ * < v * /N e , while the second does so only for large initial velocity v * and at small scales kξ * 1.

Comparison of the SGWB spectra
In Sections 7.2.1 to 7.2.3, we have presented three different methods of estimating the GW power spectrum when turbulence is generated instantaneously. We first presented the GW spectra from simulations, then we performed an exact numerical integration and finally computed an analytical approximation assuming a constant source.
In this section, we compare the GW spectra obtained from the three different methods. The comparison is performed (i) long after the source has dissipated and (ii) for the lowest value of H * ξ * for which we performed the numerical integration, namely H * ξ * = 10 −3 . Condition (i) is necessary to perform the time average, while condition (ii) is necessary since the simulations neglect the effect of expansion. For more detail, we refer the reader to the discussions at the end of Sections 5.5 and 7.2.2.
We plot the late-time averaged GW power spectra from our simulations in Fig. 13, and compare them to those computed in the full numerical integration and in the constant source approximation. We show the GW power spectra from the numerical integration and the constant source approximation  Table 2, which are plotted using colored lines as specified in the legend. The GW power spectra from simulations has been cut off at high wavenumbers due to numerical precision noise.
for v * = 0.03, 0.1 and 0.3. There is excellent agreement between the SGWBs obtained with the three methods, for all the values of v * we simulated. They are all consistent with respect to the GW power spectrum amplitude, the breadth of the peak, the peak location at kξ * 2.7/A and the power laws at low and high wavenumbers. In particular, they all agree in the regions we expect the simulations to be accurate, namely at low wavenumbers for simulations (B), (D) and (F), and at high wavenumbers for simulations (C), (E) and (G). Simulation (A), which is our largest simulation, has close agreement with the other methods for the full range of wavenumbers depicted. We have therefore validated that the numerical integration and the constant approximation are accurate in the limit of H * ξ * 1. Furthermore, as they include the effect of expansion, we expect them to also be accurate outside this limit. As discussed in the previous section and depicted in Fig. 12, there is a discrepancy between the constant source approximation and the full numerical integration for wavenumbers close to the kink at k = 1/( √ 2N cut τ ξ * ). This can also be seen in Fig. 13. As the constant source approximation displays excellent agreement around the peak and obtains the correct asymptotic power laws, Eq. (7.16) can be used as a simple approximation to the GW power spectrum which does not require simulations.

Phase of turbulence growth: results from the numerical integration
In this section we go beyond the scenario of instantaneous turbulence generation, considered so far, and include a growth phase for the turbulence kinetic energy. Since the Reynolds number in the early Universe is very large (of the order of 10 13 at the electroweak scale, see Section 7.2.2), turbulence is expected to arise from vorticity generated during the PT [16], or due to the interaction of shocks Figure 14. Diagram for the evolution of turbulence in terms of η and ζ. The injection of kinetic energy starts at τini and turbulence develops in region 1 on a timescale τgro, the value of which depends on the model: it might correspond to the PT duration β −1 , or to the eddy turnover time τ ξ * . In region 2, the turbulence is freely decaying. The growth and free decay phases are correlated: regions 3 and 4 also contribute to the production of GWs. [29,30]. In order to model the GW source properly, one would have to model the onset of turbulence by simulating the complete system of scalar field and fluid. We leave this complicated problem for future work.
In the present analysis we model the (purely vortical) turbulence growth phase heuristically. This allows us to gauge in a simple way the importance of this phase, as far as the SGWB spectral shape is concerned. Refs. [36,37,77] have demonstrated analytically that time continuity of the GW sourcing process must be ensured, since it plays a role in shaping the SGWB signal, and Refs. [42,92,93] also obtained different SGWB spectral shapes in the forced case. As we shall see, we confirm this result in the present work, via the numerical integration of the SGWB source (see Section 7.2.2).
Starting from the initial time τ ini , we assume that turbulence is sourced on a timescale τ gro . At time τ * = τ ini + τ gro , turbulence is fully developed, and then starts decaying (see Fig. 14). The duration of the growth phase would depend on the particular mechanism sourcing the turbulence. If the sourcing process is related to the PT dynamics, one would expect the duration of the growth phase to be of the order of the duration of the PT. Alternatively, one can assume that it takes approximately one initial eddy turnover time τ ξ * = ξ * /v * to build the turbulent cascade. This is the assumption we make in the following, namely that τ gro = τ ξ * .
Concerning the correlations of the velocity field, the Gaussian decorrelation scenario of Sections 3.3.1 and 4.2 holds in principle only in the free-decay phase. However, in the absence of a better model, we assume that Eq. (4.7) is valid also during the growth phase, and in the mixed regions. Therefore, in the UETC of Eq. (4.7), the times τ and ζ can pertain either both to the growth phase, or both to the free decay phase, or to the two phases mixed (see Fig. 14).
We consider two scenarios for the generation of turbulence: a C 0 growth phase and a C 1 growth phase. First we adopt the model of Ref. [37] and suppose that the vortical kinetic energy grows linearly in time before turbulence enters the free-decay phase. Furthermore, we assume that the integral scale remains constant during the growth phase, and starts growing during the phase of free decay in the In the second scenario the kinetic energy and the integral scale evolve in such a way that they are both continuous and differentiable at time τ * . The kinetic energy is also continuous and differentiable at initial time τ ini . This scenario requires defining two extra functions. For the growth phase of the kinetic energy we use a smooth step function with the properties SStep(1) = 1, SStep (1) = 0 (and also, SStep(0) = SStep (0) = 0): x > 1. (7.20) To connect with the decay phase, we use a smooth power law with the property SPL(1, p) = 1 and SPL (1, p) = 0: where p < 1. Using the above functions, we set for the C 1 scenario: Different choices for the growth phase lead to different spectral shapes in the resulting SGWB. This is demonstrated by Fig. 15, showing the SGWB with both the C 0 and C 1 growth phases, as well as for the instantaneous generation scenario. While the spectra at large scales are comparable, inserting a continuous growth phase shifts the spectral peak position from the characteristic length scale of the source, k peak ξ −1 * , to the characteristic timescale of the source, k peak τ −1 ξ * : the peak position becomes therefore velocity-dependent [77]. This behaviour is caused by the component in the velocity field UETC Eq. (4.7) which is factorizable in time, termed the "coherent" case in Ref. [77]; see the discussion below Eq. (7.16). Furthermore, the addition of a growth phase which is continuous in time leads to a reduction in the GW power, and a steeper slope at high frequencies, as shown in Fig. 16 (in this figure we only plot the results for the C 1 growth phase). Since the peak now corresponds to the inverse eddy turnover time, the k 1 region observed in the instantaneous growth scenario never develops, as illustrated by Fig. 15. Note that the spectra in the C 0 and C 1 growth phase are almost identical, contrary to what was observed in Ref. [77]. Indeed, we would not expect Ref. [77] to capture all the complexity of the spectral changes due to the insertion of a continuous growth phase, since it studied a very simplified case in which the anisotropic stress UETC was separable in time.
At first, it may seem counter-intuitive that the addition of a growth phase decreases the energy in GWs, given Mercer's condition (see Section 4.1). We illustrate this using Fig. 14: Mercer's condition Eq. (4.2) applies to intervals of the form I × I, hence the contributions from the regions 1 and 2 are positive. However, the contributions from the regions 3 and 4 can very well be negative, lowering the total GW energy. This effect depends crucially on the correlations between the growth and the free decay phase. We leave the study of the growth of turbulence from more realistic initial conditions to a future work.

Discussion
In this paper, we have studied the GW signal generated by a phase of freely decaying vortical turbulence in a relativistic fluid. We were motivated by the fact that thermal first order phase transitions occurring in the early Universe are expected to lead to a turbulent flow, either directly from vorticity generation during bubble collisions, or from the interaction of shocks following the acoustic phase. Since the sourcing of the turbulent flow remains an active area of research, we have chosen to focus on GW production in the phase of turbulent free decay, inserting as our initial conditions a fully developed turbulent spectrum, as if the turbulence generation was instantaneous. While this is not the most realistic assumption, it simplified our analysis, allowing us to build a thorough understanding of the temporal and spatial structure of both the velocity field and the subsequent GW signal. We thereby derived a model for the UETC of the velocity field in freely decaying turbulence. The GW source, built in the context of this model, was then integrated to compute the GW spectrum. The validity of our model and of the integration technique was supported with a series of simulations using the Minkowski space relativistic hydrodynamics code SCOTTS. The simulations were prepared with equivalent initial conditions to our model, namely a fully developed turbulent spectrum.
The velocity of the turbulent fluid was assumed to follow the von Kármán power spectrum, interpolating between a causal k 5 slope at low wavenumbers and a Kolmogorov k −2/3 one at high wavenumbers. This basic choice influences the spectral shape of the anisotropic stress power spectrum and consequently of the GW power spectrum. Concerning the velocity UETC, motivated by the Kraichnan sweeping model, we assumed that the velocity field decorrelation is Gaussian in time, with a timescale set by the Eulerian eddy turnover time. We implemented the decorrelation at larger scales, around the spectral peak, by extending the decorrelation velocity according to Ref. [39]. We further addressed the modelling of decorrelation in the context of the theory of positive kernels, and propose to model the UETC of the turbulent velocity field as a Gibbs kernel. This both consistently enforces time symmetry in the UETC, and guarantees that the GW power spectrum is non-negative. We have validated our model for the velocity field decorrelation with a series of relativistic hydrodynamic simulations performed with SCOTTS, with the initial root-mean-squared velocity, v * , ranging from 0.03 to 0.3. We used our simulations also to study the time evolution of both the kinetic energy and the integral scale, allowing us to fix the number of eddy turnover times that it takes for the turbulence to decay.
We have reviewed the analytical computation of the SGWB from a purely vortical fluid, arriving at an equation for the GW power spectrum in terms of an integral over the Green's function and the anisotropic stress UETC, which is provided by the turbulence model we have developed. The computation of the GW spectrum could then be performed without any further approximation using Monte Carlo integration with importance sampling. Varying the exponents of the turbulence decay laws had limited impact on the resulting GW spectrum. We therefore conclude that our predictions are robust with respect to this unknown. For instantaneously generated turbulence, we found excellent agreement between the GW signal evaluated with direct numerical integration and the one extracted from our simulations.
To conclude our exploration of the GW signal in instantaneously generated turbulence, we studied the case of a turbulent source which is constant in time. As the GW spectrum at high wavenumbers and around the peak is built up on short timescales in comparison to the decay of the source, the constant source approximation works well in this region. The GW power spectrum derived within the constant source approximation is in very good agreement with both the simulation result and the power spectrum arising from direct numerical integration of the turbulence model we have built. We provide an analytical form for the SGWB signal from turbulence, Eq. (7.16), which has the advantage that it can be readily used without the need to perform simulations or the challenging computation of 4-dimensional integrals.
Finally, we explored how the GW spectrum was modified by inserting a growth phase for the turbulence, therefore removing the time discontinuity inherent in the instantaneous generation assumption. We considered in particular a C 0 growth phase and a C 1 growth phase. These models, although simplistic, allowed us to show how the GW signal might be modified by the turbulence generation phase. In contrast to the decay law exponents, the form of the initial growth phase has a large effect on the GW power spectrum: flows with equivalent v * can lead to GW signals with different spectral shapes and amplitudes, if the initial condition differs. This underlines the importance of understanding the turbulence formation mechanism in the context of first order phase transitions.
Future works will therefore focus on the generation of turbulence -which we showed is crucial to make precise predictions for future GW experiments -and on the inclusion of the longitudinal velocity field and a magnetic field. Academy of Finland grants 319066, 320123 and 308791. The work of PA (ORCID ID 0000-0002-4814-1406) was partially supported by the Wallonia-Brussels Federation Grant ARC №19/24 -103. We acknowledge PRACE for awarding us access to HAWK at GCS@HLRS, Germany.

A Integration variables for the evaluation of the GW power spectrum
This appendix summarizes the changes of variable we perform in order to tackle the convolution integral d 3 p in Eq. (7.1). To make the symmetries of the integrand apparent, we first set h = p − k/2, see Fig. 17. Denoting α =k ·ĥ, it follows that The projection factor becomes (note the symmetry α → −α) For the numerical integration, it is useful to define the new variable M ≡ k 2 /(k 2 + h 2 ). The volume element of the spherical coordinates (h, θ, φ) can then be written as with M ∈]0, 1] and α ∈ [0, 1]. Note that the integrand in Eq. (7.1) is invariant under rotations of the azimuthal angle φ.

B Analytical expressions for the stationary assumption
Under the stationary assumption, whose implications are described in Section 7.1, Eq. (7.4) can be rewritten as:   (3.13). In order to proceed analytically, we make some extra assumptions: (i) we push the upper limit of the integration in dτ diff to +∞; (ii) we neglect the term τ 2 diff /4 with respect to τ 2 mid in the denominator of Eq. (B.1). Indeed, since τ 2 diff /4 ≤ τ 2 mid /4, we expect that the Gaussian decorrelation favours the region τ 2 diff 1 for the generation of GW. Under these assumptions, the 4 dimensional integration reduces to a 1 dimensional one. First, the integral over the time difference becomes (B.2) Additionally, the integral over absolute time starting at τ * = H −1 * and lasting until τ fin = τ * + N cut τ ξ * where τ ξ * = ξ * /v * , gives Finally, it is possible to integrate analytically over α = cos θ, yielding the GW power spectrum in terms of the 1-dimensional integration (we define H = Aξh and K = Aξk) dΩ gw d ln k = 32 2 1/3 55 where the function f (H, K) is given by with f 1 ≡ 16(415H 2 + 171)K 6 + 768(H 2 + 1) 2 (5H 2 + 9)K 2 + 6912(H 2 + 1) 4 (B.6) Within the stationary assumption, the GW spectrum peaks at k 2v * /(Aξ * ), while in the instantaneous generation scenario the peak occurs at k 2.7/(Aξ * ). The stationary assumption also leads to a substantial reduction of power at small scales. In Section 7.1 we provide the low and high wavenumber asymptotic behaviour of the SGWB (B.4), for which the prefactors are where we have fit the factor Υ 8.

C Finite volume tests
We performed two additional test simulations (A ) and (A ) to understand how the low-wavenumber behaviour seen in simulation (A) depends on the cutoff induced by the box size. As shown in Table 2, these simulations have the same input power spectrum as simulation (A), but with different simulation volumes. By reducing the number of lattice points by a factor of 2 or 4 for simulation (A ) and (A ) respectively, there is less dynamic range above the integral scale ξ * . The small differences in v * and ξ * are explained partially by the use of a different random seed and partially by the different simulation volume.
In Fig. 19 we plot the final velocity power spectra from each simulation. It can be seen that even by the end of the simulation, there is broad agreement between the simulations for the wavenumbers that they share in common. We plot the evolution of ξ and v 2 in each of the simulations in Fig. 20b. From this we can see that the evolution of the kinetic energy is relatively similar, however for (A ) we see that the integral scale is unable to evolve to larger values due to the limited volume. Simulation  A ) has a closer agreement with (A) at early times, but still deviates at late times as the ξ continues to grow. This effect then modifies the rest of our global plots. In Fig. 21, we show the evolution of the instantaneous kinetic energy and integral scale exponents (p, q), both parametrically and as a function of time. We see that the simulations agree on p and q at early times, but diverge in particular for the exponent q at later times. This is unsurprising, as when q is positive it indicates that the integral scale is growing. Eventually ξ becomes comparable to the maximum scale within the simulation and the fluid flow becomes increasingly sensitive to the size of the box. At this point it appears the growth of ξ is restricted, leading to a reduction in q.

E Other GW power spectra from simulations
In this appendix we present the plots of the GW power spectrum from simulations (B)-(C) with v * ≈ 0.03 in Fig. 25, and (F)-(G) with v * ≈ 0.3 in Fig. 26.    Table 2 respectively. The coloured lines show the GW power spectrum with interval ∆τ as listed in the caption, with darker shades corresponding to later times. The black dashed line shows an average over the GW power spectrum in the last 50% of the simulation. We have cut off the spectrum at high wavenumbers due to numerical precision noise.    Table 2 respectively. The coloured lines show the GW power spectrum with interval ∆τ as listed in the caption, with darker shades corresponding to later times. The black dashed line shows an average over the GW power spectrum in the last 50% of the simulation. We have cut off the spectrum at high wavenumbers due to numerical precision noise.

F Vortical and longitudinal velocity in simulations
In this work we assume that the velocity field is purely vortical. While the longitudinal component can be set to zero in the analytical modelling or during the numerical integration, in the numerical simulations a longitudinal component can in principle develop even if it is removed in the initial conditions. In the following short appendix, we present the evolution of the kinetic energy, v 2 , within the simulations, split into its vortical and longitudinal components. Note that in the following we adjust our notation to explicitly indicate both the longitudinal and vortical components of the velocity field. We can decompose the spectral density into vortical and longitudinal components, v i (k, τ )v * j (q, τ ) ≡ (2π) 3 δ(k − q) ⊥ ij (k)P v ⊥ (k, τ ) + (1 − ⊥ ij (k))P v (k, τ ) . (F.1) From these we can construct the vortical and longitudinal kinetic energy analogous to Eq. (3.8), where the vortical and longitudinal power spectra are given by and the total kinetic energy is v 2 (τ ) = v 2 ⊥ (τ ) + v 2 (τ ). with v * 0.3 have the largest increase of v 2 at the start of the simulation, and interestingly this remains relatively constant throughout the simulation. By the end of simulation (F) the ratio of the longitudinal to total kinetic energy is v 2 /v 2 = 0.176, but by this time the total kinetic energy has decayed by more than a factor 20.