NORDITA-2021-XXX

Can we observe the QCD phase transition-generated gravitational waves through pulsar timing arrays?

Axel Brandenburg,^{1, 2, 3, 4}^{,} ^{∗} Emma Clarke ^{†},^{4}^{,} ^{‡} Yutong He,^{1, 2}^{,} ^{§} and Tina Kahniashvili^{4, 3, 5, 6}^{,} ^{¶}

1Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden 2Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden

3Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia

4McWilliams Center for Cosmology and Department of Physics,

Carnegie Mellon University, Pittsburgh, PA 15213, USA 5Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia

6Department of Physics, Laurentian University, Sudbury, ON P3E 2C, Canada (Dated: February 23, 2021)

We perform numerical simulations of gravitational waves (GWs) induced by hydrodynamic and hydromagnetic turbulent sources that might have been present at cosmological quantum chromodynamic (QCD) phase transitions. For turbulent energies of about 4% of the radiation energy density, the typical scale of such motions may have been a sizable fraction of the Hubble scale at that time. The resulting GWs are found to have an energy fraction of about 10^{−9} of the critical energy density in the nHz range today and may already have been observed by the NANOGrav collaboration. This is further made possible by our findings of shallower spectra proportional to the square root of the frequency for nonhelical hydromagnetic turbulence. This implies more power at low frequencies than for the steeper spectra previously anticipated. The behavior toward higher frequencies depends strongly on the nature of the turbulence. For vortical hydrodynamic and hydromagnetic turbulence, there is a sharp drop of spectral GW energy by up to five orders of magnitude in the presence of helicity, and somewhat less in the absence of helicity. For acoustic hydrodynamic turbulence, the sharp drop is replaced by a power law decay, albeit with a rather steep slope. Our study supports earlier findings of a quadratic scaling of the GW energy with the magnetic energy of the turbulence and inverse quadratic scaling with the peak frequency, which leads to larger GW energies under QCD conditions.

PACS numbers:

I. INTRODUCTION

Gravitational wave (GWs) astronomy opens a new window to study the physical processes in the early universe. Relic GWs can be sourced by violent processes such as cosmological phase transitions and after generation they propagate almost freely throughout the expansion of the universe that causes the dilution of their strain amplitude and frequency; for a review, see [1] and references therein. On the other hand, the detection of these relic GWs is a challenging task due to their small amplitudes, the specific range of the characteristic frequencies, and astrophysical foregrounds [2]. Despite tremendous advancements in GW detection techniques, the stochastic GW background of cosmological origin remained unobserved.

Recently, the NANOGrav collaboration reported strong evidence for a stochastic GW background [3]. In addition to the possibility of GWs induced by astrophysical sources such as supermassive black holes, the

†Corresponding author; the authors are listed alphabetically.

∗Electronic address: brandenb@nordita.org
‡Electronic address: emmaclar@andrew.cmu.edu

§Electronic address: yutong.he@su.se

¶Electronic address: tinatin@andrew.cmu.edu

NANOGrav data can also be understood as a possible signal from the early universe, such as inflationary GWs [4–8], cosmic strings and domain walls [9–14] phase transitions including the supercooled phase transitions [15], dark phase transitions [16, 17] and quantum chromodynamic (QCD), with axionic string network and QCD axion [18–21], and/or magnetic fields [22] and turbulence [23].

In this paper we present a self-consistent study of the GWs from turbulent sources possibly present at QCD phase transitions. In contrast to Ref. [22], we will not limit the turbulent sources to the presence of magnetic fields [24–30] but rather consider any turbulent source possibly presented at the QCD energy scale. Even if primordial fields are not dynamically strong, turbulence can still develop at QCD energy scales [24, 28, 31–35]; the latent heat they release still gives rise to pressure gradients resulting in macroscopic plasma motions. Given the very high Reynolds number of the primordial plasma, such motions will inevitably decay into turbulence [24, 35]. Particularly important for our work is the earlier finding that the separation and size of nucleation bubbles in a QCD phase transition is a sizeable fraction of the Hubble scale, see Ref. [36] for pioneering works and follow-up papers, [37–47]. Furthermore, the assumption of turbulence being driven by magnetic fields, allows us to avoid the requirement of first order QCD phase transitions, [35].

The paper is organized as follows. We first review basic properties of relic GWs (Sec. II), discuss then the NANOGrav observations (Sec. III), present our numerical approach (Sec. IV) and results (Sec. V) of our simulations, before concluding in Sec. VI. Throughout the paper, we use natural units with h̄ = c = k_{B} = 1. We also set the permeability of free space to unity, i.e., μ_{0} = 1, expressing the electromagnetic quantities in Lorentz-Heaviside units. The Latin indices run i ∈ (1, 2, 3) and define the spatial coordinates, and the Greek indices run λ ∈ (0, 1, 2, 3). We choose the metric signature as (−1, 1, 1, 1).

II. THE EARLY-UNIVERSE GRAVITATIONAL WAVE SIGNAL

GWs correspond to the tensor mode of perturbations δg_{μν} above the spatially flat, homogeneous, and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) background, in the transverse-traceless (TT) gauge^{1} defined through the spatial component h_{ij}^{phys} with a^{2}h_{ij}^{phys} = Λ_{ijlm}δg_{lm}, where a is the scale factor at the physical time t_{phys}. Here and below, super/subscript “phys” denotes physical quantities.

In order to eliminate the expansion-induced dilution from the governing hydromagnetic equations, we use rescaled quantities together with the conformal time t, defined through dt = dt_{phys}/a, which reduces the metric tensor to the Minkowski form. The background expansion of the universe during the radiation-dominated epoch is governed by the (dominant) radiation energy density E_{rad} = π^{2}g(T)T^{4}/30, where g(T) is the effective number of relativistic degrees of freedom at temperature T. In the epoch(s) of interest, the expansion of the universe is fully governed by radiation, and the Hubble parameter H ≡ a^{−1}da/dt_{phys} = a^{−2}da/dt is given through H^{2}(t) = (8πG/3) E_{rad}(t), where G is Newton’s gravitational constant and E_{rad}(t) denotes the total energy density of radiation (including all relativistic components).

In order to connect physical and comoving variables and to determine the scaling of physical quantities, we compute the ratio of the scale factor today, a_{0} = a(t_{0}) (here and below, “0” denotes the present moment), to that at the time t_{∗} (at the temperature T_{∗} at which the source becomes active and the gravitational signal is generated) corresponding to the start of the simulation. We assume the adiabatic expansion of the universe, such that g_{S}(T)T^{3}a^{3}(T) is constant, where g_{S}(T) is the

1 The TT gauge is determined by the TT projection tensor Λ_{ijkl} = P_{ik}P_{jl}− \frac{1}{2} P_{ij}P_{kl}, where the P_{ij} is a transverse operator (∂_{i}P_{ij} = 0), defined as P_{ij} = δ_{ij} − ∂_{i}∂j/∇^{2}, where δ_{ij} is the Kronecker delta, ∂_{λ} ≡ ∂/∂x^{λ} denotes the partial derivative in respect of x^{λ} coordinate, and (∇ defines the vector differential operator with the components equal to ∇_{i} ≡ ∂_{i}, i.e., ∇^{2} is the Laplacian in respect of spatial coordinates, see for more details Chapter 1 (1.2) of [48]

number of adiabatic degrees of freedom at temperature T. At high enough temperatures (T > 1 MeV), we have g_{S}(T) = g(T) [49]. Note, that our consideration below is valid for any time period during the radiation dominated epoch. However, we will be focused on the time period around the QCD energy scale (150 MeV). We also normalize the scale factor a_{∗} ≡ a(t_{∗}) = 1, which differs from the usual convention a_{0} = 1. Entropy conservation leads to

\frac{a_{0}}{a_{∗}} = 10^{12} (\frac{g_{S}(T_{∗})}{15})^{1/3}(\frac{T_{∗}}{150 MeV} ),  (1)

where we have used T_{0} = 2.73 K and g_{S}(T_{0}) = 3.91, while at the QCD energy scale we have g_{S}(T_{∗}) ≈ 15 [49]. The degrees of freedom at QCD is approximate due to uncertainty in the exact temperature of the QCD transition and knowledge of the standard model (see discussions in [45, 50]). However, as (a_{0}/a_{∗}) ∼ g_{S}(T_{∗})^{1/3}, small deviations in g_{S}(T_{∗}) will not significantly impact our results. The GW equation in physical time and space coordi-

nates is given by

(∂_{t_{phys}}^{2} + 3H∂_{t_{phys}} − ∇_{phys}^{2})h_{ij}^{phys} = 16πGT_{ij,phys}^{TT}, (2)

where the TT superscript denotes the TT projection of the stress-energy tensor such that T_{ij,phys}^{TT} = Λ_{ijlm}T_{lm}^{phys}. To make the connection with observations, we define the characteristic strain, h_{c}(t), which obeys h_{c}^{2}(t) = ⟨(h_{ij}^{phys}(x, t))^{2}⟩ /2, where angle brackets denote volume averaging in physical space, and the physical energy density E_{GW}^{phys}(t) carried by the GWs is given by [48]

E_{GW}^{phys}(t) = \frac{1}{32πG} ⟨(∂_{t_{phys}}h_{ij}^{phys}(x, t))^{2}⟩ .  (3)

It is then expressed in terms of today’s frequency f = k/(2πa_{0}) that corresponds to the time Fourier transform Q(t) = ∫_{−∞}^{∞} df Q(f)e^{−2πft} (and Q(f) = 2π ∫_{−∞}^{∞} dt Q(t)e^{−2πft}) [48].

The relic GW signal strength today is given through the normalized GW energy density parameter Ω_{GW}(f) reduced by the factor (H_{∗}/H_{0})^{2}(a_{∗}/a_{0})^{4}, where H_{∗} is the Hubble parameter at t_{∗}. This accounts for the dilution of the GW energy density parameter with the expansion of the universe and renormalizes the GW energy density by the critical energy density at the present time, E_{crit}^{0} = (3H_{0}^{2})/(8πG), where H_{0} = 100 h_{0} km s^{−1} Mpc^{−1} ≃ 3.241 × 10^{−18} h_{0} s^{−1} is the present value of the Hubble parameter. A frequency of particular interest is the frequency f_{∗} corresponding to the Hubble horizon scale at t_{∗}:

f_{∗} = \frac{a_{∗}H_{∗}}{a_{0}} ≃ (1.8×10^{−8} Hz)(\frac{g_{∗}}{15})^{1/6}(\frac{T_{∗}}{150 MeV} ). (4)

As discussed above, there is a variety of possible sources of a stochastic GW background in the nHz frequency range, accessible to Pulsar Timing Arrays (PTAs), see Section III for more details and Ref. [51] for a review and references therein, and these sources include a cosmic population of supermassive black hole binaries (SMBHBs) [51, 52], cosmic strings [53–55], inflationary GWs [4, 56],^{2} and phase transitions in the early universe [19, 60–62]. In particular, if we assume that GWs could be sourced by bubble collisions at a phase transition, we expect the frequency of the GWs to be related to the bubble size. We consider that the bubble length scale is the Hubble horizon H_{∗} at generation divided by the total number of phase transition bubbles N_{b}. Then, for the QCD phase transitions, the frequency is given by

f_{∗} ≃ (1.1 × 10^{−7} Hz)(\frac{g_{∗}}{15})^{1/6}(\frac{T_{∗}}{150 MeV} )(\frac{N_{b}}{6} ), (5)

where we have normalized to 6 bubbles expected at the QCD phase transition [45].

III. NANOGRAV DATA

A pulsar is a highly magnetized and rapidly rotating neutron star that emits a beam of electromagnetic radiation along its magnetic axis [63]. The times of arrival (TOA) of these pulses are extremely regular and can be predicted very accurately over long times [64]. The presence of a GW passing between the observer and pulsar shifts the pulse TOA proportional to the amplitude of the GW [65]. By monitoring the fluctuations in the TOA of radio pulses from millisecond pulsars (see, e.g., Ref. [66] and references therein) international PTA missions^{3} aim to probe a stochastic GW background.

The maximum sensitivity of a PTA experiment is limited by the total observation time. That is, the lowest detectable frequency is on the order of the inverse of the time span of the data (e.g., f ∼ nHz for datasets spanning ∼ 10 years) [68]. Furthermore, data sampling (i.e., pulsars are usually observed on the order of weeks [65]) limits the maximum detectable frequency. The NANOGrav 12.5-year data is sensitive to GW frequencies between approximately 2.5 nHz and 1 μHz [69].

PTA measurements typically characterize a stochastic GW background in terms of its characteristic strain spectrum h_{c}(f) fitted with a power-law dependence on frequency [3],

h_{c}(f) = A_{CP}(\frac{f}{f_{yr}} )^{α_{CP}},  (6)

2 The quantum mechanical fluctuations during the inflationary

epoch induces GWs via parametric resonance [57–59]. 3 The International Pulsar Timing Array (IPTA) is a consortium

of consortia, comprised of the European Pulsar Timing Array (EPTA), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav), and the Parkes Pulsar Timing Array (PPTA) [67].

where the subscript “CP” denotes a common-spectrum (CP) process (common to the observed pulsars), the spectral index α_{CP} depends on the source of the stochastic GW background, and A_{CP} is the strain amplitude at a reference frequency of f_{yr} = 1 yr^{−1}. This choice of reference frequency is arbitrary and does not affect the ability to detect a GW signal.

The energy density spectrum of the GW background today expressed in terms of the characteristic strain spectrum is given by [65]

Ω_{GW}(f) = \frac{2π^{2}}{3H_{0}^{2}} f^{2}h_{c}^{2}(f) = Ω_{GW}^{yr}(\frac{f}{f_{yr}} )^{5−γ_{CP}}, (7)

where we have used Eq. (6) in the second term on the right-hand side, γ_{CP} = 3 − 2α_{CP} and Ω_{GW}^{yr} ≡ 2π^{2}A_{CP}^{2}f_{yr}^{2}/(3H_{0}^{2}). The quantity h_{0}^{2}Ω_{GW}(f) is typically considered in order to remove the uncertainty in the value of H_{0}.

The NANOGrav collaboration reports joint A_{CP} −γ_{CP} posterior distributions [3]. Posteriors for a commonspectrum process in the NANOGrav 12.5-year data were recovered with four models: free-spectrum, broken power law, 5-frequency power law, and 30-frequency power law. The fits were performed for frequencies f ∈ [2.5 × 10^{−9}, 7 × 10^{−8}] Hz, with the exception of the 5-frequency power law, which was fit to the five lowest frequency bins. The four lowest frequency bins have the strongest response to the presence of a GW background (see Figure 13 of Ref. [3]). Thus, the 5-frequency power law was fit within the signal-dominated frequency range (approximately f ∈ [2.5 × 10^{−9}, 1.2 × 10^{−8}] Hz). Figure 1 of Ref. [3] shows the 1σ and 2σ posterior contours for the amplitude A_{CP} and spectral slope γ_{CP}.

Fig. 1 shows the NANOGrav detection expressed in terms of h_{0}^{2}Ω_{GW}(f) as given by Eq. (7). The shaded regions show the 2σ confidence contours of the A_{CP} − γ_{CP} parameter space in terms of f and h_{0}^{2}Ω_{GW}(f) for frequencies from 2.5–100 nHz (i.e., the NANOGrav 12.5year sensitivity range); see Ref. [3] for more detail.

IV. GRAVITATIONAL WAVE GENERATION

As mentioned in the introduction (see also Ref. [51]), low frequency GWs can originate from various astrophysical foreground sources (white dwarfs, SMBH mergers, etc), and from relic sources related to inflation and cosmic strings, for example, and in particular, from phase transition-generated turbulence and primordial magnetic fields. We focus here on the latter two. Turbulence and/or magnetic fields would only be generated during a limited amount of time before they would decay. The decay process itself remains highly turbulent and could affect GW production. Let us therefore begin with some general remarks about turbulent decay.

FIG. 1: NANOGrav 12.5-year data set 2σ confidence contours for the posteriors of a common-spectrum process (see Ref. [3] Figure 1) expressed in terms of the GW energy density h_{0}^{2}Ω_{GW}(f) and frequency f. This is shown over the NANOGrav 12.5-year sensitivity range of 2.5–100 nHz. The three models used to fit the process include a: broken power law (blue), 5-frequency power law (orange), and 30-frequency power law (green).

A. Gravitational Waves from Turbulent Sources

Using scaled quantities h_{ij} = ah_{ij}^{phys} and T_{ij}^{TT} = a^{4}T_{ij,phys}^{TT}, together with a(t) ∝ t in the radiation dominated epoch, the GW equation takes the form

(∂_{t}^{2} − ∇^{2}) h_{ij} = \frac{16πG}{a}T_{ij}^{TT}.  (8)

To obtain the GW equation in Fourier (wavenumber) space, we use the Fourier transforms and the polarization r = (+, ×) decomposition of the tensor metric perturbations and stress energy tensor projected onto the TT gauge (i.e., Q_{ij}(k, t) = ∑_{r=+,×} e_{ij}^{r}(k̂)Q_{ij}(k, t), where

e_{ij}^{+}(k̂) and e_{ij}^{×}(k̂) are the polarization tensors with k̂ the unit vector, and k = ak_{phys} is the rescaled wavenumber).^{4}

As in earlier work [72, 73], we use normalized conformal time, t̄ = t/t_{∗}, where t_{∗} = H_{∗}^{−1} is our starting time, and a_{∗} = 1 has been chosen. Therefore, a = t̄. We also use the scaled wave vector, k̄ = k/H_{∗}, and a scaled normalized stress, T̄_{+/×}^{TT} = T_{+/×}^{TT}/E_{rad}^{∗}. The GW equation

4 We use the spatial Fourier transform convention: Q(x, t) = ∫ \frac{d^{3}k}{(2π)^{3}} e^{−ik·x}Q(k, t) and Q(k, t) = ∫ d^{3}xe^{−ik·x}Q(x, t). The transverse operator P_{ij} in the Fourier space is given P_{ij}(k̂) = δ_{ij} − k̂_{i}ĵ_{k} and the TT projection operator is Λ_{ijkl}(k̂) = P_{ik}(k̂)P_{jl}(k̂) − \frac{1}{2} P_{ij}(k̂)P_{kl}(k̂), correspondingly. The polarization tensors e_{ij}^{+}(k̂) and e_{ij}^{×}(k̂) can be written as e_{ij}^{+}(k) = ê_{i}^{1}ê_{j}^{1} − ê_{i}^{2}ê_{j}^{2} and e_{ij}^{×}(k) = ê_{i}^{1}ê_{j}^{2} + ê_{i}^{2}ê_{j}^{1}, where ê^{1} and ê^{2} are unit vectors that are orthogonal to k̂ and each other; see Chapter 1 (1.2) of Ref. [48]. and Ref. [72] for further detail.

can then be written in the form [72, 73]

(∂2 ; t̄ + k̄^{2}) h_{+/×}(k, t) = \frac{6}{t̄}T̄_{+/×}^{TT}(k, t),  (9)

but from now on we omit all overbars.

Throughout this paper, all numerical results will usually be presented as the scaled variables introduced above. In particular, we quote the rms strain, h_{rms} = ⟨h^{2}⟩^{1/2}, where h^{2} = h_{+}^{2} + h_{×}^{2} = (h_{ij})^{2}/2, and likewise for the scaled GW energy, E_{GW} = ⟨ḣ^{2}⟩ /6, where ḣ_{+/×} = ∂_{t}h_{+/×} with ḣ^{2} ≡ ḣ_{+}^{2} + ḣ_{×}^{2}; see Ref. [72, 73] for additional subdominant terms that are applied in the calculations. We sometimes also quote the (frequency dependent) characteristic amplitude of the physical strain measured today, h_{c}(f) = h_{rms}/a_{0}; see Sec. II.

B. Turbulent Sources

Turbulent flows in the early universe can be modeled by solving the hydromagnetic equations for the density ρ, the velocity u, and the magnetic field B with ∇ · B = 0, adopting an ultrarelativistic equation of state in an expanding universe using conformal time and comoving variables [70, 71] with a forcing term F in the induction equation for B

\frac{∂ ln ρ}{∂t} = −\frac{4}{3} (∇ · u + u · ∇ ln ρ) + \frac{1}{ρ} [u · (J × B) + ηJ^{2}] ,
\frac{∂u}{∂t} = −u · ∇u + \frac{u}{3} (∇ · u + u · ∇ ln ρ) + \frac{2}{ρ}∇ · (ρνS)

−\frac{1}{4}∇ ln ρ − \frac{u}{ρ} [u · (J × B) + ηJ^{2}] + \frac{3}{4ρ}J × B,

\frac{∂B}{∂t} = ∇ × (u × B − ηJ + F), J = ∇ × B.

We recall that the conformal time t is normalized to unity at the time t_{∗} of magnetic field generation, ρ is in units of the initial value, u is in units of the speed of light, and the magnetic energy density B^{2}/2 is measured in units of the radiation density at the time of generation. Furthermore, S_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) − \frac{1}{3}δ_{ij}∇ · u are the components of the rate-of-strain tensor with commas denoting partial derivatives, J is the current density, ν is the kinematic viscosity, and η is the magnetic diffusivity. The electromotive force, F, is used to model magnetic field generation with

F(x, t) = Re[N f̃(k) exp(ik · x + iφ)],  (10)

where the wave vector k(t) and the phase φ(t) change randomly from one time step to the next. This forcing function is therefore white noise in time and consists of plane waves with average wavenumber k_{f} such that |k| lies in an interval k_{f} − δk/2 ≤ |k| < k_{f} + δk/2 of width δk. Here, N = f_{0}/δt^{1/2} is a normalization factor, where δt is the time step and f_{0} is varied to achieve a certain magnetic field strength after a certain time, and f̃(k) = (k × e)/[k^{2} − (k · e)^{2}]^{1/2} is a nonhelical forcing function. Here, e is an arbitrary unit vector that is not aligned with k. Note that |f|^{2} = 1. Following earlier work, the forcing is only enabled during the time interval 1 ≤ t ≤ 2. The kinetic and magnetic energy densities are defined as E_{K}(t) = ⟨ρu^{2}⟩ /2 and E_{M}(t) = ⟨B^{2}⟩ /2, respectively.

The vigor of turbulence is characterized by the Reynolds number, Re = u_{rms}/νk_{f}, where u_{rms} is the maximum rms velocity. It can only be determined a posteriori from the velocity resulting from the magnetic field through the Lorentz force. For all our runs, we use η = ν.

C. Turbulent Decay Laws

Turbulence is known to decay in power-law fashion [74, 75] such that the magnetic energy E_{M}(t) decays with time t like Δt^{−p} and the correlation length ξ_{M}(t) increases like Δt^{q}, where Δt = t−t_{off} is the time interval after the forcing has been turned off. The exponents p and q are positive and depend on the physical circumstances (magnetically or kinetically dominated turbulence), and whether or not there is magnetic helicity. In helical turbulence, for example, one finds p = q = 2/3, while for non-helical magnetically dominated turbulence one finds p = 1 and q = 1/2, although other variants are sometimes possible [76, 77].

In this paper, we are specifically interested in the dependence of the decay behavior on the forcing wave number k_{f} of the turbulence while it was still being driven. The parameter k_{f} enters through the prefactor in the decay law.

Furthermore, Δt^{−p} would become infinite for p > 0 and Δt = 0 (when t = t_{off}). The singularity of Δt^{−p} at t = 0 is a consequence of a simplified description at the initial time moment. For this reason, it is convenient to express the decay laws as

E_{M}(t) = E_{M}^{max} (1 + Δt/τ)^{−p},  (11)

where τ is the turnover time, which we will treat as an empirical parameter that we expect to be of the order of (v_{A}k_{f})^{−1}, where v_{A} = (3E_{M}^{max}/2)^{1/2} is the Alfvén speed, evaluated at the time when E_{M} reaches its maximum value E_{M}^{max}. In some simulations of purely hydrodynamic turbulence, we replace E_{M} by E_{K} in Eq. (11) and use τ = (u_{rms}k_{f})^{−1} with u_{rms} = (2E_{K})^{1/2} as the nominal turnover time. We recall here that we are using nondimensional variables where the radiation energy density is unity.

For all our simulations, we choose t_{off} = 2, i.e., turbulence is being driven for one Hubble time during 1 ≤ t ≤ 2. In the following, we vary k_{f} between 2 and 60. For k_{f} = 60, we find that τ is shorter than a Hubble time, but in all other cases, it exceeds it by up to factors between ten (in the nonhelical cases) and a hundred (in helical cases).

We arrange the strength of the forcing f_{0} such that E_{M} is similar for different values of k_{f}. This allows us then to determine the resulting GW energy solely as a function of k_{f}. For small values of k_{f}, the turbulence may not be able to reach a statistically steady state by the time t_{off}, when the driving is turned off.

It is therefore necessary to adjust f_{0} for each value of k_{f} separately. Once we have two values of E_{M}^{max} that are close enough to the target strength, we determine the desired forcing strength through linear interpolation. We also consider the case of different values of f_{0} for a fixed value of k_{f} (Runs noh5,6 and Runs hel5,6).

V. NUMERICAL SIMULATIONS

We solve the governing equations using the Pencil Code [78], where the GW solver has already been implemented [72]. We consider a cubic domain of side length 2π/k_{1}, where k_{1} is the smallest wave number in the domain. We choose k_{1} = k_{f}/6, so that the scale separation between the initial spectral peak and the lowest wave number in the domain is six. In the following, we discuss the results for different values of k_{f}. The temporal growth of E_{M}(t) is similar for small values of k_{f}; see Fig. 2(a) and (b), where we compare the evolution of E_{M} and E_{GW} for the nonhelical and helical cases. The parameters of those runs are listed in Tables I and II (for nonhelical and helical runs). The numerical resolution is 512^{3} mesh points, except for run noh1, where we use 1024^{3} mesh points. Unless specified otherwise, we use ν = η = 5 × 10^{−5}.

In Table I, we have quoted the values of E_{GW}^{sat} and h_{rms}^{sat} obtained at the end of the simulation. To compute the relic observable h_{0}^{2}Ω_{GW} at the present time, we have to multiply E_{GW}^{sat} by a factor (H_{∗}/H_{0})^{2}(a_{∗}/a_{0})^{4}; see Refs. [72, 73] for details. Using g_{∗} = 15 and T_{∗} = 150 MeV, we find H_{∗} = 1.8×10^{4} s^{−1}, and thus this factor is ≈ 3×10^{−5}. The largest value of E_{GW}^{sat} quoted in Table I is 3.5×10^{−4} and corresponds therefore to h_{0}^{2}Ω_{GW} ≈ 10^{−8}. Likewise, the values of h_{rms}^{sat} in Table I have to be multiplied by a_{0}^{−1} ≈ 10^{−12} to obtain the observable h_{c} at the present time; see Eq. (1). Again, the largest value of h_{rms}^{sat} = 5 × 10^{−2} corresponds therefore to the observable h_{c} = 5 × 10^{−14}.

To simplify the comparisons, we have arranged the forcing amplitude f_{0} such that E_{M}^{max} is similar in certain cases. The values of E_{M}^{max} listed in the upper block of Tables I and II (for nonhelical and helical hydromagnetic turbulence, respectively) are around 0.038 and correspond to 0.8 μG. The growth phase of E_{M}(t) is similar, but the decay is significantly slower when k_{f} is smaller. The GW energy saturates at a value E_{GW}^{sat} some time after E_{M}(t) has reached its maximum, and is smaller for larger values of k_{f}.

It is important to realize that in all four cases, the decay of the magnetic energy follows an approximate power law decay, as given by Eq. (11). To see this, we the plot in Fig. 3 the evolution of E_{M} versus t − 2 in a doublelogarithmic representation. The parameters p and τ describe the decay and are also listed in Table I.

FIG. 2: Evolution of E_{M}(t) and E_{GW}(t) for nonhelical (left) and helical (right) cases. Orange, black, blue, and red are for k_{f} = 2, 6, 20, and 60, respectively.

TABLE I: Summary of runs with nonhelical turbulence.

Run k_{f} k_{1}  f_{0}  p τ  E_{M}^{max}  E_{GW}^{sat}  h_{rms}^{sat}  B [μG] h_{0}^{2}Ω_{GW}(f)  h_{c}
noh1 2 0.3 1.9 × 10^{−1} 1.0 16 3.83 × 10^{−2} 3.53 × 10^{−4} 4.83 × 10^{−2} 0.78 1.09 × 10^{−8} 4.83 × 10^{−14}
noh2 6 1 6.0 × 10^{−2} 1.0 4.5 3.75 × 10^{−2} 5.61 × 10^{−5} 7.06 × 10^{−3} 0.78 1.73 × 10^{−9} 7.07 × 10^{−15}
noh3 20 3 2.3 × 10^{−2} 1.3 2.0 3.81 × 10^{−2} 1.11 × 10^{−5} 1.15 × 10^{−3} 0.78 3.44 × 10^{−10} 1.15 × 10^{−15}
noh4 60 10 1.0 × 10^{−2} 1.4 0.43 3.93 × 10^{−2} 2.62 × 10^{−6} 1.65 × 10^{−4} 0.79 8.10 × 10^{−11} 1.65 × 10^{−16}
noh5 2 0.3 1.0 × 10^{−1} — — 1.06 × 10^{−2} 2.70 × 10^{−5} 1.40 × 10^{−2} 0.41 8.37 × 10^{−10} 1.40 × 10^{−14}
noh6 2 0.3 3.0 × 10^{−1} — — 9.48 × 10^{−2} 2.08 × 10^{−3} 1.02 × 10^{−1} 1.2 6.42 × 10^{−8} 1.02 × 10^{−13}
noh7 6 1 2.0 × 10^{−2} — — 4.63 × 10^{−3} 6.56 × 10^{−7} 8.10 × 10^{−4} 0.27 2.03 × 10^{−11} 8.11 × 10^{−16}

TABLE II: Similar to Table I, but for helical turbulence.

Run k_{f} k_{1}  f_{0}  p τ  E_{M}^{max}  E_{GW}^{sat}  h_{rms}^{sat}  B [μG] h_{0}^{2}Ω_{GW}(f)  h_{c}
hel1 2 0.3 1.9 × 10^{−1} 0.67 100 3.90 × 10^{−2} 4.85 × 10^{−4} 4.33 × 10^{−2} 0.79 1.50 × 10^{−8} 4.33 × 10^{−14}
hel2 6 1 5.6 × 10^{−2} 0.67 20 3.81 × 10^{−2} 5.05 × 10^{−5} 4.69 × 10^{−3} 0.78 1.56 × 10^{−9} 4.69 × 10^{−15}
hel3 20 3 2.0 × 10^{−2} 0.67 4.0 3.96 × 10^{−2} 7.26 × 10^{−6} 6.66 × 10^{−4} 0.80 2.24 × 10^{−10} 6.66 × 10^{−16}
hel4 60 10 6.5 × 10^{−3} 0.67 0.50 3.76 × 10^{−2} 8.15 × 10^{−7} 7.18 × 10^{−5} 0.78 2.52 × 10^{−11} 7.18 × 10^{−17}
hel5 2 0.3 1.0 × 10^{−1} — — 1.06 × 10^{−2} 3.61 × 10^{−5} 1.08 × 10^{−2} 0.41 1.12 × 10^{−9} 1.08 × 10^{−14}
hel6 2 0.3 3.0 × 10^{−1} — — 9.85 × 10^{−2} 3.07 × 10^{−3} 1.12 × 10^{−1} 1.3 9.49 × 10^{−8} 1.12 × 10^{−13}
hel7 6 1 2.0 × 10^{−2} — — 4.93 × 10^{−3} 8.33 × 10^{−7} 6.26 × 10^{−4} 0.28 2.58 × 10^{−11} 6.26 × 10^{−16}
hel8 6 1 1.0 × 10^{−1} — — 1.20 × 10^{−1} 5.09 × 10^{−4} 1.59 × 10^{−2} 1.4 1.57 × 10^{−8} 1.59 × 10^{−14}

Our results confirm that the turbulence decays more slowly for large values of τ, or small values of k_{f}. As already found from earlier simulations [73], the GW energy generally decreases with increasing k_{f}. This is seen more

FIG. 3: Similar to Fig. 2(a) and (b), but in a double-logarithmic representation where E_{M} is now plotted versus Δt ≡ t − 2, the time after which the electromagnetic source is turned off.

TABLE III: Comparison of nonhelical magnetic turbulence (mag) with irrotational (irro) and vortical (vort) turbulence.

Type  f_{0}  ν  E_{M}^{max}  E_{GW}^{sat}  h_{rms}^{sat}  B [μG] h_{0}^{2}Ω_{GW}(f)  h_{c}
magnetic 1.9 × 10^{−1} 5.0 × 10^{−5} 3.83 × 10^{−2} 3.53 × 10^{−4} 4.83 × 10^{−2} 0.78 1.09 × 10^{−8} 4.83 × 10^{−14}
vortical  3.8 × 10^{−1} 1.0 × 10^{−2} 4.21 × 10^{−2} 8.81 × 10^{−4} 8.26 × 10^{−2} 0.82 2.73 × 10^{−8} 8.27 × 10^{−14}
irrotational 7.0 × 10^{−1} 2.0 × 10^{−2} 4.26 × 10^{−2} 8.30 × 10^{−4} 7.95 × 10^{−2} 0.83 2.57 × 10^{−8} 7.96 × 10^{−14}

FIG. 4: (a) E_{GW} versus E_{M}/k_{f}; the straight line shows E_{GW} = 5.2 × 10^{−4} (E_{M}/k_{f})^{1/2}. (b) Positions of our runs in a diagram showing E_{GW}^{sat} versus E_{M}^{max}. For orientation the old data points of the Ref. [73] are shown as gray symbols. The open red (filled blue) symbols are for the helical (nonhelical) runs. The green symbols refer to the two hydromagnetic runs of Table III.

clearly in a diagram of E_{GW} versus E_{M}/k_{f}; see Fig. 4(a).

For k_{f} = 2, we have performed additional simulations with smaller and larger values f_{0}, both with and without helicity. The resulting values of E_{GW}^{sat} obey quadratic scaling of the form

E_{GW}^{sat} = (qE_{M}^{max}/k_{f})^{2}  (12)

with a coefficient q = 1.1; see the straight line in Fig. 4(a). Only the data point for k_{f} = 60 is slightly above the line represented by Eq. (12). This could be an artefact of our Reynolds numbers still not being large enough in our simulations, especially for large value of k_{f}.

To compare with earlier work, we show in Fig. 4(b) the positions of our runs in a E_{GW}^{sat} versus E_{M}^{max} diagram. For orientation, we also show the data points from Ref. [73]. We see that the new data points are well above the older ones of Ref. [73]. This is mainly a consequence of using here smaller values of k_{f} (2–60, compared to 600 in Ref. [73]). For k_{f} = 2 and k_{1} = 0.3, we show here the results for hydrodynamic runs using irrotational and vortical forcings; see the green symbols in Fig. 4(b). Those runs are listed in Table III and compared with the nonhelical magnetic turbulence run ‘noh1’.

In Fig. 5, we plot the resulting present-day GW energy and strain spectra for our four runs with k_{f} = 2, 6,

FIG. 5: h_{0}^{2}Ω_{GW}(f) and h_{c}(f) at the present time for all four runs presented in Table I, for the nonhelical (left) and helical (right) runs. The 2σ confidence contour for the 30-frequency power law of the NANOGrav 12.5-year data set is shown in gray.

FIG. 6: Magnetic energy spectra for the (a) nonhelical and (b) helical cases.

20, and 60, both without and with helicity in the driving function F. The first two cases with k_{f} = 2 and 6 lie well within the frequency and amplitude range accessible to NANOGrav. In all cases, the spectra show a sharp drop slightly above the peak frequency. This is a consequence of the rapid temporal growth of the spectra, which leads to a correspondingly large growth at the peak frequency, while at higher frequencies, the spectrum settled at values that were determined by somewhat earlier times when the energy was still weaker.

At frequencies below the peak, we now find a spectrum that is even shallower than the h_{0}^{2}Ω_{GW}(f) ∝ f spectrum found already earlier [73]. A spectrum shallower than proportional to f, such as the present f^{1/2} spectrum, could perhaps be explained by the finite size of the computational domain; see Ref. [79], who found

FIG. 7: Similar to Fig. 5, but comparing vortical (red) and irrotational turbulence (blue) with MHD turbulence (black).

even a f^{−1/2} spectrum for Ω_{GW}(f). Alternatively, the shallower spectrum might well be physical, or at least significantly extended over a substantial frequency interval below the peak frequency, for example due to inverse cascading in helical [80] and nonhelical [81] cases.

In the absence of sources, a Ω_{GW}(f) ∝ f^{α} spectrum implies h_{c}(f) ∝ f^{α/2−1} for arbitrary spectral indices α. For α = 1/2, we would thus expect h_{c}(f) ∝ f^{−3/4}. However, the observed strain spectrum, h_{c}(f) ∝ f^{−1/2}, seems to agree with that found previously from numerical simulations [73]. However, looking more carefully at the strain spectrum for k_{f} = 60, we see a h_{c}(f) ∝ f^{−3/4} spectrum is actually compatible with the simulation; see the corresponding dashed-dotted line in Fig. 5(c). This agreement is probably related to the fact that the turnover time is shorter for the run with k_{f} = 60, compared with those at smaller values (i.e., longer turbulence driving time will allow for more efficient inverse cascading).

In the runs with helicity, we do find h_{0}^{2}Ω_{GW}(f) ∝ f, together with a slight enhancement just before reaching a maximum. The subsequent decay for larger values of f is much steeper in the case with helicity than without. Furthermore, in h_{c}(f) we see a sharper drop to the right of the maximum than in simulations without helicity. These differences in the spectra for helical and nonhelical cases are surprisingly strong and might allow us to infer the presence of magnetic helicity once such a spectrum is detected.

It is important to note that the h_{0}^{2}Ω_{GW} ∝ f spectra in Figs. 4(a) and (c) show an increase towards smaller k_{f}. This is to be expected from Eq. (12), but it was not included in the sketch of Ref. [22]; see their Fig. (1). By contrast, in their Eq. (4), an effectively cubic dependence on the magnetic energy was motivated.

The underlying magnetic energy spectrum is shown in Fig. 6(a) for nonhelical and in Fig. 6(b) helical cases where k_{f} is ranging from 2 to 60. Those are averaged spectra obtained by averaging over the time interval 15 ≤ t ≤ 20. In the nonhelical case, the amplitude of the spectrum is smaller for larger values of k_{f}, because here the energy has decayed more rapidly. In the helical case, the spectra have approximately the same height for all values of k_{f}. This is because the height of the spectrum is related to the helicity, which is conserved. For small values of k_{f}, the spectrum has a more extended subinertial range. This is because the turnover time is larger and there was not enough time for the inverse cascade to produce energy and small values of k.

Finally, we compare the results for two types of purely hydrodynamic turbulence with vortical and irrotational forcings of Table III. The result is shown in Fig. 7. All these cases are for plane wave forcings. For irrotational forcing, we do not see the sharp drop-off of spectral power for frequencies above the peak value as in the vortical case. This suggests that in the inertial range of irrotational turbulence, there is still some power to contribute to GW driving compared with the vortical case, where this is almost not possible at all. However, the spectrum in the irrotational case shows a fairly steep spectrum proportional to f^{−7}, so the effect on GW production is here also rather weak. Nevertheless, the spectral form of the peak might give interesting diagnostic clues about the nature of turbulent driving at the time of GW production.

We mention in passing that in earlier work, it was found that irrotational turbulence is much more efficient in driving GWs than vortical turbulence [73]. Remarkably, here this is no longer the case and vortical and irrotational turbulence have rather similar GW energies. This could be related to the small value of k_{f}, possibly combined with a comparatively short time of driving. However, to clarify this further, more targeted numerical experiments would need to be performed.

VI. CONCLUSIONS

In the present work, we have shown that the magnetic stress from MHD turbulence with scales comparable to the cosmological horizon scale at the time of the QCD phase transition can drive GWs in the range accessible to NANOGrav, if the magnetic energy density is 3–10% of the radiation energy density. The low-frequency tail below the peak frequency at 10 nHz or so is shallower in

FIG. 8: Upper left: NANOGrav 12.5-year contours, gray shaded region corresponds to possible peak frequencies of Ω_{GW} for source described in this paper, for 1-6 bubbles at QCD; Upper Right: cosmic strings: four average power spectrum models described in [82] are mono (orange), kink (green), cusp (blue), and a spectrum computed from a simulated gravitational backreaction model (‘BOS’) (magenta); Lower Left: Inflation [4] n_{T} − r parameter space corresponding within parameter space of NANOGrav 12.5-year results; Lower Right: GWs from SMBHBs are expected to have a spectral index of γ_{CP} = 13/3, which corresponds to A_{CP} = [1.4, 2.7] × 10^{−15} from the NANOGrav 5-frequency and broken power laws.

the nonhelical case than in the helical one, i.e., ∝ f^{1/2} compared to ∝ f. Both scaling are, however, shallower than what was expected based on earlier analytical calculations. Also the inertial range spectrum above the peak is shallower without helicity than with, but here, both spectra are steeper than what is expected if the GW spectrum was a direct consequence of the MHD turbulence spectrum. The reason for this is primarily the relatively short time of turbulent driving (one Hubble time). This is short compared with the turnover time which, for our runs with the smallest k_{f} of two, is much longer: 16 (100) Hubble times for our runs without (with) helicitly. Therefore, there was not enough time to fully establish the GW spectrum at high wave numbers. For our earlier runs with larger values of k_{f}, this effect was less pronounced than for smaller values of k_{f}, but it is still quite noticeable, especially in the helical case where forward cascading is weaker than in the nonhelical case.

Our work has led to new insights regarding the possibility of using an observed GW spectrum for making statements about the nature of the underlying turbulence in the early universe. One is the already mentioned slope of the subinertial range spectrum. Another is the position of the peak of the spectrum. Finally, there is the strength of the drop of the spectral power for frequencies above the peak frequency, and the subsequent slope after the drop, which is most likely too small to be detectable. This, however, depends on the duration of turbulent driving and could be higher if the driving time was longer. The specific features of the spectrum near the peak are different for helical and nonhelical turbulence. This could, in principle, give information about the presence of parity violation, when would also lead to circularly polarized GWs.

Data availability—The source code used for the simulations of this study, the Pencil Code, is freely available from Ref. [78]. The simulation setups and the corresponding data are freely available from Ref. XXX.

Acknowledgments

Support through the Swedish Research Council, grant 2019-04234, and Shota Rustaveli GNSF (grant FR/18-1462) are gratefully acknowledged. 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.

[1] C. Caprini and D. G. Figueroa, “Cosmological Backgrounds of Gravitational Waves,” Class. Quant. Grav. 35, 163001 (2018) [arXiv:1801.04268 [astro-ph.CO]].

[2] J. D. Romano and N. J. Cornish, “Detection methods for stochastic gravitational-wave backgrounds: a unified treatment,” Living Rev. Rel. 20, 2 (2017) [arXiv:1608.06889 [gr-qc]].

[3] Z. Arzoumanian et al. [NANOGrav], “The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background,” Astrophys. J. Lett. 905, L34 (2020) [arXiv:2009.04496 [astro-ph.HE]].

[4] S. Vagnozzi, “Implications of the NANOGrav results for inflation,” Mon. Not. Roy. Astron. Soc. 502, L11 (2021) doi:10.1093/mnrasl/slaa203 [arXiv:2009.13432 [astroph.CO]].

[5] Z. Zhou, J. Jiang, Y. F. Cai, M. Sasaki and S. Pi, “Primordial black holes and gravitational waves from resonant amplification during inflation,” Phys. Rev. D 102, 103527 (2020) [arXiv:2010.03537 [astro-ph.CO]].

[6] H. W. H. Tahara and T. Kobayashi, “Nanohertz gravitational waves from NEC violation in the early universe,” Phys. Rev. D 102, 123533 (2020) [arXiv:2011.01605 [grqc]].

[7] S. Kuroyanagi, T. Takahashi and S. Yokoyama, “Bluetilted inflationary tensor spectrum and reheating in the light of NANOGrav results,” [arXiv:2011.03323 [astroph.CO]].

[8] Y. Cai and Y. S. Piao, “Intermittent NEC violations during inflation and primordial gravitational waves,” [arXiv:2012.11304 [gr-qc]].

[9] J. Ellis and M. Lewicki, “Cosmic String Interpretation of NANOGrav Pulsar Timing Data,” [arXiv:2009.06555 [astro-ph.CO]].

[10] S. Blasi, V. Brdar and K. Schmitz, “Has NANOGrav found first evidence for cosmic strings?,” Phys. Rev. Lett. 126, 041305 (2021) [arXiv:2009.06607 [astro-ph.CO]].

[11] W. Büchmüller, V. Domcke and K. Schmitz, “From NANOGrav to LIGO with metastable cosmic strings,” Phys. Lett. B 811, 135914 (2020) [arXiv:2009.10649 [astro-ph.CO]].

[12] R. Samanta and S. Datta, “Gravitational wave complementarity and impact of NANOGrav data on gravitational leptogenesis: cosmic strings,” [arXiv:2009.13452 [hep-ph]].

[13] J. Liu, R. G. Cai and Z. K. Guo, “Large anisotropies of the stochastic gravitational wave background from cosmic domain walls,” [arXiv:2010.03225 [astro-ph.CO]].

[14] A. Paul, U. Mukhopadhyay and D. Majumdar, “Gravitational Wave Signatures from Domain Wall and Strong First-Order Phase Transitions in a Two Complex Scalar extension of the Standard Model,” [arXiv:2010.03439 [hep-ph]].

[15] M. Lewicki and V. Vaskonen, “Gravitational waves from colliding vacuum bubbles in gauge theories,” [arXiv:2012.07826 [astro-ph.CO]].

[16] Y. Nakai, M. Suzuki, F. Takahashi and M. Yamada, “Gravitational Waves and Dark Radiation from Dark Phase Transition: Connecting NANOGrav Pulsar Timing Data and Hubble Tension,” [arXiv:2009.09754 [astroph.CO]].

[17] A. Addazi, Y. F. Cai, Q. Gan, A. Marciano and K. Zeng, “NANOGrav results and Dark First Order Phase Transitions,” [arXiv:2009.10327 [hep-ph]].

[18] N. Kitajima, J. Soda and Y. Urakawa, “Nano-Hz gravitational wave signature from axion dark matter,” [arXiv:2010.10990 [astro-ph.CO]].

[19] N. Ramberg and L. Visinelli, “The QCD Axion and Gravitational Waves in light of NANOGrav results,” [arXiv:2012.06882 [astro-ph.CO]].

[20] V. S. H. Lee, A. Mitridate, T. Trickle and K. M. Zurek, “Probing Small-Scale Power Spectra with Pulsar Timing Arrays,” [arXiv:2012.09857 [astro-ph.CO]].

[21] M. Gorghetto, E. Hardy and H. Nicolaescu, “Observing Invisible Axions with Gravitational Waves,” [arXiv:2101.11007 [hep-ph]].

[22] A. Neronov, A. Roper Pol, C. Caprini and D. Semikoz, “NANOGrav signal from MHD turbulence at QCD phase transition in the early universe,” [arXiv:2009.14174 [astro-ph.CO]].

[23] K. T. Abe, Y. Tada and I. Ueda, “Induced gravitational waves as a cosmological probe of the sound speed during the QCD phase transition,” [arXiv:2010.06193 [astroph.CO]].

[24] J. M. Quashnock, A. Loeb and D. N. Spergel, “Magnetic Field Generation During the Cosmological QCD Phase Transition,” Astrophys. J. Lett. 344, L49 (1989)

[25] G. Sigl, A. V. Olinto and K. Jedamzik, “Primordial magnetic fields from cosmological first order phase transitions,” Phys. Rev. D 55, 4582 (1997) [arXiv:astroph/9610201 [astro-ph]].

[26] M. M. Forbes and A. R. Zhitnitsky, “Primordial galactic magnetic fields from domain walls at the QCD phase transition,” Phys. Rev. Lett. 85, 5268 (2000) [arXiv:hepph/0004051 [hep-ph]].

[27] D. Boyanovsky, H. J. de Vega and M. Simionato, “Large scale magnetogenesis from a nonequilibrium phase transition in the radiation dominated era,” Phys. Rev. D 67, 123505 (2003) [arXiv:hep-ph/0211022 [hep-ph]].

[28] L. S. Kisslinger, S. Walawalkar and M. B. Johnson, “Basic treatment of QCD phase transition bubble nucleation,” Phys. Rev. D 71, 065017 (2005) [arXiv:hepph/0503157 [hep-ph]].

[29] P. V. Buividovich, M. N. Chernodub, E. V. Luschevskaya and M. I. Polikarpov, “Numerical evidence of chiral magnetic effect in lattice gauge theory,” Phys. Rev. D 80, 054503 (2009) [arXiv:0907.0494 [hep-lat]].

[30] F. R. Urban and A. R. Zhitnitsky, “Large-Scale Magnetic Fields, Dark Energy and QCD,” Phys. Rev. D 82, 043524 (2010) [arXiv:0912.3248 [astro-ph.CO]].

[31] D. H. Rischke, “The Quark gluon plasma in equilibrium,” Prog. Part. Nucl. Phys. 52, 197 (2004) [arXiv: nuclth/0305030].

[32] J. Kiskis, R. Narayanan and H. Neuberger, “Does the crossover from perturbative to nonperturbative physics in QCD become a phase transition at infinite N?,” Phys. Lett. B 574, 65 (2003) [arXiv:hep-lat/0308033 [hep-lat]].

[33] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, “The Order of the quantum chromodynamics transition predicted by the standard model of particle physics,” Nature 443, 675 (2006) [arXiv:hep-lat/0611014].

[34] J. P. Blaizot, E. Iancu and Y. Mehtar-Tani, “Mediuminduced QCD cascade: democratic branching and wave turbulence,” Phys. Rev. Lett. 111, 052001 (2013) [arXiv:1301.6102 [hep-ph]].

[35] F. Miniati, G. Gregori, B. Reville and S. Sarkar, “Axion-Driven Cosmic Magnetogenesis during the QCD Crossover,” Phys. Rev. Lett. 121, 021301 (2018) [arXiv:1708.07614 [astro-ph.CO]].

[36] E. Witten, “Cosmic Separation of Phases,” Phys. Rev. D 30, 272 (1984).

[37] J. H. Applegate and C. J. Hogan, “Relics of Cosmic Quark Condensation,” Phys. Rev. D 31, 3037 (1985).

[38] S. Midorikawa, “Bubble Collisions in the Cosmological Quark - Hadron Phase Transition,” Phys. Lett. B 158, 107 (1985).

[39] M. Hindmarsh, “Axions and the QCD phase transition,” Phys. Rev. D 45, 1130 (1992).

[40] W. N. Cottingham, D. Kalafatis and R. Vinh Mau, “Bubble nucleation rates in first order phase transitions,” Phys. Rev. B 48, 6788(1993).

[41] J. Ignatius, K. Kajantie, H. Kurki-Suonio and M. Laine, “The growth of bubbles in cosmological phase transitions,” Phys. Rev. D 49, 3854 (1994) [arXiv:astroph/9309059 [astro-ph]].

[42] J. Ignatius, “Early stages of growth of QCD and electroweak bubbles,” [arXiv:hep-ph/9708383 [hep-ph]].

[43] A. Strumia and N. Tetradis, “Bubble nucleation rates for cosmological phase transitions,” JHEP 11, 023 (1999) [arXiv:hep-ph/9904357 [hep-ph]].

[44] A. Strumia and N. Tetradis, “A Consistent calculation of bubble nucleation rates,” Nucl. Phys. B 542, 719 (1999) [arXiv:hep-ph/9806453 [hep-ph]].

[45] D. J. Schwarz, ”The first second of the universe,” Annalen Phys. 12, 220 (2003) [arXiv:astro-ph/0303574 [astro-ph]].

[46] A. Tawfik, “The Hubble parameter in the early universe with viscous QCD matter and finite cosmological constant,” Annalen Phys. 523, 423 (2011) [arXiv:1102.2626 [gr-qc]].

[47] A. Tawfik and T. Harko, “Quark-Hadron Phase Transitions in Viscous Early Universe,” Phys. Rev. D 85, 084032 (2012) [arXiv:1108.5697 [astro-ph.CO]].

[48] M. Maggiore, ”Gravitational Waves. Vol. 1: Theory and Experiments” (Oxford University Press, 2007).

[49] E. W. Kolb and M. S. Turner, “The Early Universe,” Addison-Wesley Publishing Company (1990).

[50] L. Husdal, “On Effective Degrees of Freedom in the Early Universe,” Galaxies 4, 78 (2016) [arXiv:1609.04979 [astro-ph.CO]].

[51] S. Burke-Spolaor, S. R. Taylor, M. Charisi, T. Dolch, J. S. Hazboun, A. M. Holgado, L. Z. Kelley, T. J. W. Lazio, D. R. Madison and N. McMann, et al. “The Astrophysics of Nanohertz Gravitational Waves,” Astron. Astrophys. Rev. 27, 5 (2019) [arXiv:1811.08826 [astro-ph.HE]].

[52] A. Sesana, A. Vecchio and C. N. Colacino, “The stochastic gravitational-wave background from massive black hole binary systems: implications for observations with Pulsar Timing Arrays,” Mon. Not. Roy. Astron. Soc. 390, 192 (2008) [arXiv:0804.4476 [astro-ph]].

[53] S. A. Sanidas, R. A. Battye and B. W. Stappers, “Constraints on cosmic string tension imposed by the limit on the stochastic gravitational wave background from the European Pulsar Timing Array,” Phys. Rev. D 85, 122003 (2012) [arXiv:1201.2419 [astro-ph.CO]].

[54] C. Cutler, S. Burke-Spolaor, M. Vallisneri, J. Lazio and W. Majid, “The Gravitational-Wave Discovery Space of Pulsar Timing Arrays,” Phys. Rev. D 89, 042003 (2014) [arXiv:1309.2581 [gr-qc]].

[55] J. J. Blanco-Pillado, K. D. Olum and X. Siemens, “New limits on cosmic strings from gravitational wave observation,” Phys. Lett. B 778, 392 (2018) [arXiv:1709.02434 [astro-ph.CO]].

[56] P. Campeti, E. Komatsu, D. Poletti and C. Baccigalupi, “Measuring the spectrum of primordial gravitational waves with CMB, PTA and Laser Interferometers,” JCAP 01, 012 (2021) [arXiv:2007.04241 [astro-ph.CO]].

[57] L. P. Grishchuk, “Amplification of gravitational waves in an istropic universe,” Sov. Phys. JETP 40, 409 (1975).

[58] V. A. Rubakov, M. V. Sazhin and A. V. Veryaskin, “Graviton Creation in the Inflationary Universe and the Grand Unification Scale,” Phys. Lett. B 115, 189 (1982).

[59] A. A. Starobinsky, “Spectrum of relict gravitational radiation and the early state of the universe,” JETP Lett. 30, 682(1979).

[60] C. Caprini, R. Durrer and X. Siemens, “Detection of gravitational waves from the QCD phase transition with pulsar timing arrays,” Phys. Rev. D 82, 063511 (2010) [arXiv:1007.1218 [astro-ph.CO]].

[61] S. Capozziello, M. Khodadi and G. Lambiase, “The quark chemical potential of QCD phase transition and the stochastic background of gravitational waves,” Phys. Lett. B 789, 626 (2019) [arXiv:1808.06188 [gr-qc]].

[62] A. Kobakhidze, C. Lagger, A. Manning and J. Yue, “Gravitational waves from a supercooled electroweak phase transition and their detection with pulsar timing arrays,” Eur. Phys. J. C 77, 570 (2017) [arXiv:1703.06552 [hep-ph]].

[63] M. A. Ruderman and P. G. Sutherland, “Theory of pulsars: Polar caps, sparks, and coherent microwave radiation,” Astrophys. J. 196, 51 (1975).

[64] T. Gold, ”Rotating Neutron Stars as the Origin of the Pulsating Radio Sources”, Nature, 218, 731 (1968).

[65] M. Maggiore, ”Gravitational Waves. Vol. 2: ”Astrophysics and Cosmology” (Oxford University Press, 2017).

[66] B. Goncharov, D. J. Reardon, R. M. Shannon, X. J. Zhu, E. Thrane, M. Bailes, N. D. R. Bhat, S. Dai, G. Hobbs and M. Kerr, et al. “Identifying and mitigating noise sources in precision pulsar timing data sets,” [arXiv:2010.06109 [astro-ph.HE]].

[67] http://ipta4gw.org/

[68] G. Hobbs, A. Archibald, Z. Arzoumanian, D. Backer, M. Bailes, N. D. R. Bhat, M. Burgay, S. Burke-Spolaor, D. Champion and I. Cognard, et al. “The international pulsar timing array project: using pulsars as a gravitational wave detector,” Class. Quant. Grav. 27, 084013 (2010) [arXiv:0911.5206 [astro-ph.SR]].

[69] A. Brazier, S. Chatterjee, T. Cohen, J. M. Cordes, M. E. DeCesar, P. B. Demorest, J. S. Hazboun, M. T. Lam, R. S. Lynch and M. A. McLaughlin, et al. “The NANOGrav Program for Gravitational Waves and Fundamental Physics,” [arXiv:1908.05356 [astro-ph.IM]].

[70] A. Brandenburg, K. Enqvist and P. Olesen, “Large scale magnetic fields from hydromagnetic turbulence in the very early universe,” Phys. Rev. D 54, 1291 (1996) [arXiv:astro-ph/9602031 [astro-ph]].

[71] A. Brandenburg, T. Kahniashvili, S. Mandal, A. R. Pol, A. G. Tevzadze and T. Vachaspati, “Evolution of hydromagnetic turbulence from the electroweak phase transition,” Phys. Rev. D 96, 123528 (2017). [arXiv:1711.03804

[astro-ph.CO]].

[72] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky and S. Mandal, “The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence,” Geophys. Astrophys. Fluid Dynamics 114, 130 (2020) [arXiv:1807.05479 [physics.fludyn]].

[73] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky, “Numerical Simulations of Gravitational Waves from Early-Universe Turbulence,” Phys. Rev. D 102, 083512 (2020) [arXiv:1903.08585 [astro-ph.CO]].

[74] G. K. Batchelor and I. Proudman, “The large-scale structure of homogeneous turbulence,” Phil. Trans. Roy. Soc. Lond. A, 248, 369 (1956).

[75] P. G. Saffman, “Note on decay of homogeneous turbulence,” Phys. Fluids 10, 1349 (1967).

[76] A. Brandenburg, T. Kahniashvili, S. Mandal, A. R. Pol, A. G. Tevzadze and T. Vachaspati, “The dynamo effect in decaying helical turbulence,” Phys. Rev. Fluids. 4, 024608 (2019) [arXiv:1710.01628 [physics.flu-dyn]].

[77] D. N. Hosking and A. A. Schekochihin, “Reconnectioncontrolled decay of magnetohydrodynamic turbulence and the role of invariants,” arXiv:2012.01393 (2021).

[78] A. Brandenburg, A. Johansen, P. A. Bourdin, W. Dobler, W. Lyra, M. Rheinhardt, S. Bingert, N. E. L. Haugen, A. Mee, F. Gent, N. Babkovskaia, C.-C. Yang, T. Heinemann, B. Dintrans, D. Mitra, S. Candelaresi, J. Warnecke, P. J. Käpylä, A. Schreiber, P. Chatterjee, M. J. Käpylä, X.-Y. Li, J. Krüger, J. R. Aarnes, G. R. Sarson, J. S. Oishi, J. Schober, R. Plasson, C. Sandin, E. Karchniwy, L. F. S. Rodrigues, A. Hubbard, G. Guerrero, A. Snodin, I. R. Losada, J. Pekkilä, and C. Qian, “The Pencil Code, a modular MPI code for partial differential equations and particles: multipurpose and multiusermaintained,” Journal of Open Source Software 6, 2807 (2021).

[79] A. Brandenburg, Y. He, T. Kahniashvili, M. Rheinhardt, and J. Schober, “Gravitational waves from the chiral magnetic effect,” Astrophys. J., in press, arXiv:2101.08178 (2021).

[80] U. Frisch, A. Pouquet, J. Léorat, and A. Mazure, J. Fluid Mech. 68, 769 (1975). 778Possibility of an inverse cascade of magnetic helicity in hydrodynamic turbulence

[81] A. Brandenburg and T. Kahniashvili, Phys. Rev. Lett. 118, 055102 (2017). Classes of hydrodynamic and magnetohydrodynamic turbulent decay

[82] J. J. Blanco-Pillado, K. D. Olum and J. M. Wachter, “Comparison of cosmic string and superstring models to NANOGrav 12.5-year results,” [arXiv:2102.08194 [astroph.CO]].