\documentclass[a4paper,11pt]{article}
\pdfoutput=1 % if your are submitting a pdflatex (i.e. if you have
             % images in pdf, png or jpg format)

\usepackage{jcappub} % for details on the use of the package, please
                     % see the JCAP-author-manual

\usepackage[T1]{fontenc} % if needed
\graphicspath{{./fig/}}
\usepackage{bm}
\def\red{\textcolor{red}}
\def\blue{\textcolor{blue}}
\newcommand{\tildeH}{\tilde{H}}
\input macros

%\title{\boldmath Steepness of the gravitational wave spectrum from phase transitions}
%\title{\boldmath Shallow gravitational wave spectrum with sharp peak from phase transitions}
\title{\boldmath Shallow gravitational wave spectrum with a sharp peak from phase transitions}
%\title{\boldmath Shallow gravitational wave spectrum with peak from phase transitions}
%\title{\boldmath The peak on a shallow gravitational wave spectrum from phase transitions}

\author{\today}
%% %simple case: 2 authors, same institution
%% \author{A. Uthor}
%% \author{and A. Nother Author}
%% \affiliation{Institution,\\Address, Country}

% more complex case: 4 authors, 3 institutions, 2 footnotes
%% \author[1]{Jani Dahl,}
%% \author[2,3]{Ramkishor Sharma,}
%% \author[2,3,4,5]{Axel Brandenburg}
%plus others in an order to be decided upon later

% The "\note" macro will give a warning: "Ignoring empty anchor..."
% you can safely ignore it.

%% \affiliation[1]{Department of Physics and Helsinki Institute of Physics, PL 64, FI-00014 University of Helsinki, Finland}
%% \affiliation[2]{Nordita, KTH Royal Institute of Technology and Stockholm University,
%% Hannes Alfv\'ens v\"ag 12, 10691 Stockholm, Sweden}
%% \affiliation[3]{The Oskar Klein Centre, Department of Astronomy,
%% Stockholm University, 10691 Stockholm, Sweden}
%% \affiliation[4]{McWilliams Center for Cosmology \& Department of Physics,
%% Carnegie Mellon University, Pittsburgh, PA 15213, USA}
%% \affiliation[5]{School of Natural Sciences and Medicine, Ilia State University,
%% 3-5 Cholokashvili Avenue, 0194 Tbilisi, Georgia}

% e-mail addresses: one for each author, in the same order as the authors
%% \emailAdd{jani.dahl@helsinki.fi}
%% \emailAdd{ramkishor.sharma@su.se}
%% \emailAdd{brandenb@nordita.org}

\abstract{
We confirm that acoustic turbulence has a tendency to produce a
short frequency range with a steep gravitational wave (GW) spectrum,
$\OmGW(f)\propto d\rhoGW/d\ln f\propto f^9$.
However, there is a background with a shallower $\OmGW(f)\propto f$
subinertial range spectrum that was not obtained in a particular
analytic approximation used previously.
The relative amplitude of the $\OmGW(f)\propto f^9$ spectrum
over the background may therefore be small in practice.
The background is shown to result from contributions that were previously
thought negligible.
In the absence of cosmological expansion, the steep subrange grows
linearly in time, while the height of the background stays constant,
although it continues toward progressively lower frequencies.
When cosmological expansion is included, the $f^9$ spectrum reduces to
just a small enhancement of the spectral peak near the frequency where
most of the turbulent kinetic energy resides.
This enhancement may still be an observable feature of acoustic
turbulence, which is absent in vortical turbulence.
}

\begin{document}
\maketitle
\flushbottom

\section{Introduction}
\label{sec:intro}

Relic gravitational waves (GWs) from the early universe can reveal
valuable information about the underlying generation mechanism and thus
the physical conditions at that time.
A particularly interesting epoch is the electroweak (EW) era, which may
have involved a first-order phase transition.
This transition is characterized by the formation, expansion, and
subsequent merging of bubbles representing the new phase, leading to a
transition of the entire universe.
In the context of the EW phase transition, the scalar field parameter
undergoes a shift from a symmetric phase to a broken phase, resulting
in the appearance of broken-phase bubbles within the hot plasma.
These bubbles of the broken phase expand and merge with each other,
causing the universe to transit into this broken phase of the scalar
field parameter.
During this process, the energy transferred to the plasma can be a
considerable fraction of the total available energy and would therefore
be an important source of GWs.
The generation of GWs during this phase transition occurs in three
distinct stages: (i) the bubble collision phase, (ii) the sound wave
phase, and (iii) the turbulence phase.
However, the contribution from the bubble collision phase is typically
negligible compared to the other two stages.
It should be noted that in the case of a vacuum phase transition,
the bubble collision phase becomes the only source of GW production
\cite{PhysRevD.45.4514,Kosowsky:1992vn,Caprini:2009fx}.
Extensive research has been conducted to study the GW production in the
sound shell and turbulence phases through both numerical simulations
\cite{2015PhRvD..92l3009H,Caprini+16,2017PhRvD..96j3520H,RoperPol+20}
and semi-analytical or analytical models
\cite{Caprini:2009yp,Hindmarsh_2016,Hindmarsh_2019,RoperPol:2022iel,Auclair:2022jod,SB22}.
For weak and intermediate phase transitions, the contribution from
the turbulence phase is considered subdominant compared to the
sound waves phase.
This is primarily because the sound wave phase remains active for a
longer duration, thus serving as the primary source of GWs production.
Additionally, the effect of the stress from the velocity fields on the GWs production significantly dilutes before entering
the turbulent phase due to the expansion of the Universe.
Overall, understanding the various stages of the phase transition and
their relative contributions to generating GWs is crucial for gaining
insights into the nature of the first-order EW phase transition.
This paper focuses specifically on the sound wave phase and aims
to provide further insights into the analytical model proposed in
refs.~\cite{Hindmarsh_2016,Hindmarsh_2019}.

Relic GWs are characterized primarily by their normalized
energy spectra, $\Omega_{\rm GW}(f)$, where $f$ is the frequency, and
$\int\Omega_{\rm GW}\,\dd\ln f$ is the total GW energy
fraction of the critical energy density of the universe.
For GWs, from the electroweak phase transition, the relevant frequencies
lie in the mHz range, which is accessible to the Laser Interferometer
Space Antenna (LISA).
Studies of GWs from that epoch have therefore attracted considerable
attention.
Since the dispersion relation for GWs is linear, frequency spectra
translate directly into spatial spectra, where $k=2\pi f/c$ is the
wave number.
A finite graviton mass, however, would violate this correspondence, but
in practice this only affects much lower frequency ranges \cite{He+21},
which we are not interested in here.

In the study of turbulence, kinetic energy 
spectra are characterized by the wave number of the peak
$k_{\rm p}$ (also referred to as the energy-carrying wave number),
a declining inertial range spectrum for $k>k_{\rm p}$, and a rising
subinertial range spectrum for $k<k_{\rm p}$.
When GWs are sourced by an instantaneously occurring turbulent background
with a certain spectrum, the GW spectrum depends primarily on the spectrum
of the hydrodynamic stress \cite{RoperPol+20}.
In the inertial range, where the spectrum is red, i.e., dominated by
power from low wave numbers, the spectrum of the stress is also that
of the turbulence.
For sufficiently steep subinertial range spectra, this correspondence
between turbulence and stress spectra is no longer valid, unless the
spectrum is still red.
However, when the subinertial range spectrum is blue, the spectrum of
the stress is always white -- regardless how blue the turbulence spectrum
\cite{BB20}.
This can have consequences for the resulting GW energy spectrum.
Since the GW energy spectrum is proportional to the square of the first
time derivative of the strain, this corresponds to a linearly increasing
spectrum, $\Omega_{\rm GW}(f)\propto f$.
Such a spectrum is then expected for all frequencies below the peak,
but larger than the frequency corresponding to the horizon scale at the
time of generation.
Until recently, however, GW spectra were believed to have a white
noise spectrum for frequencies below the peak value, which means
that $\Omega_{\rm GW}(f)\propto f^3$.
Our considerations leading to the $\Omega_{\rm GW}(f)\propto f$ spectrum
apply when the source turns on instantaneously, so the aforementioned
white stress spectrum in the subinertial range translates to a white
spectrum for the second time derivative of the strain.
Of course, GW turbulence from the electroweak phase transition cannot
be regarded as an instantaneously established turbulence field.
Instead, the actual generation mechanism must play a role, and the
resulting velocity field may more realistically be described as a random
ensemble of sound waves.

In a recent analytical calculation by Hindmarsh and Hijazi \cite{Hindmarsh_2019}
studying the GW production during the sound wave phase, it was found that
sound (or acoustic) waves produce a very steep $f^9$ spectrum at low
frequencies for causally generated kinetic spectrum.
This result appears to be inconsistent with the findings from direct
numerical simulations of the sound wave phase of the phase transition
conducted in ref.~\cite{2017PhRvD..96j3520H,Jinno+23}.
However, it remains unclear whether this discrepancy arises due to the
inadequate resolution of the simulation in capturing the infrared part
of the spectrum accurately.
It is worth noting that the resulting GW spectra in the sound wave
phase differ significantly from those generated by turbulent processes,
regardless of whether the turbulence is vortical or acoustic in nature.
The purpose of the present work is therefore to reconsider the model
of ref.~\cite{Hindmarsh_2019} with a range of different approaches to
clarify the origin of the apparent conflict between different results.
We also discuss the possibility of identifying features of acoustic
turbulence in the GW spectra that might be observable in future.
Furthermore, we compare the GW spectra generated by acoustic and vortical
turbulence.

In \Sec{diff_approaches}, we discuss the different approaches used to
calculate the GW spectra originating from sound waves.
In \Sec{results}, we present our findings and compare the results
obtained using different approaches discussed in \Sec{diff_approaches}.
Finally, we conclude our findings in \Sec{conclusion}.
Throughout this paper, we adopt units, 
where the speed of light is unity.

\section{Our approaches}\label{diff_approaches}

\subsection{Gravitational waves spectrum}
\label{GravitationalWavesSpectrum}

Gravitational waves are represented by the transverse and traceless
components of metric perturbations.
In our analysis, we consider the background to be a homogeneous,
isotropic, and spatially flat expanding universe.
The metric describing such a universe with tensor perturbations can be
expressed as
\begin{equation}
ds^2 = a^2(\eta)\left[-d\eta^2 + (\delta_{ij} + h_{ij})dx^idx^j\right],
\end{equation}
where $a(\eta)$ represents the scale factor and $\eta$ denotes the conformal time.
The evolution of $h_{ij}$ is determined by the Einstein equation.
By employing this equation, we obtain the following linearized equation of motion for the Fourier space representation of $h_{ij}$,
\begin{equation}
\tilde{h}_{ij}'' + \frac{2a'}{a}\tilde{h}_{ij}' + k^2 \tilde{h}_{ij} = 16\pi G a^2 \Pi^{TT}_{ij}.
\end{equation}
Here, tildes symbolize quantities in Fourier space; this
convention is consistently used throughout this paper.
In the above equation, $a^2 \Pi^{TT}_{ij}$ represents the transverse
traceless part of the energy-momentum tensor $\Pi_{ij}$ of a GW source.
In terms of the stress tensor normalized by the energy density at the
initial epoch ($\rho_*$) and $\tildeH_{ij}=a \tildeh_{ij}$, the above equation
reduces to,
\begin{equation}
\tildeH''_{ij}+ \left(k^2-\frac{a''}{a}\right) \tildeH_{ij} =\frac{6 H_*^2}{a} \tilde{T}_{ij},
\end{equation}
where $\tilde{T}_{ij}=a^4 \Pi_{ij}/\rho_*$ is the normalized stress. 
In the radiation-dominated era, using $a=\eta/\eta_*$, and after
replacing $k \rightarrow k/H_*$ and $\eta \rightarrow \eta/\eta_*$
($\eta_*$ and $H_*$ denote the conformal time and Hubble parameter at
the initial epoch, respectively), the above equation reduces to
\begin{equation}\label{hijeq}
\tildeH_{ij}''+k^2\tildeH_{ij} ={\cal G}(\eta) \tilde{T}_{ij}.
\end{equation}
Here, ${\cal G}(\eta)=6/\eta$.
This symbol is used because the evolution equation for metric tensor
perturbations ($\tilde{h}_{ij}$) in the non-expanding universe case is
analogous to the equation mentioned above, but with ${\cal G}(\eta)=6$
(see ref.~\cite{RoperPol+20b} for details).
We will discuss the non-expanding case in the next section. 

Usually, one is interested in estimating the GW spectrum, $\OmGW(k)$.
It represents the GW energy density per logarithmic wave number interval normalized
by the present-day critical energy density of the Universe.
For the case when the source is active during the interval
$1<\eta<\eta_e$, where $\eta_e$ denotes the time until the turbulence
is active, and $\OmGW(k)$ is given by
\begin{align}
\OmGW(k)\Big|_0&=\frac{3k^3}{ 4\pi^2  a_0^4}\frac{H_*^2}{H_0^2}\int_1^{\eta_e} 
\int_1^{\eta_e} \frac{d\eta_1 d\eta_2}{\eta_1 \eta_2} \cos k(\eta_1-\eta_2)
|\tilde{T}|^2(k,\eta_1,\eta_2)\label{GWspec}.
\end{align}
Here $H_0$ and $a_0$ denote the Hubble parameter and the scale factor at
the present epoch, respectively and $|\tilde{T}|^2(k,\eta_1,\eta_2)$
is defined through $\langle \Tilde{T}_{ij}(\kk,\eta_1)
\tilde{T}^{ij}(\qq,\eta_2)\rangle=(2\pi)^3\delta(\kk-\qq)|\tilde{T}|^2(k,\eta_1,\eta_2)$.

\subsection{The non-expanding universe case}

When the effective duration of the GW source is shorter than the Hubble
expansion time, it is possible to neglect the expansion of the universe
in estimating the GW energy spectrum.
In such cases, $\OmGW$ is given by
\begin{equation}\label{GWstatic}
\OmGW(k)\Big|_0=\frac{3k^3}{ 4\pi^2  a_0^4}\frac{H_*^2}{H_0^2}\int_1^{\eta_e} 
\int_1^{\eta_e} d\eta_1 d\eta_2 \cos k(\eta_1-\eta_2)|\tilde{T}|^2(k,\eta_1,\eta_2).
\end{equation}
The above expression is similar to the expression obtained in
ref.~\cite{Hindmarsh_2019} by substituting their equation (3.11)
in their (3.6).
The only difference is that in our case we have normalized the stress
tensor with the energy density of the initial epoch and the GW energy
density with the present-day critical energy density.

Another reason for considering the non-expanding case is that we want to assess the validity of the approximations used in
ref.~\cite{Hindmarsh_2019}, where a non-expanding universe was assumed.
This allows us to provide a more detailed understanding of the origin
of particular features in the GW spectra.

\subsection{Calculation based on the sound shell model}
\label{Calculation}

The evaluation of the GW energy spectrum requires determining the fluid shear
stress unequal time correlator $|\tilde{T}|^2(k,\eta_1,\eta_2)$
resulting from colliding sound waves generated by an electroweak phase
transition in the early universe.
The stress tensor for the fluid is given by
\begin{equation}
\Pi_{ij}=w \gamma^2 u_i u_j,
\end{equation}
where $w = \rho + p$ is the enthalpy density consisting of the
energy density $\rho$ and the pressure of the fluid $p$.
Here, $u_i$ and $\gamma=(1-\uu^2)^{-1/2}$ are the components of the
fluid velocity and Lorentz factor, respectively.
The sound shell model of ref.~\cite{Hindmarsh_2019} assumes that the
spherically expanding phase transition fronts, often called bubbles,
continue to propagate after the phase transition is completed, each of
them acting as a driver for a sound wave.
The velocity field generated by a phase transition involving a large
amount of expanding bubbles is then at each point a collection of sound
waves resulting from many of these sound shells.
This can be treated as a Gaussian random field.
Therefore, in calculating the shear stress unequal time correlator for
a non-relativistic fluid using the standard method of expanding the
resulting connected four-point correlator using the Wick expansion,
any non-Gaussianities are assumed to be small and to give negligible
contributions.
The four-point correlator then reduces to a sum containing products of
two-point velocity correlators.
In the model, these are written in terms of the fluid power spectrum and
a stationary cosine function that is found by using the properties of
sound waves, such as their longitudinality, fulfilment of the linearized
%AB: "longitudinality" is not clear
wave equation, and the assumptions made by the sound shell model under
bubble collisions.
The result for $|\tilde{T}|^2(k,\eta_1,\eta_2)$ is equation~(3.34)
in ref.~\cite{Hindmarsh_2019}, which can be written as
\begin{equation}\label{stressinssm}
|\tilde{T}|^2(k,\eta_1,\eta_2)=\frac{16 \pi^2\bar{w}^2}{ k}
\int_0^{\infty} dq \int_{|q-k|}^{q+k} d\tilde{q}~\frac{q}{\tilde{q}^3} (1-\mu^2)^2 \EK(q) \EK(\tilde{q}) \cos(\omega(\eta_1-\eta_2))  \cos(\tilde{\omega}(\eta_1-\eta_2)).
\end{equation}
Here $\EK$ is the energy spectrum per linear wave number interval with
the normalization $\int \EK dk =\bra{\uu^2}/2$, and
$\omega=\cs q$, $\tilde{\omega}=\cs \tilde{q}$ are frequencies,
where $\cs=1/\sqrt{3}$ is the sound speed.
Furthermore, $\bar{w}=\bar{\rho}+\bar{p}$ is the mean enthalpy
density and
\begin{equation}
\label{mu}
\mu=\frac{q^2+k^2-\tilde{q}^2}{2kq}.
\end{equation}

By substituting $|\tilde{T}|^2(k,\eta_1,\eta_2)$ from \Eq{stressinssm}
in \Eq{GWstatic}, the GW spectrum can be obtained.
%AB: Jani, you wanted to add here something about the numerical method!
%AB: Mention accuracy and random sampling. (->RS)
%rs Axel, I realise that I use the adative monte carlo integration method only for the expanding case, not for non expanding case.
In their analytical study of the GW spectrum, Hindmarsh and Hijazi
\citep{Hindmarsh_2019} solved one of the time integrals by assuming a $\delta$
function; see their equation~(3.39), which is valid
in the limit of time approaching infinity.
The $\delta$ function in their approach ensures that the GW wave number,
$k$, is equal to the sum of $\omega$ and $\tilde{\omega}$.
In the following, we refer to this as the `$\delta$ function approach'.
Out of the four pairs of terms identified originally in
Ref.~\cite{Hindmarsh_2019}, only one contributes to the
steep $f^9$ spectrum.
However, as we will point out, two additional pairs of terms,
which were previously thought negligible, also contribute.
In their analysis, they found that the GW spectrum is proportional
to $k^{2n-1}$ below the spectral peak.
Here, $n$ is the spectral index of the kinetic spectrum defined as
$k\EK \propto k^n$.
However, as alluded to above, since the interval during which sound
waves contribute to GW production is finite, we numerically evaluate
the integrals with a finite time integration limit to estimate the
GW spectrum.
Our numerical analysis demonstrates that, below the peak, there is indeed a
slope of $k^{2n-1}$, although it does not extend to very low wavenumbers.
Instead, there is a transition to a linear spectrum at a point that
depends on the time integration limit.
We will provide a detailed discussion of these findings in \Sec{result1}.

%AB: we should explain why in their approach the k^9 spectrum occurs only for acoustic turbulence and not for vortical.
%rs Do not understand this fully. Let us keep this comment here to remind ourselves of this point.
%AB: Also, did we explain the expected modifications if the sub-inertial range is not proportional to k^4?
%rs I checked the steeper slope is dependent on the spectral index of the kinetic spectrum. Like for E_K \propto k^2, I got the slope of the GW spectrum propto k^5. Therefore for the general case, E_K \propto k^{n-1}, the \Omega_GW \propto k^{2n-1} as also said in H & H 2019.


\subsection{Simulations with a kinematic velocity field}
\label{KinematicVelocityField}

To compute the gravitational wave field numerically, we need to compute
the hydrodynamic stress in regular time intervals, $\delta\eta$.
In this section, we construct a three-dimensional, time-dependent, random
irrotational velocity field in Fourier space, $\tilde{\uu}(\kk,\eta)$, as
\begin{equation}
\tilde{\uu}(\kk,\eta)=\ii\kk\,\tilde{\phi}(k)\,e^{-\ii\omega(k)\eta},
\end{equation}
where $\tilde{\phi}(\kk)$ is the Fourier transformed scalar potential
of the velocity, $\omega(k)=\cs k$ is the dispersion relation for sound
waves, $k=|\kk|$ is the wave number.
We construct $\tilde{\phi}(\kk)$ such that $k\EK(k)$ has a $k^5$
subinertial range for $k<k_{\rm p}$ and a $k^{-1}$ inertial range for
$k>k_{\rm p}$.
The horizon wave number at the time of generation is $k=1$.
In our calculations, we sometimes include wave numbers below the
horizon wave number in order to see the subsequent evolution from
the $\OmGW\propto k^3$ spectrum at early times to a linearly increasing
one at later times for $k\eta>1$.

At each time step, we Fourier transform $\tilde{\uu}(\kk,\eta)$
into real space, $\uu(\xx,\eta)$, and compute $\Pi_{ij}(\xx,\eta)$.
From this, we compute the transverse tracefree (TT) projection
$\tilde{T}_{ij}^{\rm TT}(\kk,t)$ in Fourier space and express it in terms
of the plus and cross polarization modes, 
$\tilde{T}_{+/\times}(\kk,\eta)$.
We then compute the strain field $\tilde{H}(\kk,\eta)$ by solving \cite{RoperPol+20}
\begin{equation}
\left(\frac{\partial^2}{\partial\eta^2}+k^2\right)\tilde{H}_{+/\times}(\kk,\eta)=
{\cal G}\,\tilde{T}_{+/\times}(\kk,\eta),
\end{equation}
where ${\cal G}=6$ in all our calculations without expansion, and
${\cal G}=6/\eta$ in our calculations with expansion.
For the simulations, we use the {\sc Pencil Code} \cite{JOSS}, which is
a massively parallel public domain code, where the relevant equations
have already been implemented \cite{RoperPol+20b}.

Although $\tilde{\uu}(\kk,\eta)$ can be computed accurately, regardless
of the choice of $\delta\eta$, we cannot choose it to large, because
otherwise the accuracy of $\tilde{H}(\kk,\eta)$ will become poor.
The error in the solution for $\tilde{H}(\kk,\eta)$ is $O(\delta\eta^2)$,
and the coefficient in the error is smaller than the third-order accurate
solution of the hydrodynamic equations, when those are computed; see
\Sec{Direct} below.

\subsection{Direct numerical simulation with the Navier Stokes equation}
\label{Direct}

The purpose of the kinematic velocity fields discussed in
\Sec{KinematicVelocityField} is to bridge the gap between actual
turbulence, which is nonlinear, and the linear random sound waves
considered in ref.~\cite{Hindmarsh_2019}.
To model turbulence more realistically, we also solve the compressible
Navier-Stokes equations directly, i.e.,
\begin{align}
\frac{\DD\uu}{\DD\eta}&=\frac{2}{\rho}\nab\cdot\left(\rho\nu\SSSS\right)
-\frac{1}{4}\nab\ln\rho+\frac{\uu}{3}\left(\nab\cdot\uu
+\uu\cdot\nab\ln\rho\right),
\label{dudt} \\
\frac{\partial\ln\rho}{\partial\eta}
&=-\frac{4}{3}\left(\nab\cdot\uu+\uu\cdot\nab\ln\rho\right),
\end{align}
where $\DD/\DD\eta\equiv\partial/\partial\eta+\uu\cdot\nab$ is the
advective derivative, ${\sf S}_{ij}=(\partial_i u_j+\partial_j u_i)/2
-\delta_{ij}\nab\cdot\uu/3$ are the components of the rate-of-strain
tensor, $\nu$ is the viscosity, and $p=\rho\cs^2$ is the pressure.

We construct our initial condition in Fourier space by multiplying
a random vector field with a superposition of a vortical and an
irrotational contribution,
\begin{equation}
Q_{ij}=Q_0\,\left[
(1-q)(\delta_{ij}-\hat{k}_i\hat{k}_j)+q\hat{k}_i\hat{k}_j\right],
\end{equation}
where $0\leq q\leq1$ quantifies the irrotational fraction.
This tensor can also be written as
\begin{equation}
Q_{ij}=Q_0\,\left[(1-q)\delta_{ij}-(1-2q)\hat{k}_i\hat{k}_j\right].
\end{equation}
In this work, we consider the extreme cases $q=0$ and $q=1$
for vortical and irrotational flows, respectively.

The relevant control parameters are the root-mean-square (rms)
Mach and Reynolds numbers,
\begin{equation}
\Ma=\urms/\cs,\quad
\Rey=\sigma\urms/\kf,
\end{equation}
as well as the wave number $k_0$ of the spectral peak of the
initial velocity spectrum.
We recall that $k=1$ corresponds to the horizon wave number at the time
of generation.

In all our calculations without expansion, our conformal time is equal
to cosmic time.
However, as in the solutions with expansion, which start with $\eta=1$,
our initial time is still unity when there is no expansion.
%AB: when we normalize like so, those are no longer "natural units", as advertised above

\section{Results}\label{results}

\subsection{Spectral evolution with time}\label{result1}

In \Fig{figure1}, we compare $\OmGW(k,\eta)$ for $\eta=2$, 10, 50,
and 250 in the non-expanding universe case.
These spectra were obtained by numerically integrating \Eq{GWstatic}
using $|\tilde{T}|^2(k,\eta_1,\eta_2)$ given by \Eq{stressinssm}.
In our calculations, we assumed a kinetic energy spectrum of the form
$E(k) \propto (k/k_p)^4/[1+(k/k_p)^6]$ with $k_p=30$ and normalized such
that the total kinetic energy normalized by the initial background energy
density is about 0.1.
We see that the peak of $\OmGW(k,\eta)$ increases approximately linearly
with time.
It does indeed have a $k^9$ slope below the peak over a short range
$10<k<30$.
Above the peak, $\OmGW(k,\eta)$ falls off approximately like $k^{-3}$.

\begin{figure*}\begin{center}
\includegraphics[scale=0.7]{fig/spectrum_wo_exp_kp30.eps}
\includegraphics[scale=0.5]{fig/spectrum2_wo_exp_kp30.eps}
\end{center}\caption[]{GW spectrum without expansion at different times.
The amplitude of the spectrum in the range $1<k<10$ is unchanged.
However, it continues increasing at $k=30$ with $\propto \eta$.
By assuming that the GW spectra below the peak maintain the $k^8$
spectrum, the transition from $k^8$ to the flat regime moves towards
left with a speed proportional to $\eta^{1/8}$.
The transition to $k^2$ propagates $\propto \eta$ in the horizontal
direction.
}\label{figure1}\end{figure*}

%rs added the figure below
\begin{figure*}\begin{center}
\includegraphics[scale=0.6]{fig/GW_spec_peak_with_kp.eps}
\includegraphics[scale=0.6]{fig/linear_scaling_with_kp.eps}
\end{center}\caption[]{Scaled GW spectrum for different $k_p$ at time, $\eta=10$. In the left panel, the spectrum is normalized by $k_p$ and by $k_p^2$ in the right panel.
}\end{figure*}
%rs

\begin{table}[htb]\caption{
Summary of the turbulence runs discussed in the paper.
All runs have $\kp=10$ and a numerical resolution of $1024^3$ mesh points.
}\vspace{12pt}\centerline{\begin{tabular}{lccccccc}
%Run & $\sigma$ & $A$ & $k_1$ & $\nu$ & $\EEK$ & $\EEGW$ & $\hrms$ \\
%AB: sigma -> q
%rs We have used q for the wavenumber but I think it is okay.
Run & $q$ & $A$ & $k_1$ & $\nu$ & $\EEK$ & $\EEGW$ & $\hrms$ \\
\hline
%Run  sig        A           k1        nu              EEK                    EEGW             hrms          run
A1 &$1$&0.1 &$0.1$&$5\times10^{-3}$&$3.8\times10^{-2}$&$2.7\times10^{-5}$&$2.4\times10^{-3}$\\% P1024_k01_kf10e
V1 &$0$&0.1 &$0.1$&$5\times10^{-3}$&$7.5\times10^{-2}$&$1.1\times10^{-4}$&$6.9\times10^{-3}$\\% V1024_k01_kf10e
A2 &$1$&0.1 &$1.0$&$2\times10^{-3}$&$4.8\times10^{-2}$&$4.4\times10^{-5}$&$2.0\times10^{-3}$\\% P1024_k1_kf10b
V2 &$0$&0.1 &$1.0$&$2\times10^{-3}$&$9.7\times10^{-2}$&$1.3\times10^{-4}$&$6.4\times10^{-3}$\\% V1024_k1_kf10b
A3 &$1$&0.05&$1.0$&$5\times10^{-4}$&$1.2\times10^{-2}$&$7.7\times10^{-6}$&$6.5\times10^{-4}$\\% P1024_k1_kf10c
V3 &$0$&0.05&$1.0$&$5\times10^{-4}$&$2.4\times10^{-2}$&$8.4\times10^{-6}$&$2.0\times10^{-3}$\\% V1024_k1_kf10c
\label{Ttimescale}\end{tabular}}\end{table}

The reason why we see the $k^9$ subrange only for a relatively short $k$
range is the presence of an approximately stationary linearly varying
background $\OmGW(k)\propto k$ for $1<k<10$.
This background was not present in the earlier work of ref.~\cite{Hindmarsh_2019}.
As demonstrated in \App{OriginBackground}, the background emerges
because of additional contributions to the time integral that are
not negligible in practice.
As discussed in the introduction, a linearly increasing GW spectrum,
such as this background, is expected when the stress has a white-noise
spectrum in that $k$ range.
%AB: Ramkishor, have you looked at the stress spectrum?
Furthermore, as time goes on, this linear GW spectrum
continues to extend to progressively lower $k$ values as the horizon
grows, and smaller $k$ values become causally connected \cite{SB22}.

With time, the $\OmGW(k)\propto k^9$ part of the spectrum begins to
grow, and the intersection with the shallower background spectrum
moves toward lower wave numbers, but this does not happen as quickly
as the horizon wave number, which moves toward lower $k$ like $1/t$.
Therefore, the linearly growing background $\OmGW(k)\propto k$ broadens
with time.
However, since the height of the $\OmGW(k)\propto k^9$ feature continues
to grow, this steep part of the spectrum does indeed constitute an
important contribution to the spectrum when the expansion of the universe
is neglected.

To determine the ratio between the peak value of the spectrum
and the background, it is necessary to estimate the duration
during which the source remains active.
In this article, we primarily focus on the contribution from the
sound mode-dominated stage of the phase transition.
The sound mode contribution remains active until the development of
turbulence, and the relevant timescale for this is the eddy turnover
time corresponding to the peak of the spectrum.
This timescale is equal to the ratio between the peak of the kinetic
spectrum and the root-mean-square (rms) velocity of the fluid.

However, if this timescale exceeds the Hubble expansion timescale,
then considering the effect of expansion on GW production becomes
important.
The Hubble expansion time may determine the effective duration
of the sound mode contribution, and we discuss this further
in \Sec{result3}.
This scenario may be applicable to weak phase transitions and may be
of observational relevance.
For instance, when the peak of the kinetic spectrum occurs at a wavenumber
$10$ times greater than the Hubble scale and the rms fluid velocity
is 0.1, the nonlinear timescale relevant for turbulence development,
denoted as $\eta_{\rm NL}=1/(\kp\urms)$, becomes comparable to the
Hubble expansion time.
In \Sec{result3}, we provide a comparison considering the expansion for
this specific case.
%AB: in our discussions on the white board, we talked about a factor of 4.
%AB: I do now actually mention this, but maybe not yet strongly enough.

%AB: Ramkishor, we should write here something about the Adaptive Monte Carlo Method. (->RS)
%AB: In this connection, we also talked about the Vegas code; do you remember what that was?
For the semi-analytic approach, we use the Adaptive Monte Carlo Method
in Mathematica to perform the integration.
It is considerably more efficient than a straightforward implementation
of a numerical integration.

\subsection{Comparison with synthesized random sound waves}\label{result2}

We now compare with the results from our model where the stress is
constructed from synthetic sound waves; see \Fig{Comparison}.
We should emphasize here that GW energies from the semianalytical method
and the three-dimensional simulations agree with each other rather well.
We see that the kinetic energy spectra have the expected $\EK\propto k^4$
subrange for $k<k_{\rm p}$, followed by a $k^{-2}$ subrange for $k>k_{\rm p}$.
(This $k^{-2}$ subrange is expected for acoustic turbulence, especially
in the presence of shocks, although in the present case, shocks are only
present when we solve the Navier-Stokes equation.)

\begin{figure*}\begin{center}
\includegraphics[scale=1.2]{fig/spectrum_512a_2nd_dt4c.eps}
%AB: the k on the abscissa should be in italics
%AB: please also check the y-axis.
\end{center}\caption[]{Left: kinetic spectra, Right: GW spectra.
Red curves show the analytical spectra and blue curves represent numerical spectra.
In simulations, we get a broader peak. Oscillations at large $k$ are an artifact of large time steps.
}\label{Comparison}\end{figure*}

\begin{figure*}\begin{center}
\includegraphics[scale=0.6]{fig/spectrum_with_exp_kp10.eps}
\includegraphics[scale=0.6]{fig/comparison_with_and_wo_expansion.eps}
\end{center}\caption[]{
GW spectrum with expansion at different times obtained from the analytic model.
Note that the height of the peak over its background is only about one
order of magnitude.
The gray lines indicate various slopes $\propto k^4$ and $k^5$ for orientation,
but are not meant to refer to any particular theoretical prediction.
}\label{figure4}\end{figure*}

%rs added the figure below
\begin{figure*}\begin{center}
\includegraphics[scale=0.7]{fig/spectrum_with_exp_kp30.eps}
\end{center}\caption[]{GW spectrum for $k_p=30$ at times $\eta=2,4,6,8$.}\end{figure*}
%rs 


\begin{figure*}\begin{center}
\includegraphics[scale=1.25]{fig/acoustic_vs_vortical.eps}
\end{center}\caption[]{
Kinetic, stress and mean GW spectra.
The black and red curves show the kinetic and stress spectra 
at the initial and final times of simulation, respectively.
Note the presence of the bump in all irrotational runs.
It is absent in all vortical runs.
}\label{figure3}\end{figure*}

In \Fig{Comparison}, the GW spectra are shown as energy per linear wave
number interval, so as to see more clearly the transition from the white
noise $\OmGW/k\propto k^2$ spectrum for superhorizon scales ($k\eta<1$)
and the $k^0$ range for larger wave numbers, but still below $k_{\rm p}$.
For early times ($\eta=0.1$), the spectra are entirely dominated by the
superhorizon scaling $\propto k^2$, and there is neither a $k^0$ subrange
nor a $\OmGW\propto k^9$ spectrum yet.
At late times, however, a short $k^9$ subrange does emerge; see
\Fig{Comparison}(d).
In comparison with the analytic model, however,
the $k^9$ spectrum begins already for smaller $k$ ($k\approx5$ instead of 10) and
is then followed by a small plateau in the range $20<k<40$, before
the spectrum declines with a $k^{-2}$ subrange at larger $k$.

\subsection{Effect of expansion}\label{result3}
%AB: added text
As discussed in \Sec{GravitationalWavesSpectrum}, when the expansion of the
universe is taken into account, the effective duration of the turbulence
is limited and the steep part of the GW spectrum is much smaller.
This is shown in \Fig{figure4}, where we show $\OmGW(k)/k$ for $\eta=5$,
15, and 30.
%AB: Ramkishor, is the following correct? We discussed this last time. (->RS)
%rs YEs. However, this is about our full numerical simulation. I think this is not the right place for the below text.
%AB: So not the kinematic module? But how can we have so much freedom with the direct numerical simulations?
%rs No, it is about the kinetic module. However, this section is about the effect of expansion. The below text should go in the previous section.
%AB: ok, but this has not been moved yet, right?
%rs not yet.
We mention in this context that our results are insensitive to the
length of the time step and that the same results have been produced
with a time step that was up to 50 times shorter.
We see that the spectra do not change much at late times after $\eta=4$,
and that the results for $\eta=8$ and 15 are almost indistinguishable
from each other and close to those for $\eta=4$.
Furthermore, the height of the bump is limited to about one order of
magnitude.

\subsection{Results from turbulence simulations}

In \Fig{figure3}, we show spectra of $\dot{h}$, $\uu$, and $\TT$ at
different times for turbulence starting from vortical and irrotational
initial velocity fields.
In our present simulations, we have $\Ma\approx0.2$ and $\Rey\approx40...300$.
We see that $\Sp(\dot{h})$ has a bump at $k\approx10...20$ for
irrotational turbulence, but not for vortical turbulence.
This is potentially an important characteristic of acoustic turbulence.

We also see that the inclusion of low wave numbers ($k<1$ for Runs~A1 and
V1) only results in a clear representation of the $\OmGW(k)\propto k^3$
range, which is not visible for the other runs.
This aspect is independent of whether the turbulence is acoustic or vortical.
In all cases, the height of the bump is about one order of magnitude.
The height is also independent of the overall amplitude of the turbulence;
compare Runs~A2 and A3.

\section{Conclusions}\label{conclusion}

The present work has confirmed that there is indeed the very steep
$\OmGW\propto k^9$ subrange in the GW spectrum for acoustic turbulence,
as originally proposed in ref.~\cite{Hindmarsh_2019}.
However, it is usually not very prominent, because it is superimposed on
a stationary linearly rising background for all subhorizon
wave numbers above $k>1/\eta$ and below the spectral peak at $k_{\rm p}$.
The absence of this background in earlier analytical calculations is a
consequence of having adopted the $\delta$ function approach discussed
in \Sec{Calculation}.
Furthermore, the $k^9$ subrange exists only in the immediate proximity
of the peak of the GW spectrum.
It grows linearly in time, so it becomes strong only at
late times.

It must also be pointed out that the $k^9$ subrange only becomes prominent
when the expansion of the universe is neglected.
This is obviously not a reasonable approximation, especially at late
times, which is when the $k^9$ subrange has otherwise the best chance
of becoming prominent.
%AB: added this; is this correct?
In addition, the height of the bump depends on the strength of the phase
transition and occurs only for sufficiently strong phase transitions.
%AB.

Given that the $k^9$ subrange is now expected to be a relatively weak
phenomenon, we must ask ourselves whether indications of it have still
been seen in earlier work on acoustically generated GWs.
Acoustic turbulence has been considered in several recent direct numerical
simulations of turbulence \cite{RoperPol+20,BGKMRPS21}.
Contrary to the case of vortical turbulence, the present simulations
tend to show a slightly more pronounced peak.
It was not recognized as a particularly obvious feature, but now that we
know of the possibility of a subdominant $k^9$ subrange, it is tempting
to associate this slight enhancement near the peak, seen in some earlier
papers, with a specific feature of acoustic turbulence.
Obviously, more comparative studies are needed before this can be used
as a diagnostic tool in future observational studies with LISA.

\paragraph{Data availability.}
The source code used for the numerical solutions of this study, the
{\sc Pencil Code}, along with the additions included for the present
study, is freely available~\cite{JOSS}; see also ref.~\cite{DATA}
for the numerical data.

\acknowledgments
Support through the grant 2019-04234 from the Swedish Research Council
(Vetenskapsr{\aa}det) is gratefully acknowledged.
Nordita is sponsored by Nordforsk.
We acknowledge the allocation of computing resources provided by the
Swedish National Allocations Committee at the Center for Parallel
Computers at the Royal Institute of Technology in Stockholm.

\appendix
\section{Origin of the background contribution}
\label{OriginBackground}

\begin{figure*}\begin{center}
\includegraphics[scale=.9]{fig/contribution_by_parts.eps}
\end{center}\caption[]{
Contributions to $\OmGW(k) h_0^2/k$ from the four terms in \Eq{FourTerms}
for $\eta_e=5$.
Note that the steep $\OmGW(k)\propto k^9$ contribution comes from the
$S_{--}$ term, while the background results from the 
$S_{+-}$ and $S_{-+}$ terms.
}\label{contributions}\end{figure*}

As shown in Ref.~\cite{Hindmarsh_2019}, the expression for $\OmGW(k)$
can be written as
\begin{equation}\label{GWstatic_appendix}
\OmGW(k)\Big|_0=\frac{24 k^2 \bar{w}^2}{  a_0^4}\frac{H_*^2}{H_0^2}\int_0^{\eta_e} 
\int_0^{\eta_e} d\eta_1 d\eta_2
\int_0^{\infty} dq \int_{|q-k|}^{q+k} d\tilde{q}~\frac{q}{\tilde{q}^3} (1-\mu^2)^2 \EK(q) \EK(\tilde{q}) \Delta(k,\omega,\tilde{\omega}),
\end{equation}
where
\begin{equation}
\Delta(k,\omega,\tilde{\omega})=\frac{1}{2}\int_{0}^{\eta_e} d\eta_1 \int_{0}^{\eta_e}d\eta_2 \cos(k(\eta_1-\eta_2)) \cos(\omega(\eta_1-\eta_2))  \cos(\tilde{\omega}(\eta_1-\eta_2))
\end{equation}
contains the unequal time correlator derived in Ref.~\cite{Hindmarsh_2019}
for sound waves.
It has four pairs of contributions:
\begin{equation}
 \Delta(k,\omega,\tilde{\omega})= \frac{1}{2}(S_{+-}^2+S_{-+}^2+S_{--}^2+S_{++}^2)
\label{FourTerms}
\end{equation}
with
\begin{equation}
    S_{\pm\pm}=\frac{\sin ((k\pm \omega\pm\tilde{\omega})\eta_e/2)}{k\pm\omega\pm\tilde{\omega}}.
\end{equation}
%rs I think we need not to write the line below.
%AB: but Hindmarsh_2019 used 8 terms; it would have been sufficient to write 4, but they didn't. How else would you clarify this?
Each of the four terms in \Eq{FourTerms} contributes twice for a
corresponding term with opposite sign in front of $k$.
%rs

The term relevant for the steep $f^9$ spectrum results from $S_{--}^2$;
see \Fig{contributions}.
However, there are still contributions from $S_{+-}^2$ and $S_{-+}^2$,
which result equally in the formation of a flat background, while that
from $S_{-+}^2$ is indeed negligible.

\bibliographystyle{JCAP}
\bibliography{ref}
\end{document}