{'fraction_non_alphanumeric': 0.025182048567670742, 'fraction_numerical': 0.007987681647579636, 'mean_word_length': 4.434797768479777, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gravitational radiation from a bulk flow model

We perform simulations in a simple model that aims to mimic the hydrodynamic evolution of a relativistic fluid during a cosmological first-order phase transitions. The observable we are concerned with is hereby the spectrum of gravitational radiation produced by colliding fluid shells. We present simple parameterizations of our results as functions of the wall velocity, the duration of the phase transition and the latent heat. We also improve on previous results in the envelope approximation and compare with hydrodynamic simulations.


Introduction
It is very intriguing that gravitational wave observatories like LISA [1] might be able to probe the sub-atomic world. Nonetheless, several of these links exist: the impact of the QCD equation of state on neutron star mergers [2], gravitational waves from inflation right after the big bang [3], gravitational signals from topological defects [4] or cosmological phase transitions in the early universe [5]. The last examples are particularly interesting, sincein principle -energy scales can be probed that are even beyond the reach of today's particle colliders.
In case a cosmological phase transition is of first order, it proceeds via bubble nucleation. The latent heat is then partially converted into kinetic energy of the Higgs field as well as into a bulk flow of the plasma that fills the early Universe. The contribution to gravitational waves from the Higgs field is probably well captured by the envelope approximation [6] while the contribution from the bulk flow is much harder to quantify. The state-of-the-art calculations are hereby large scale hydrodynamic simulations [7][8][9] 1 .
Even though these hydrodynamic simulations are invaluable and shaped our understanding of GW production from phase transitions, they are not without flaws. First, these simulations are very costly and hence typically simulate a relatively small volume with a limited number of bubble nucleations (typically below 100). Often, all bubbles are nucleated at the same time in order to avoid that the first bubble dominates the simulation or even extends beyond the simulation volume. The high costs also limit the systematic study of the parametric dependence on the wall velocity and the latent heat. Another weakness is that the simulation can only cope with relatively small fluid velocities which stems from the general numerical stability of hydrodynamic codes in the non-linear regime. Finally, the dynamic range is constrained by the number of grid points on the lattice. This also prohibits the simulation of highly relativistic bubble wall velocities due to the Lorentz contraction of the Higgs bubble walls. This motivates the search for simplified models as the analytic approach in [10], the sound shell model [11] or the bulk flow model [12].
The aim of the present work is to present light-weight simulations in a bulk flow model [12] that bypass some of above short-comings. Our analysis is based on the algorithm that was previously used for simulations in the envelope approximation as presented in [13]. This approach assumes a vanishing wall thickness and allows for the study of highly relativistic bubble walls. The dynamic range is only limited by numerical accuracy and we present four octaves of frequencies for the gravitational wave spectrum. Ultimately, we want to assess to what extent this simple model can reproduce the results from hydrodynamic simulations.
The plan of the manuscript is as follows. In Section 2 we describe the model and the general formalism to calculate the GW spectrum. In Section 3 we present a more detailed account how the algorithm works that tracks bubble collisions and models the fluid. Our results are presented in condensed form in Sec. 4 and are discussed in Sec. 5.

Methodology
In this section, we present the model we employ to mimic the dynamics of the hydrodynamic fluid during (and after) the first-order phase transition. We also review the basic formalism of gravitational wave production and how different histories of nucleating bubbles are generated in our framework.

The bulk flow model
The first simulations of gravitational waves have been performed in the so-called envelope approximation [6,14] (see also [13,15,16]). The main assumption was hereby that the majority of the energy density resides in the uncollided regions. This was motivated by simulations of the scalar field -however without the relativistic fluid that filled the Universe at early times. In the recent years, it became apparent [7][8][9] that the fluid might actually be the main source of gravitational radiation. The reason being that the sound waves in the fluid persist for a much longer time. Once all initially nucleated bubbles have collided and percollation is finished, the source in the envelope approximation vanishes, while the fluid still stores most of the latent heat in form of bulk motion.
The main purpose of the present work is to devise a simple model that will mimic this effect and still allows for simplified simulations along the lines of the envelope approximation. As a guiding principle we use hereby energy conservation. Consider a set of bubbles that nucleates at the beginning of the phase transition. The bubble walls reach a terminal expansion velocity v b quickly and we assume in our model that the bubble walls are infinitely thin.
During the expansion of the uncollided bubble, the latent heat drives the expansion of the bubble. Since the volume of the broken phase scales as R 3 (R is the bubble radius) and the surface of the bubble scales with R 2 , energy conservation requires that the energy density per surface element increases linear with the bubble radius R. In contrast, once the bubble surface collided with some surrounding bubbles, no further latent heat is injected into the fluid. We assume that in this case, the bubble surface continues to propagates with velocity v b and accordingly the energy density per surface element has to decay as 1/R 2 . In this model, the first collision of a surface element determines at what point the energy density will start to fade away (in the envelope approximation, the anisotropic stress is abruptly removed at this point). Subsequent collisions with other bubble surface elements will not change the dynamics anymore. For the lack of a better name, we call this model in the following the bulk flow model.
This picture was first put forward in the work [12] and the resulting energy density is sketched in Fig. 1. The model has clear limitations that have been partially already addressed in [12] and will be part of the discussion below. The main advantage -compared to the envelope approximation -is that it captures the fact that the fluid motion sources gravitational waves long after the phase transition is finished. In this respect it is very similar to what is observed in hydrodynamic simulations. The main advantage compared to hydrodynamic simulations is that the simulations are much more economic. This allows to study a wider dynamic range (we will study the frequency range f /β ∈ [10 −2 , 10 2 ]) and also simulations with a larger number of bubbles and/or a large number of nucleation histories. We provide results for a wide range of wall velocities v b ∈ [0.01, 1].

Nucleation histories
Besides the bulk flow model, we also will improve in some aspects on the results in envelope approximation that have been reported in [13] and [16]. One main difference between the results here and the ones in [13] is that we apply periodic boundary conditions while [13] mirrored the bubble configurations at a large spherical boundary bubble. Periodic boundary conditions are also used in the hydrodynamic simulations which makes the results more compatible with this approach. On the other hand, the methods we use below will allow us to only calculate the produced GW spectrum along the three symmetry axes of the system. With a spherical boundary, the GW spectrum along arbitrary directions could be determined. In turn, this means that with periodic boundary conditions, more simulations have to be run or simulations with a larger bubble count. Figure 2 shows the distribution of the number of bubbles during the phase transition for 1024 random nucleation histories. The probability of nucleating a bubble follows the simple exponential law dP = dt exp(β t) . (1) In practice, we choose a time step dt such that dP < 1% and randomly determine if a bubble is potentially nucleated in this time step or not. The nucleation position is also randomly chosen and the bubble is injected in the simulation in case the position is still in the symmetric phase. The volume of the simulation was hereby V = (20/β) 3 and the average number of bubbles is 316.8. There is a clear correlation between the number of bubbles nucleated with the time of the first bubble that nucleates which typically happens somewhere in the range t ∈ [−2, 2]/β. The nucleation of the last bubbles is rather sharply peaked around t 7.5/β. Towards the end of the phase transition, the first bubble typically fills 15% of the total volume. So far, most simulations only chose a handful of nucleation histories. For example Ref. [13] chose eight scenarios while hydrodynamic simulations often only chose a single one for a specific parameter set and/or nucleate all bubbles at the same time. In the present work, we perform a more careful analysis and study a larger number of histories and some trends in the prediction of GW spectra. Our strategy to is to randomly chose 3 × 16 histories that come in three groups. The first group consists of histories with a number of bubbles in the range

Gravitational wave production
For the determination of the produced GW radiation we follow the approach from Refs. [6] and [13] that is based on Weinberg's master formula [17]. The energy fraction in GW radiation per octave is given as where we used the definition of the Hubble parameter and α = ρ vac /ρ rad denotes the latent heat in units of the radiation energy. The star indicates that this is the fraction in energy density at the time of production. The energy fraction today is still subject to red-shifting. The dimensionless function ∆ is defined as and the functions C ij contain the geometric information of the nucleation history and Λ projects on the transverse traceless subspace alongk. It is worth mentioning that for a large volume V (in terms of β −3 ) with N bubbles, one obtains for small wall velocities the scaling and hence we recover the well known result [18] At this stage it is advantageous to use cylindrical coordinates where the two polarizations of the GWs are denoted as + and × following the standard notation. Finally, the geometry of the nucleation history is encoded in the functions and Here, we defined the time since nucleation ∆t n = t − t n and the quantity ∆t n that keeps track of the collided regions of the surfaces of the individual bubbles. For the bulk flow model, we want to weight the surface elements according to where t c (ζ, φ) denotes the time of the first collision of the surface element. This ensures that the energy density decays as R −2 after collision as discussed in Section 2.1. For the envelope approximation, the energy vanishes after the collision and accordingly, we use using the Heaviside step function θ H . Notice that in case of the envelope approximation, the function ∆t n is piece-wise constant. One only has to keep track of the collided region and can then pull the factor ∆t 3 n out of the integrals B n,p and A n,p . In the bulk flow model, the case is more complicated since now the integral over the surface elements are weighted by the time of collision. The feasibility of our algorithm hinges on the fact that the integration over the azimuth angle φ can be performed analytically. That this is really possible is the topic of the next section.
Furthermore, notice that the energy fraction Ω GW * contains in principle a double sum over all bubbles. At first sight, it seems strange that this leads to the scaling mentioned in (5). The reason for this scaling is that the collision region of every bubble is only correlated with the surrounding bubbles. All bubbles that are rather far apart are uncorrelated and their contribution average out in the double sum. This is essential to obtain the scaling (5) and hence a volume independent result for the energy fraction. And this is exactly what the simulation yields.
Finally, let us comment on the velocity dependence in this formalism. The explicit dependence on the wall velocity v b in (8) is inherited from the physical bubble radius R n = v b ∆t n . But there is also an inherent dependence in the nucleation history. Still, it is possible to a large extent to only perform one simulation with bubble wall velocities of the speed of light and then rescale the result to other wall velocities. Since the nucleation history is invariant under a rescaling of space and time, there are two equivalent ways of obtaining the correct scaling.
First, consider a simulation with bubble wall velocities of the speed of light and then rescale the physical coordinates by a factor v b . In such a simulation, the expansion velocity of the bubble walls is v b . The nucleation history is still valid. The only difference between the two scenarios is that the volume for v b < 1 is smaller which results in an overall factor in the probability (1) and a shift in the time t n → t n − 3 log v b . This shift is unobservable in the final GW spectral density. The rescaling of the spatial coordinates reproduced the correct factor in the exponent of A n,p and also the prefactor v But we also see, that the relative phase between the bubbles should be replaced by Using this rescaling, one can produce the spectra for different wall velocities easily. In particular, the functions B n,p will not depend on the wall velocities. Moreover, A n,p depends only on the product v b ω. If wall velocities v b and frequencies ω are distributed on a logarithmic scale, only a relatively small number of evaluations of A n,p are necessary. In the following we use 60 values of v b ω spread over six orders of magnitude (two orders in v b and four orders in ω).
Compare this to 20 * 40 = 800 evaluations that would be required without this simplification. As mentioned before, there is a second way of scaling the simulation, namely scaling in time. Consider a simulation with bubble wall velocities of the speed of light and then rescale the physical time by a factor t n → t n /v b . Again, bubble walls expand with the speed v b in this rescaled simulation. However, this time the simulation would be a genuinely different simulation since it corresponds to a different bubble nucleation timescale β. Rescaling the time t and the frequency ω appropriately, we recover again the formulas (8) and the relative phase (14). In particular, the overall factor is dt 2 /β 5 → v 3 b dt 2 /β 5 . In conclusion, we only have to consider a simulation where the bubble walls expand with the speed of light. Spectra for finite wall velocity are obtained using above rescaling.

The algorithm
In this section, details of the algorithm are discussed. Readers without interest in the technical details of the numerical implementation can skip to the next section.

Geometric considerations
We want to evaluate the integrals in (10) in an efficient way. The main problem is that the collision time t c (ζ, φ) depends on the azimuth angle φ and the polar angle ζ = cos θ in a nontrivial way. In a first step, we devise a data structure that stores which is the first bubble a specific surface element collides with. An example is shows in Fig. 4. The data structure is divided into rows with a specific range of the polar angle cos θ. Every row stores several (at least one) ordered cells with the information which bubbles are relevant for collision.
In order to generate this data structure, we initialize the data structure with the two mirror images at z = z n ± L. So the total surface is split into two rows with ζ ∈ [−1, 0] and ζ ∈ [0, 1] and the rows contain only the corresponding neighbor in one cell. We then successively add more bubbles to the data structure. Now, consider two bubbles. The first one is nucleated at x µ a and is the bubble we are currently evaluating in A a,p . The second one at x µ b is one of the surrounding bubbles. Consider a volume element x µ that is on the surface of both bubbles. Being on the surface of the first bubble implies |x µ − x µ a | 2 = 0 and hence wherex is an arbitrary normalized direction. Being on the surface of the second bubble then implies Hence for a specific direction X µ = (x, 1) = (sin θ cos φ, sin θ sin φ, cos θ, 1) , The point C separates the surface of the bubble at x a into two regions according to (19) depending on which bubble the surface element collided first with. For example, points in the segment between C and B are inside both surrounding bubbles but collided with the bubble at x b before they collided with the bubble at x c . on bubble a the time of collision with bubble b can be calculated as Notice that this expression can lead to negative ∆t a which is clearly not the solution we are after. The surface regions of positive and negative ∆t a are split by the lineX µ δ b µ = 0, so ∆t a has a discontinuity from +∞ to −∞ there. Luckily, this will never be an issue. The square δ µ b δ b µ has to be positive, otherwise the bubble b would be nucleated inside the bubble a. Hence, the region that can be potentially covered by bubble b on the surface of bubble a is given by X µ δ b µ > 0. However, the region that is covered is not determined by this equation but by the requirement that the current bubble collides with the bubble a before the bubbles that are already in the data structure, let's say c. One then obtains the constraint that can be rewritten asX ζ = cos θ φ Figure 4: The plot shows an example for the data structure that stores which is the first bubble a specific surface element collides with. Different colors correspond to different neighboring bubbles.
In the following we normalize all constraints c µ such that the spatial components are normalized, c · c = 1. Hence we can parameterize the constraint as that allows for a straightforward evaluation ofX µ c µ > 0. For example, the extremal values in terms of the polar angle cos θ are given by Notice that only for |c t | ≤ 1 the constraint is non-trivial, since for c t < −1 the constraint is always fulfilled while for c t > 1 never. Every row of our data structure covers a specific range in terms of the polar angle cos θ and the row contains a set of cells. Every cell is bounded by neighboring cells. So if a row contains the set b, c, d, the boundary of the cells are given by the constraints c b,c , c c,d and c d,b . Adding a new bubble to the data structure then amounts to checking which cells (and rows) have to be split and which cells are completely covered and can be removed.
This data structure has to be constructed once for all times. This is beneficial, since it will contain only a relatively small number of surrounding bubbles which saves the time to check overlap with all bubbles for all time steps. For a fixed time, the cells in the data structure are eventually surrounded by uncollided regions. For every cell, one has to check, according to (18) /∆t a that can again be transformed into a normalized constraint (using the fact that the time component ofX is 1). Using this constraint, rows and cells have to be split again to account for uncollided regions. An example for the resulting energy density in spherical coordinates is shows in Fig. 5 ζ = cos θ ζ = cos θ ζ = cos θ φ φ

Integration over the azimuth angle φ
The heart of the algorithm is the based on the fact that the innermost integration can actually be performed analytically. In the envelope case this is trivial using The case for the bulk fluid model is a bit more involved since an additional factor ∆t 3 n is contained in the integrand of (10). According to (18) all relevant integrals can be brought to the form and For |α| < 1 the expression I s is generally valid while the expression for I c is valid in the domain φ 0 ∈ (−π, π). Even though the integrand is periodic, I c is only periodic up to a linear term.
For |α| ≥ 1 the integrands have singularities and the expressions are only valid up to the singularities in the integrand, cos(φ 0 ) = −1/α. Using the relation one can map the desired integrals to the regions that are well behaved. Furthermore, there is an numerical instability in I c near α 1 due to cancellation between the two terms in (26). As long as the argument of arctan is relatively small, one can expand the expression around α = 1 to obtain an approximation that is more stable when evaluated numerically.

Remaining integrations
The two remaining integrals in (8) are of the form However, we cannot use the standard fast Fourier transform algorithms to evaluate these integrals since ω is logarithmically distributed in our approach. We approximate these integrals with two methods. In the first one, we approximate the function f (x) by a second order polynomial and then perform the integration. This leads to accurate results when ω dx is large. However, when ω dx is small, this leads to numerical instabilities due to cancellations between terms with different inverse powers of ω dx. In this case, we approximate the full integrand exp(i ω x)f (x) by a second order polynomial and integrate. Due to the precision we are using ( 10 −16 for floating point numbers of long double type) these two methods are used for ω dx ≷ 3 × 10 −4 .
We evaluate the integrals for C p and A n,p in an adaptive way. First, we calculate the integrals for ω = 0 and further refine the integral until the local relative error is below 10 −7 and 10 −10 respectively. Once this accuracy is reached, the resulting grid of support points are used to evaluate the integral for all values of ω. This high precision is required to obtain meaningful results at the high frequency end of the spectrum.
One peculiarity of the bulk flow model is that the energy density in the surface of the expanding bubbles persists even beyond the end of the phase transition (see the right panels of Fig. 5). For these times, the function B n,p (ζ, t) becomes constant in time. The corresponding contribution to C p is proportional to In particular, the time integration is finite. This contribution scales as ω −1 for small ω and as ω −2 for large ω. However, the leading contribution for large ω in fact cancels against the contribution from the time integration during the phase transition, ∆t n < t < ∆t end n , such that C p scales as ω −3 for large ω. This serves as another important cross check for the stability of our numerical integration.

Results
In this section, we present the final results of our simulations. As described before, we simulate in total 48 nucleation histories in three batches with different bubble count. The simulations are performed in a box of size (20/β) 3 which leads to around ∼ 300 bubbles in each simulation. Using these three batches, one can study if the observed spectrum follows some trends with the bubble count. A strong dependence on the bubble count could distort the result if only nucleation histories with average bubble count are studied. We did not observe such strong correlations. For a bubble wall velocity of the speed of light, no clear correlation is observed while for small bubble wall velocities, the amplitude seems to be anti-correlated with the bubble count. However, this effect is rather weak (the amplitude is approximately inversely proportional to the bubble count). We weight the three batches according to the distribution in Fig. 2 which should take these correlations into account.
The 48 spectra are averaged and we assign an error to the spectrum that results from the variance of this averaging. The results are shown in Figs. 6 and 7. We parameterize the spectrum as   wheref andΩ GW denote the peak frequency and amplitude of the spectrum. The spectrum observed today is obtained by red-shiftinḡ where g * denotes the effective number of degrees of freedom in the plasma of the early Universe at the time of production. The peak amplitude∆ and frequencyω also have a dependence on the wall velocity that is shown in the right panels of Figs. 6 and 7. The figures also show fits to the data that are given in Tables 1 and 2.

Discussion
We presented novel results for the gravitational wave spectrum produced by a simplified model of bulk flow (see Sec. 2.1 for a description of the model). Besides, we also provide improved GW spectra for the envelope approximation that also serve as a validation for our code.
Compared to former analysis [13], the present simulations are performed with periodic boundary conditions. This has the advantage that they are directly compatible with the results from hydrodynamic simulations on a grid. However, with the present method, it is only possible to determine the GW radiation along the three symmetry axes of the box. This makes it harder to disentangle the variation in the GW signal due to the quadrupole nature of the radiation from the variation due to the nucleation history. In particular, the variance in our results most probably overestimates the true error.
The present simulation uses 48 scenarios with ∼ 300 bubbles averaged over six directions. In comparison, the simulation from Ref. [13] used eight scenarios with ∼ 100 bubbles averaged over 32 directions. In particular, the eight simulations have been chosen with a bubble count around 100 in mind which introduces a certain bias and a systematic error that was not quantified. Overall, the current simulations should be more precise and also give a more realistic estimate of the error. Numerically, our findings agree well with the former results in the envelope approximation (see Fig. 6). We observe a somewhat smaller GW signal for a wall velocity close to the speed of light and somewhat larger signals for small wall velocities. These effects are of order 50% which is well within the uncertainty of our method. The peak frequency is slightly (∼ 25%) smaller than previously reported.
Concerning the bulk flow model, our results also agree well with the soft tail of the spectrum (ω β) that was already presented in Ref. [12]. Unfortunately, the method used there was not able to predict the hard part of the spectrum due to oscillations in the multidimensional integrals.
Before we contrast our results with hydrodynamic simulations, we qualitatively compare the results from the bulk flow model and the envelope approximation. The main difference between the envelope approximation and the bulk flow model is the asymptotic behavior away from the peak of the spectrum. In the ultra-violet (ω β), the bulk flow model leads to a scaling ∼ k −3 , while the envelope shows only a slow decent ∼ k −1 . This is due to the fact that the envelope approximation leads to many kinks and cusps in the energy distribution and also in the time profile of the anisotropic stress. The bulk flow model, leads to a rather smooth anisotropic stress that in turn leads to a suppression of GW production for large frequencies 2 . In the infrared (ω β), the bulk flow model leads to a scaling ∼ k 1 , while the envelope shows a steeper slope, ∼ k 3 . Naively, the power spectrum has to scale at least as ∼ k 3 for causally generated spectra [19], so at first sight the scaling ∼ k 1 is surprising. Ultimately, The reason is that the source is lasting forever in the bulk flow model. At the same time, the nucleated bubbles keep expanding thus shifting the produced GW signal to smaller frequencies. In effect, this is the reason for the factor 1/ω in (30) that leads to the enhancement of the infrared part of the spectrum. For frequencies k H, the GW spectrum will scale as k 3 which is however not recovered explicitly since we neglect the expansion of the Universe. Finally, the bulk flow model leads to slightly smaller peak frequencies compared to the envelope approximation, which is due to the longer lasting source.
Some behavior we find is in agreement with hydrodynamic simulations. For example, we observe the scaling k −3 for large wave numbers that is also seen in simulations. Nevertheless, we also find qualitative differences. The main finding of the hydrodynamic simulations is that after the phase transition, a constant anisotropic stress persists in the plasma. This is also true in our bulk flow model. However, in the bulk flow model, the stress slowly moves to smaller wave numbers as just explained. In the hydrodynamic simulations the typical correlation length seems to stay constant after the phase transition. This leads to an additional enhancement of the GW power spectrum of order β/H (assuming that the Hubble expansion will ultimately damp the anisotropic stress). On the other hand, the produced GW spectrum in our model is finite even without this assumption. Besides, the most recent simulations [8,9] show that the peak of the GW spectrum in hydrodynamic simulations is due to the thickness of the sound shells and not the average bubble size. In the bulk flow model, the spectrum decreases beyond the peak that is related to the average bubble size.
How can the hydrodynamic simulations lead to a qualitatively different result than the bulk flow model? Most approximations in the bulk flow model seem easily justified, at least in the detonation regime [12], and should not have a major impact on the GW spectrum. One realistic possibility is that after the first collision, the sound shells lose energy by reducing the propagation speed. The bulk flow model explicitly assumes that the propagation speed of the sound shell is unchanged and only the energy density is reduced. In particular, the spherical shape of the nucleating bubbles is preserved while the energy density varies on the surface. Ultimately, a change in propagation speed could lead to a fragmentation of the bubbles and a constant correlation length at late times. Another possibility is that collisions of sound shells that happen after the phase transition is finished play an important role. In the bulk flow model, the energy density is a linear superposition of the energy densities of the individual bubbles. Besides, the wall thickness is shrunk to zero. This renders collisions after the phase transition irrelevant. In the hydrodynamic simulation, the velocity field is (in the weak regime) a superposition of the velocity fields of the individual bubbles. Since the anisotropic stress is quadratic in the velocity field, the regions where sound shells overlap might be of importance. The last point could also explain why the peak of the GW spectrum is related to the sound shell thickness and not the average bubble size. This line of argument is also supported by the sound shell model [11] that (assuming a Gaussian velocity field) reproduces many qualitative features of the hydrodynamic simulation.