Can we observe the QCD phase transition-generated gravitational waves through pulsar timing arrays?
aa r X i v : . [ a s t r o - ph . C O ] F e b NORDITA-2021-016
Can we observe the QCD phase transition-generated gravitational waves throughpulsar timing arrays?
Axel Brandenburg,
1, 2, 3, 4, ∗ Emma Clarke † , ‡ Yutong He,
1, 2, § and Tina Kahniashvili
4, 3, 5, 6, ¶ Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia McWilliams Center for Cosmology and Department of Physics,Carnegie Mellon University, Pittsburgh, PA 15213, USA Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia Department of Physics, Laurentian University, Sudbury, ON P3E 2C, Canada (Dated: February 25, 2021)We perform numerical simulations of gravitational waves (GWs) induced by hydrodynamic andhydromagnetic turbulent sources that might have been present at cosmological quantum chromody-namic (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 − of the critical energy densityin 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 ofthe frequency for nonhelical hydromagnetic turbulence. This implies more power at low frequenciesthan for the steeper spectra previously anticipated. The behavior toward higher frequencies dependsstrongly 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 ofhelicity, and somewhat less in the absence of helicity. For acoustic hydrodynamic turbulence, thesharp drop is replaced by a power law decay, albeit with a rather steep slope. Our study supportsearlier findings of a quadratic scaling of the GW energy with the magnetic energy of the turbulenceand inverse quadratic scaling with the peak frequency, which leads to larger GW energies underQCD conditions. I. INTRODUCTION
Gravitational wave (GWs) astronomy opens a newwindow to study the physical processes in the early uni-verse. Relic GWs can be sourced by violent processessuch as cosmological phase transitions and after genera-tion they propagate almost freely throughout the expan-sion of the universe that causes the dilution of their strainamplitude and frequency; for a review, see [1] and refer-ences therein. On the other hand, the detection of theserelic GWs is a challenging task due to their small ampli-tudes, the specific range of the characteristic frequencies,and astrophysical foregrounds [2]. Despite tremendousadvancements in GW detection techniques, the stochas-tic GW background of cosmological origin remained un-observed.Recently, the NANOGrav collaboration reportedstrong evidence for a stochastic GW background [3]. Inaddition to the possibility of GWs induced by astro-physical sources such as supermassive black holes, theNANOGrav data can also be understood as a possiblesignal from the early universe, such as inflationary GWs † Corresponding author; the authors are listed alphabetically. ∗ Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected] [4–8], cosmic strings and domain walls [9–14] phase tran-sitions including the supercooled phase transitions [15],dark phase transitions [16, 17] and quantum chromody-namic (QCD), with axionic string network and QCD ax-ion [18–21], and/or magnetic fields [22] and turbulence[23].In this paper we present a self-consistent study of theGWs from turbulent sources possibly present at QCDphase transitions. In contrast to Ref. [22], we considermuch larger typical length scales of the turbulent mo-tions. Also, will not limit the turbulent sources to thepresence of magnetic fields [24–30] but rather considerany turbulent source possibly presented at the QCD en-ergy scale. Even if primordial fields are not dynami-cally strong, turbulence can still develop at QCD energyscales [24, 28, 31–35]; the latent heat they release stillgives rise to pressure gradients resulting in macroscopicplasma motions. Given the very high Reynolds numberof the primordial plasma, such motions will inevitablydecay into turbulence [24, 35]. As already alluded toabove, particularly important for our work is the earlierfinding that the separation and size of nucleation bub-bles in a QCD phase transition is a sizeable fraction ofthe Hubble scale, see Ref. [36] for pioneering works andfollow-up papers, [37–47]. Furthermore, the assumptionof turbulence being driven by magnetic fields, allows usto avoid the requirement of first order QCD phase tran-sitions, [35].The paper is organized as follows. We first review ba-sic properties of relic GWs (Sec. II), discuss then theNANOGrav observations (Sec. III), present our numeri-cal approach (Sec. IV) and results (Sec. V) of our sim-ulations, before concluding in Sec. VI. Throughout thepaper, we use natural units with ¯ h = c = k B = 1. Wealso set the permeability of free space to unity, i.e., µ =1, expressing the electromagnetic quantities in Lorentz-Heaviside units. The Latin indices run i ∈ (1 , ,
3) anddefine the spatial coordinates, and the Greek indicesrun λ ∈ (0 , , , − , , , II. THE EARLY-UNIVERSE GRAVITATIONALWAVE SIGNAL
GWs correspond to the tensor mode of perturbations δg µν above the spatially flat, homogeneous, and isotropicFriedmann-Lemaˆıtre-Robertson-Walker (FLRW) back-ground, in the transverse-traceless (TT) gauge definedthrough the spatial component h phys ij with a h phys ij =Λ ijlm δg lm , where a is the scale factor at the physical time t phys . Here and below, super/subscript “phys” denotesphysical quantities.In order to eliminate the expansion-induced dilutionfrom the governing hydromagnetic equations, we userescaled quantities together with the conformal time t ,defined through dt = dt phys /a , which reduces the met-ric tensor to the Minkowski form. The background ex-pansion of the universe during the radiation-dominatedepoch is governed by the (dominant) radiation energydensity E rad = π g ( T ) T /
30, where g ( T ) is the effectivenumber of relativistic degrees of freedom at temperature T . In the epoch(s) of interest, the expansion of the uni-verse is fully governed by radiation, and the Hubble pa-rameter H ≡ a − da/dt phys = a − da/dt is given through H ( t ) = (8 πG/ E rad ( t ), where G is Newton’s gravita-tional constant and E rad ( t ) denotes the total energy den-sity of radiation (including all relativistic components).In order to connect physical and comoving variablesand to determine the scaling of physical quantities, wecompute the ratio of the scale factor today, a = a ( t )(here and below, “0” denotes the present moment), tothat at the time t ∗ (at the temperature T ∗ at whichthe source becomes active and the gravitational signalis generated) corresponding to the start of the simula-tion. We assume the adiabatic expansion of the universe,such that g S ( T ) T a ( T ) is constant, where g S ( T ) is the The TT gauge is determined by the TT projection tensor Λ ijkl = P ik P jl − P ij P kl , where the P ij is a transverse operator ( ∂ i P ij =0), defined as P ij = δ ij − ∂ i ∂j/ ∇ , where δ ij is the Kroneckerdelta, ∂ λ ≡ ∂/∂x λ denotes the partial derivative in respect of x λ coordinate, and ( ∇ defines the vector differential operatorwith the components equal to ∇ i ≡ ∂ i , i.e., ∇ is the Laplacianin 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 > g S ( T ) = g ( T ) [49]. Note, that our consideration below isvalid for any time period during the radiation dominatedepoch. However, we will be focused on the time periodaround the QCD energy scale (150 MeV). We also nor-malize the scale factor a ∗ ≡ a ( t ∗ ) = 1, which differs fromthe usual convention a = 1. Entropy conservation leadsto a a ∗ = 10 (cid:18) g S ( T ∗ )15 (cid:19) / (cid:18) T ∗
150 MeV (cid:19) , (1)where we have used T = 2 .
73 K and g S ( T ) = 3 .
91, whileat the QCD energy scale we have g S ( T ∗ ) ≈
15 [49]. Thedegrees of freedom at QCD is approximate due to uncer-tainty in the exact temperature of the QCD transitionand knowledge of the standard model (see discussions in[45, 50]). However, as ( a /a ∗ ) ∼ g S ( T ∗ ) / , small devia-tions in g S ( T ∗ ) will not significantly impact our results.The GW equation in physical time and space coordi-nates is given by (cid:0) ∂ t phys + 3 H∂ t phys − ∇ (cid:1) h phys ij = 16 πGT TT ij, phys , (2)where the TT superscript denotes the TT projection ofthe stress-energy tensor such that T TT ij, phys = Λ ijlm T phys lm .To make the connection with observations, we de-fine the characteristic strain, h c ( t ), which obeys h c ( t ) = h ( h phys ij ( x , t )) i /
2, where angle brackets denote volumeaveraging in physical space, and the physical energy den-sity E physGW ( t ) carried by the GWs is given by [48] E physGW ( t ) = 132 πG h ( ∂ t phys h phys ij ( x , t )) i . (3)It is then expressed in terms of today’s frequency f = k/ (2 πa ) that corresponds to the time Fouriertransform Q ( t ) = R ∞−∞ df Q ( f ) e − πft (and Q ( f ) =2 π R ∞−∞ dt Q ( t ) e − πft ) [48].The relic GW signal strength today is given throughthe normalized GW energy density parameter Ω GW ( f )reduced by the factor ( H ∗ /H ) ( a ∗ /a ) , where H ∗ is theHubble parameter at t ∗ . This accounts for the dilutionof the GW energy density parameter with the expansionof the universe and renormalizes the GW energy den-sity by the critical energy density at the present time, E = (3 H ) / (8 πG ), where H = 100 h km s − Mpc − ≃ . × − h s − is the present value of the Hubbleparameter. A frequency of particular interest is the fre-quency f ∗ corresponding to the Hubble horizon scale at t ∗ : f ∗ = a ∗ H ∗ a ≃ (1 . × − Hz) (cid:18) g ∗ (cid:19) / (cid:18) T ∗
150 MeV (cid:19) . (4)As discussed above, there is a variety of possiblesources of a stochastic GW background in the nHzfrequency range, accessible to Pulsar Timing Arrays(PTAs), see Section III for more details and Ref. [51]for a review and references therein, and these sourcesinclude a cosmic population of supermassive black holebinaries (SMBHBs) [51, 52], cosmic strings [53–55], infla-tionary GWs [4, 56], and phase transitions in the earlyuniverse [19, 60–62]. We also present upper limits on therelic (prior to recombination) GW background strengthbased on Big Bang nucleosynthesis (BBN) and the cosmicmicrowave background (CMB), as well as theoreticallyestimated strength and characteristic frequencies for dif-ferent sources; see Sec. V. The estimated characteristicfrequency and wave number of GWs are related to eachother through 2 πf = k , and can be expressed in terms ofthe characteristics (length and time scales) of the source.In particular, if we assume that GWs could be sourcedby bubble collisions at a phase transition, we expect thefrequency of the GWs to be related to the bubble size.We consider that the bubble length scale is the Hubblehorizon H ∗ at generation divided by the total number ofphase transition bubbles N b . Then, for the QCD phasetransitions, the frequency is given by f ∗ ≃ (1 . × − Hz) (cid:18) g ∗ (cid:19) / (cid:18) T ∗
150 MeV (cid:19)(cid:18) N b (cid:19) , (5)where we have normalized to 6 bubbles expected at theQCD phase transition [45]. III. NANOGRAV DATA
A pulsar is a highly magnetized and rapidly rotatingneutron star that emits a beam of electromagnetic radi-ation along its magnetic axis [63]. The times of arrival(TOA) of these pulses are extremely regular and can bepredicted very accurately over long times [64]. The pres-ence of a GW passing between the observer and pulsarshifts the pulse TOA proportional to the amplitude of theGW [65]. By monitoring the fluctuations in the TOA ofradio pulses from millisecond pulsars (see, e.g., Ref. [66]and references therein) international PTA missions aimto probe a stochastic GW background.The maximum sensitivity of a PTA experiment is lim-ited by the total observation time. That is, the lowestdetectable frequency is on the order of the inverse of thetime span of the data (e.g., f ∼ nHz for datasets spanning ∼
10 years) [68]. Furthermore, data sampling (i.e., pul-sars are usually observed on the order of weeks [65]) lim-its the maximum detectable frequency. The NANOGrav The quantum mechanical fluctuations during the inflationaryepoch induces GWs via parametric resonance [57–59]. The International Pulsar Timing Array (IPTA) is a consortiumof consortia, comprised of the European Pulsar Timing Array(EPTA), the North American Nanohertz Observatory for Grav-itational Waves (NANOGrav), and the Parkes Pulsar TimingArray (PPTA) [67]. . µ Hz [69].PTA measurements typically characterize a stochas-tic GW background in terms of its characteristic strainspectrum h c ( f ) fitted with a power-law dependence onfrequency [3], h c ( f ) = A CP (cid:18) ff yr (cid:19) α CP , (6)where the subscript “CP” denotes a common-spectrum(CP) process (common to the observed pulsars), the spec-tral index α CP depends on the source of the stochasticGW background, and A CP is the strain amplitude at areference frequency of f yr = 1 yr − . This choice of refer-ence frequency is arbitrary and does not affect the abilityto detect a GW signal.The energy density spectrum of the GW backgroundtoday expressed in terms of the characteristic strain spec-trum is given by [65]Ω GW ( f ) = 2 π H f h c ( f ) = Ω yrGW (cid:18) ff yr (cid:19) − γ CP , (7)where we have used Eq. (6) in the second term onthe right-hand side, γ CP = 3 − α CP and Ω yrGW ≡ π A f / (3 H ). The quantity h Ω GW ( f ) is typicallyconsidered in order to remove the uncertainty in the valueof H .The NANOGrav collaboration reports joint A CP − γ CP posterior distributions [3]. Posteriors for a common-spectrum process in the NANOGrav 12.5-year data wererecovered with four models: free-spectrum, broken powerlaw, 5-frequency power law, and 30-frequency powerlaw. The fits were performed for frequencies f ∈ [2 . × − , × − ] Hz, with the exception of the 5-frequencypower law, which was fit to the five lowest frequencybins. The four lowest frequency bins have the strongestresponse to the presence of a GW background (see Fig-ure 13 of Ref. [3]). Thus, the 5-frequency power law wasfit within the signal-dominated frequency range (approx-imately f ∈ [2 . × − , . × − ] Hz). Figure 1 ofRef. [3] shows the 1 σ and 2 σ posterior contours for theamplitude A CP and spectral slope γ CP .Fig. 1 shows the NANOGrav detection expressed interms of h Ω GW ( f ) as given by Eq. (7). The shadedregions show the 2 σ confidence contours of the A CP − γ CP parameter space in terms of f and h Ω GW ( f ) forfrequencies from 2.5–100 nHz (i.e., the NANOGrav 12.5-year 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 astrophys-ical foreground sources (white dwarfs, SMBH mergers,etc), and from relic sources related to inflation and cos-mic strings, for example, and in particular, from phase - (cid:1)(cid:2)(cid:3) - (cid:1) - (cid:4)(cid:2)(cid:3) - (cid:4) - (cid:5)(cid:6) - (cid:5)(cid:7) - (cid:1) - (cid:8) - (cid:9) (cid:1)(cid:2)(cid:3) (cid:1)(cid:2) ( (cid:1) / (cid:4)(cid:5) ) (cid:1) (cid:2)(cid:3) (cid:1)(cid:2) ( (cid:1) (cid:2) (cid:3) Ω (cid:4) (cid:5) ) FIG. 1: NANOGrav 12.5-year data set 2 σ confidence con-tours for the posteriors of a common-spectrum process (seeRef. [3] Figure 1) expressed in terms of the GW energy den-sity h Ω GW ( f ) and frequency f . This is shown over theNANOGrav 12.5-year sensitivity range of 2.5–100 nHz. Thethree models used to fit the process include a: broken powerlaw (blue), 5-frequency power law (orange), and 30-frequencypower law (green). transition-generated turbulence and primordial magneticfields. We focus here on the latter two. Turbulenceand/or magnetic fields would only be generated duringa limited amount of time before they would decay. Thedecay process itself remains highly turbulent and couldaffect GW production. Let us therefore begin with somegeneral remarks about turbulent decay. A. Gravitational Waves from Turbulent Sources
Using scaled quantities h ij = ah phys ij and T TT ij = a T TT ij, phys , together with a ( t ) ∝ t in the radiation domi-nated epoch, the GW equation takes the form (cid:0) ∂ t − ∇ (cid:1) h ij = 16 πGa T TT ij . (8)To obtain the GW equation in Fourier (wavenumber)space, we use the Fourier transforms and the polariza-tion r = (+ , × ) decomposition of the tensor metric per-turbations and stress energy tensor projected onto theTT gauge (i.e., Q ij ( k , t ) = P r =+ , × e rij (ˆ k ) Q ij ( k , t ), where e + ij (ˆ k ) and e × ij (ˆ k ) are the polarization tensors with ˆ k theunit vector, and k = a k phys is the rescaled wavenum-ber). We use the spatial Fourier transform convention: Q ( x , t ) = R d k (2 π ) e − i k · x Q ( k , t ) and Q ( k , t ) = R d x e − i k · x Q ( x , t ). Thetransverse operator P ij in the Fourier space is given P ij (ˆ k ) = δ ij − ˆ k i ˆ j k and the TT projection operator is Λ ijkl (ˆ k ) = P ik (ˆ k ) P jl (ˆ k ) − P ij (ˆ k ) P kl (ˆ k ), correspondingly. The polariza-tion tensors e + ij (ˆ k ) and e × ij (ˆ k ) can be written as e + ij ( k ) = ˆe i ˆe j − As in earlier work [72, 73], we use normalized conformaltime, ¯ t = t/t ∗ , where t ∗ = H − ∗ is our starting time,and a ∗ = 1 has been chosen. Therefore, a = ¯ t . Wealso use the scaled wave vector, ¯ k = k /H ∗ , and a scalednormalized stress, ¯ T TT+ / × = T TT+ / × / E ∗ rad . The GW equationcan then be written in the form [72, 73] (cid:16) ∂ t + ¯ k (cid:17) h + / × ( k , t ) = 6¯ t ¯ T TT+ / × ( k , t ) , (9)but from now on we omit all overbars.Throughout this paper, all numerical results will usu-ally be presented as the scaled variables introducedabove. In particular, we quote the rms strain, h rms = h h i / , where h = h + h × = ( h ij ) /
2, and like-wise for the scaled GW energy, E GW = h ˙ h i /
6, where˙ h + / × = ∂ t h + / × with ˙ h ≡ ˙ h + ˙ h × ; see Ref. [72, 73]for additional subdominant terms that are applied in thecalculations. We sometimes also quote the (frequency de-pendent) characteristic amplitude of the physical strainmeasured today, h c ( f ) = h rms /a ; see Sec. II. B. Turbulent Sources
Turbulent flows in the early universe can be modeledby 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 anexpanding universe using conformal time and comovingvariables [70, 71] with a forcing term F in the inductionequation for B ∂ ln ρ∂t = −
43 ( ∇ · u + u · ∇ ln ρ ) + 1 ρ (cid:2) u · ( J × B ) + η J (cid:3) ,∂ u ∂t = − u · ∇ u + u ∇ · u + u · ∇ ln ρ ) + 2 ρ ∇ · ( ρν S ) − ∇ ln ρ − u ρ (cid:2) u · ( J × B ) + η J (cid:3) + 34 ρ J × B ,∂ B ∂t = ∇ × ( u × B − η J + F ) , J = ∇ × B . We recall that the conformal time t is normalized to unityat the time t ∗ of magnetic field generation, ρ is in unitsof the initial value, u is in units of the speed of light,and the magnetic energy density B / S ij = ( u i,j + u j,i ) − δ ij ∇ · u are thecomponents of the rate-of-strain tensor with commas de-noting partial derivatives, J is the current density, ν isthe kinematic viscosity, and η is the magnetic diffusivity. ˆe i ˆe j and e × ij ( k ) = ˆe i ˆe j + ˆe i ˆe j , where ˆe and ˆe are unit vectorsthat are orthogonal to ˆ k and each other; see Chapter 1 (1.2) ofRef. [48]. and Ref. [72] for further detail. The electromotive force, F , is used to model magneticfield generation with F ( x , t ) = Re[ N ˜ f ( k ) exp( i k · x + iϕ )] , (10)where the wave vector k ( t ) and the phase ϕ ( t ) changerandomly from one time step to the next. This forcingfunction is therefore white noise in time and consists ofplane waves with average wavenumber k f such that | k | lies in an interval k f − δk/ ≤ | k | < k f + δk/ δk . Here, N = f /δt / is a normalization factor, where δt is the time step and f is varied to achieve a certainmagnetic field strength after a certain time, and ˜ f ( k ) =( k × e ) / [ k − ( k · e ) ] / is a nonhelical forcing function.Here, e is an arbitrary unit vector that is not alignedwith k . Note that | f | = 1. Following earlier work, theforcing is only enabled during the time interval 1 ≤ t ≤ E K ( t ) = h ρ u i / E M ( t ) = h B i /
2, respectively.The vigor of turbulence is characterized by theReynolds number, Re = u rms /νk f , where u rms is the max-imum rms velocity. It can only be determined a poste-riori from the velocity resulting from the magnetic fieldthrough 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 withtime t like ∆ t − p and the correlation length ξ M ( t ) in-creases like ∆ t q , where ∆ t = t − t off is the time interval af-ter the forcing has been turned off. The exponents p and q are positive and depend on the physical circumstances(magnetically or kinetically dominated turbulence), andwhether or not there is magnetic helicity. In helical tur-bulence, for example, one finds p = q = 2 /
3, while fornon-helical magnetically dominated turbulence one finds p = 1 and q = 1 /
2, although other variants are sometimespossible [76, 77].In this paper, we are specifically interested in the de-pendence of the decay behavior on the forcing wave num-ber k f of the turbulence while it was still being driven.The parameter k f enters through the prefactor in the de-cay law.Furthermore, ∆ t − p would become infinite for p > t = 0 (when t = t off ). The singularity of ∆ t − p at t = 0 is a consequence of a simplified description at theinitial time moment. For this reason, it is convenient toexpress the decay laws as E M ( t ) = E maxM (1 + ∆ t/τ ) − p , (11)where τ is the turnover time, which we will treat as anempirical parameter that we expect to be of the order of( v A k f ) − , where v A = (3 E maxM / / is the Alfv´en speed,evaluated at the time when E M reaches its maximumvalue E maxM . In some simulations of purely hydrodynamicturbulence, we replace E M by E K in Eq. (11) and use τ = ( u rms k f ) − with u rms = (2 E K ) / as the nominalturnover time. We recall here that we are using nondi-mensional variables where the radiation energy density isunity.For all our simulations, we choose t off = 2, i.e., turbu-lence 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, butin all other cases, it exceeds it by up to factors betweenten (in the nonhelical cases) and a hundred (in helicalcases).We arrange the strength of the forcing f such that E M is similar for different values of k f . This allows us then todetermine the resulting GW energy solely as a functionof k f . For small values of k f , the turbulence may not beable to reach a statistically steady state by the time t off ,when the driving is turned off.It is therefore necessary to adjust f for each valueof k f separately. Once we have two values of E maxM thatare close enough to the target strength, we determine thedesired forcing strength through linear interpolation. Wealso consider the case of different values of f for a fixedvalue of k f (Runs noh5,6 and Runs hel5,6). V. NUMERICAL SIMULATIONS
We solve the governing equations using the
PencilCode [78], where the GW solver has already been imple-mented [72]. We consider a cubic domain of side length2 π/k , where k is the smallest wave number in the do-main. We choose k = k f /
6, so that the scale separationbetween the initial spectral peak and the lowest wavenumber in the domain is six. In the following, we discussthe results for different values of k f . The temporal growthof 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 forthe nonhelical and helical cases. The parameters of thoseruns are listed in Tables I and II (for nonhelical and heli-cal runs). The numerical resolution is 512 mesh points,except for run noh1, where we use 1024 mesh points.Unless specified otherwise, we use ν = η = 5 × − .In Table I, we have quoted the values of E satGW and h satrms obtained at the end of the simulation. To com-pute the relic observable h Ω GW at the present time,we have to multiply E satGW by a factor ( H ∗ /H ) ( a ∗ /a ) ;see Refs. [72, 73] for details. Using g ∗ = 15 and T ∗ =150 MeV, we find H ∗ = 1 . × s − , and thus this factoris ≈ × − . The largest value of E satGW quoted in Table Iis 3 . × − and corresponds therefore to h Ω GW ≈ − .Likewise, the values of h satrms in Table I have to be mul-tiplied by a − ≈ − to obtain the observable h c atthe present time; see Eq. (1). Again, the largest value of h satrms = 5 × − corresponds therefore to the observable h c = 5 × − .To simplify the comparisons, we have arranged theforcing amplitude f such that E maxM is similar in cer-tain cases. The values of E maxM listed in the upper block 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 f p τ E maxM E satGW h satrms B [ µ G] h Ω GW ( f ) h c noh1 2 0 . . × − . . × − . × − . × − .
78 1 . × − . × − noh2 6 1 6 . × − . . . × − . × − . × − .
78 1 . × − . × − noh3 20 3 2 . × − . . . × − . × − . × − .
78 3 . × − . × − noh4 60 10 1 . × − . .
43 3 . × − . × − . × − .
79 8 . × − . × − noh5 2 0 . . × − — — 1 . × − . × − . × − .
41 8 . × − . × − noh6 2 0 . . × − — — 9 . × − . × − . × − . . × − . × − noh7 6 1 2 . × − — — 4 . × − . × − . × − .
27 2 . × − . × − noh8 6 1 1 . × − — — 8 . × − . × − . × − . . × − . × − TABLE II: Similar to Table I, but for helical turbulence.Run k f k f p τ E maxM E satGW h satrms B [ µ G] h Ω GW ( f ) h c hel1 2 0 . . × − .
67 100 3 . × − . × − . × − .
79 1 . × − . × − hel2 6 1 5 . × − .
67 20 3 . × − . × − . × − .
78 1 . × − . × − hel3 20 3 2 . × − .
67 4 . . × − . × − . × − .
80 2 . × − . × − hel4 60 10 6 . × − .
67 0 .
50 3 . × − . × − . × − .
78 2 . × − . × − hel5 2 0 . . × − — — 1 . × − . × − . × − .
41 1 . × − . × − hel6 2 0 . . × − — — 9 . × − . × − . × − . . × − . × − hel7 6 1 2 . × − — — 4 . × − . × − . × − .
28 2 . × − . × − hel8 6 1 1 . × − — — 1 . × − . × − . × − . . × − . × − of Tables I and II (for nonhelical and helical hydromag- netic turbulence, respectively) are around 0.038 and cor- FIG. 3: Similar to Fig. 2(a) and (b), but in a double-logarithmic representation where E M is now plotted versus ∆ t ≡ t − f ν E maxM E satGW h satrms B [ µ G] h Ω GW ( f ) h c magnetic 1 . × − . × − . × − . × − . × − .
78 1 . × − . × − vortical 3 . × − . × − . × − . × − . × − .
82 2 . × − . × − irrotational 7 . × − . × − . × − . × − . × − .
83 2 . × − . × − respond to 0 . µ 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 satGW some time af-ter E M ( t ) has reached its maximum, and is smaller forlarger values of k f .It is important to realize that in all four cases, the de-cay of the magnetic energy follows an approximate powerlaw decay, as given by Eq. (11). To see this, we the plotin Fig. 3 the evolution of E M versus t − p and τ de-scribe the decay and are also listed in Table I.Our results confirm that the turbulence decays moreslowly for large values of τ , or small values of k f . As al-ready found from earlier simulations [73], the GW energygenerally decreases with increasing k f . This is seen moreclearly in a diagram of E GW versus E M /k f ; see Fig. 4(a).For k f = 2, we have performed additional simulationswith smaller and larger values f , both with and with-out helicity. The resulting values of E satGW obey quadraticscaling of the form E satGW = ( q E maxM /k f ) (12)with a coefficient q = 1 .
1; see the straight line inFig. 4(a). Only the data point for k f = 60 is slightlyabove the line represented by Eq. (12). This could bean artefact of our Reynolds numbers still not being largeenough in our simulations, especially for large value of k f .To compare with earlier work, we show in Fig. 4(b) thepositions of our runs in a E satGW versus E maxM diagram. Fororientation, we also show the data points from Ref. [73]. We see that the new data points are well above the olderones of Ref. [73]. This is mainly a consequence of us-ing here smaller values of k f (2–60, compared to 600 inRef. [73]). For k f = 2 and k = 0 .
3, we show here theresults for hydrodynamic runs using irrotational and vor-tical forcings; see the green symbols in Fig. 4(b). Thoseruns are listed in Table III and compared with the non-helical magnetic turbulence run ‘noh1’.In Fig. 5, we plot the resulting present-day GW en-ergy and strain spectra for our four runs with k f = 2, 6,20, and 60, both without and with helicity in the driv-ing function F . The first two cases with k f = 2 and6 lie well within the frequency and amplitude range ac-cessible to NANOGrav. In all cases, the spectra show asharp drop slightly above the peak frequency. This is aconsequence of the rapid temporal growth of the spectra,which leads to a correspondingly large growth at the peakfrequency, while at higher frequencies, the spectrum set-tled at values that were determined by somewhat earliertimes when the energy was still weaker.At frequencies below the peak, we now find a spec-trum that is even shallower than the h Ω GW ( f ) ∝ f spectrum found already earlier [73]. A spectrum shal-lower than proportional to f , such as the present f / spectrum, could perhaps be explained by the finite sizeof the computational domain; see Ref. [79], who foundeven a f − / spectrum for Ω GW ( f ). Alternatively, theshallower spectrum might well be physical, or at leastsignificantly extended over a substantial frequency inter-val below the peak frequency, for example due to inversecascading in helical [80] and nonhelical [81] cases. FIG. 4: (a) E GW versus E M /k f ; the straight line shows E GW = 5 . × − ( E M /k f ) / . (b) Positions of our runs in a diagramshowing E satGW versus E maxM . For orientation the old data points of the Ref. [73] are shown as gray symbols. The open red (filledblue) symbols are for the helical (nonhelical) runs. The green symbols refer to the two hydromagnetic runs of Table III.FIG. 5: h Ω 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. In the absence of sources, a Ω GW ( f ) ∝ f α spectrumimplies h c ( f ) ∝ f α/ − for arbitrary spectral indices α .For α = 1 /
2, we would thus expect h c ( f ) ∝ f − / . How-ever, the observed strain spectrum, h c ( f ) ∝ f − / , seemsto agree with that found previously from numerical simu-lations [73]. However, looking more carefully at the strainspectrum for k f = 60, we see a h c ( f ) ∝ f − / spectrumis actually compatible with the simulation; see the corre-sponding dashed-dotted line in Fig. 5(c). This agreement is probably related to the fact that the turnover time isshorter for the run with k f = 60, compared with those atsmaller values (i.e., longer turbulence driving time willallow for more efficient inverse cascading).In the runs with helicity, we do find h Ω GW ( f ) ∝ f ,together with a slight enhancement just before reachinga maximum. The subsequent decay for larger values of f is much steeper in the case with helicity than with-out. Furthermore, in h c ( f ) we see a sharper drop to the FIG. 6: Magnetic energy spectra for the (a) nonhelical and (b) helical cases.FIG. 7: Similar to Fig. 5, but comparing vortical (red) andirrotational turbulence (blue) with MHD turbulence (black). right of the maximum than in simulations without he-licity. These differences in the spectra for helical andnonhelical cases are surprisingly strong and might allowus to infer the presence of magnetic helicity once such aspectrum is detected.It is important to note that the h Ω GW ∝ f spectrain Figs. 4(a) and (c) show an increase towards smaller k f . This is to be expected from Eq. (12), but it was notincluded in the sketch of Ref. [22]; see their Fig. (1). Bycontrast, in their Eq. (4), an effectively cubic dependenceon the magnetic energy was motivated.The underlying magnetic energy spectrum is shown inFig. 6(a) for nonhelical and in Fig. 6(b) helical caseswhere k f is ranging from 2 to 60. Those are averagedspectra obtained by averaging over the time interval15 ≤ t ≤
20. In the nonhelical case, the amplitude of thespectrum is smaller for larger values of k f , because herethe energy has decayed more rapidly. In the helical case,the spectra have approximately the same height for allvalues of k f . This is because the height of the spectrumis related to the helicity, which is conserved. For smallvalues of k f , the spectrum has a more extended subin-ertial range. This is because the turnover time is largerand there was not enough time for the inverse cascade toproduce energy and small values of k .Finally, we compare the results for two types of purelyhydrodynamic turbulence with vortical and irrotationalforcings of Table III. The result is shown in Fig. 7. Allthese cases are for plane wave forcings. For irrotationalforcing, we do not see the sharp drop-off of spectral powerfor frequencies above the peak value as in the vorticalcase. This suggests that in the inertial range of irrota-tional turbulence, there is still some power to contributeto GW driving compared with the vortical case, wherethis is almost not possible at all. However, the spectrumin the irrotational case shows a fairly steep spectrum pro-portional to f − , so the effect on GW production is herealso rather weak. Nevertheless, the spectral form of thepeak might give interesting diagnostic clues about thenature of turbulent driving at the time of GW produc-tion.We mention in passing that in earlier work, it wasfound that irrotational turbulence is much more efficientin driving GWs than vortical turbulence [73]. Remark-ably, here this is no longer the case and vortical andirrotational turbulence have rather similar GW energies.0 - (cid:1)(cid:2)(cid:3) - (cid:1) - (cid:4)(cid:2)(cid:3) - (cid:4) - (cid:5)(cid:2)(cid:3) - (cid:6)(cid:6) - (cid:6)(cid:7) - (cid:8)(cid:1) - (cid:8)(cid:5) - (cid:8)(cid:9) - (cid:8)(cid:6) - (cid:8)(cid:7) - (cid:1) - (cid:5) - (cid:9) (cid:10)(cid:11)(cid:12) (cid:1)(cid:2) ( (cid:1) / (cid:13)(cid:14) ) (cid:1) (cid:2)(cid:3) (cid:1)(cid:2) ( (cid:1) (cid:2)(cid:3) Ω (cid:4) (cid:5) ) (cid:1)(cid:2)(cid:3)(cid:3)(cid:3)(cid:4) - (cid:1)(cid:2)(cid:3) - (cid:1) - (cid:4)(cid:2)(cid:3) - (cid:4) - (cid:5)(cid:2)(cid:3) - (cid:8)(cid:6) - (cid:8)(cid:7) - (cid:1) - (cid:5) - (cid:9) (cid:10)(cid:11)(cid:12) (cid:1)(cid:2) ( (cid:1) / (cid:13)(cid:14) ) (cid:1) (cid:2)(cid:3) (cid:1)(cid:2) ( (cid:1) (cid:2)(cid:3) Ω (cid:4) (cid:5) ) - (cid:1)(cid:2)(cid:3) - (cid:1) - (cid:4)(cid:2)(cid:3) - (cid:4) - (cid:5)(cid:2)(cid:3) - (cid:8)(cid:6) - (cid:8)(cid:7) - (cid:1) - (cid:5) - (cid:9) (cid:10)(cid:11)(cid:12) (cid:1)(cid:2) ( (cid:1) / (cid:13)(cid:14) ) (cid:1) (cid:2)(cid:3) (cid:1)(cid:2) ( (cid:1) (cid:2)(cid:3) Ω (cid:4) (cid:5) ) - (cid:1)(cid:2)(cid:3) - (cid:1) - (cid:4)(cid:2)(cid:3) - (cid:4) - (cid:5)(cid:2)(cid:3) - (cid:8)(cid:6) - (cid:8)(cid:7) - (cid:1) - (cid:5) - (cid:9) (cid:10)(cid:11)(cid:12) (cid:1)(cid:2) ( (cid:1) / (cid:13)(cid:14) ) (cid:1) (cid:2)(cid:3) (cid:1)(cid:2) ( (cid:1) (cid:2)(cid:3) Ω (cid:4) (cid:5) ) FIG. 8: Upper left: Runs noh5,6 (yellow) and hel5,6 (magenta), which correspond to k f = 2, and noh7,8 (purple), correspondingto k f = 6, compared with the NANOGrav 12.5-year 2 σ contours (as shown in Fig. 1). Additionally, the gray and black dashedhorizontal lines show the CMB and BBN integrated bounds on h Ω GW ( f ) (see [65] for details). Upper right: Contoursrepresenting four different cosmic string average power spectrum models with different tensions, as described in [82]: mono(orange), kink (green), cusp (blue), and a spectrum computed from a simulated gravitational backreaction model (magenta).Lower left: Contour representing the n T − r (tensor spectral index and tensor-to-scalar ratio, respectively) parameter spaceconsistent with the NANOGrav 12.5-year 5-frequency power-law 2 σ confidence contour considering the Bicep2-Keck Array and
Planck constraint that r < .
07 [4]. Lower right: NANOGrav 12.5-year contours for SMBHBs, which are expected to have aspectral index of γ CP = 13 /
3, corresponding to A CP = [1 . , . × − from the NANOGrav 12.5-year 5-frequency and brokenpower law 2 σ contours. This could be related to the small value of k f , possi-bly combined with a comparatively short time of driving.However, to clarify this further, more targeted numericalexperiments would need to be performed.To put our results into perspective, we compare inFig. 8 with the contours for possible sources of GWs inthe nHz range in terms of h Ω GW and f . In the upperleft, we show limits on GWs from MHD turbulence atQCD for k f from 2 to 6 (corresponding to Runs noh5,6,hel5,6, and noh7,8) around the NANOGrav sensitivityrange. Additionally, we show the integrated bounds on h Ω GW from the CMB and BBN [65], noting that theactual bound on the peak of h Ω GW can fall above theselines. Contours in the upper right correspond to four dif-ferent models for the average power spectrum of GWsfrom a network of cosmic strings of different tensions[82]. The bottom left contour corresponds to the param-eter space of n T − r (tensor spectral index and tensor-to-scalar ratio) consistent with the 2 σ contours of theNANOGrav 12.5-year 5-frequency power law. This cor-responds to Figs. 1 and 2 of Ref. [4]. The lower rightcontour represents the NANOGrav 12.5-year 2 σ contours for the 5-frequency and broken power laws for a popu-lation of SMBHBs, expected to have a spectral index of γ CP = 13 /
3. In this work we found Ω GW ( f ) ∼ f / fornonhelical turbulence, corresponding to a spectral indexof γ CP = 4 .
5, which falls at the edge of the 1 σ confidencecontours for the NANOGrav 5-frequency and brokenpower laws. The scaling Ω GW ( f ) ∼ f , which we foundin the case of helical turbulence, corresponds to γ CP = 4and falls within the 2 σ contours from NANOGrav. Thespectral index for SMBHBs, γ CP = 13 / VI. CONCLUSIONS
In the present work, we have shown that the mag-netic stress from hydrodynamic and MHD turbulencewith scales comparable to the cosmological horizon scaleat the time of the QCD phase transition can drive GWsin the range accessible to NANOGrav, if the magneticenergy density is 3–10% of the radiation energy den-sity. The low-frequency tail below the peak frequency1at 10 nHz or so is shallower in the nonhelical case thanin the helical one, i.e., ∝ f / compared to ∝ f . Bothscalings are, however, shallower than what was expectedbased on earlier analytical calculations. Also the iner-tial range spectrum above the peak is shallower withouthelicity than with, but here, both spectra are steeperthan what is expected if the GW spectrum was a directconsequence of the MHD turbulence spectrum [73, 83].The reason for this is primarily the relatively short timeof turbulent driving (one Hubble time). It is short com-pared with the turnover time which, for our runs with thesmallest k f of two, is much longer: 16 (100) Hubble timesfor our runs without (with) helicity. Therefore, there wasnot enough time to fully establish the GW spectrum athigh wave numbers. For our earlier runs with larger val-ues of k f , this effect was less pronounced than for smallervalues of k f , but it is still quite noticeable, especially inthe helical case where forward cascading is weaker thanin the nonhelical case.Our work has led to new insights regarding the pos-sibility of using an observed GW spectrum for makingstatements about the nature of the underlying turbu-lence in the early universe. One is the already mentionedslope of the subinertial range spectrum. Another is theposition of the peak of the spectrum. Finally, there isthe strength of the drop of the spectral power for fre- quencies above the peak frequency, and the subsequentslope after the drop, which is most likely too small to bedetectable. This, however, depends on the duration ofturbulent driving and could be higher if the driving timewas longer. The specific features of the spectrum nearthe peak are different for helical and nonhelical turbu-lence. This could, in principle, give information aboutthe presence of parity violation, when would also lead tocircularly polarized GWs.Data availability—The source code used for the simu-lations of this study, the Pencil Code , is freely avail-able from Ref. [78]. The simulation setups and the cor-responding data are freely available from Ref. [84].
Acknowledgments
Support through the Swedish Research Council, grant2019-04234, and Shota Rustaveli GNSF (grant FR/18-1462) are gratefully acknowledged. We acknowledgethe allocation of computing resources provided by theSwedish National Allocations Committee at the Centerfor Parallel Computers at the Royal Institute of Technol-ogy in Stockholm. [1] C. Caprini and D. G. Figueroa, “Cosmological Back-grounds of Gravitational Waves,” Class. Quant. Grav. , 163001 (2018) [arXiv:1801.04268 [astro-ph.CO]].[2] J. D. Romano and N. J. Cornish, “Detection meth-ods for stochastic gravitational-wave backgrounds: aunified treatment,” Living Rev. Rel. , 2 (2017)[arXiv:1608.06889 [gr-qc]].[3] Z. Arzoumanian et al. [NANOGrav], “The NANOGrav12.5 yr Data Set: Search for an Isotropic Stochas-tic Gravitational-wave Background,” Astrophys. J. Lett. , L34 (2020) [arXiv:2009.04496 [astro-ph.HE]].[4] S. Vagnozzi, “Implications of the NANOGrav results forinflation,” Mon. Not. Roy. Astron. Soc. , L11 (2021)doi:10.1093/mnrasl/slaa203 [arXiv:2009.13432 [astro-ph.CO]].[5] Z. Zhou, J. Jiang, Y. F. Cai, M. Sasaki and S. Pi, “Pri-mordial black holes and gravitational waves from reso-nant amplification during inflation,” Phys. Rev. D ,103527 (2020) [arXiv:2010.03537 [astro-ph.CO]].[6] H. W. H. Tahara and T. Kobayashi, “Nanohertz gravita-tional waves from NEC violation in the early universe,”Phys. Rev. D , 123533 (2020) [arXiv:2011.01605 [gr-qc]].[7] S. Kuroyanagi, T. Takahashi and S. Yokoyama, “Blue-tilted inflationary tensor spectrum and reheating in thelight of NANOGrav results,” JCAP , 071 (2021)[arXiv:2011.03323 [astro-ph.CO]].[8] Y. Cai and Y. S. Piao, “Intermittent NEC violationsduring inflation and primordial gravitational waves,”[arXiv:2012.11304 [gr-qc]].[9] J. Ellis and M. Lewicki, “Cosmic String Interpretation of NANOGrav Pulsar Timing Data,” Phys. Rev. Lett. ,041304 (2021) [arXiv:2009.06555 [astro-ph.CO]].[10] S. Blasi, V. Brdar and K. Schmitz, “Has NANOGravfound first evidence for cosmic strings?,” Phys. Rev. Lett. , 041305 (2021) [arXiv:2009.06607 [astro-ph.CO]].[11] W. Buchm¨uller, V. Domcke and K. Schmitz, “FromNANOGrav to LIGO with metastable cosmic strings,”Phys. Lett. B , 135914 (2020) [arXiv:2009.10649[astro-ph.CO]].[12] R. Samanta and S. Datta, “Gravitational wave comple-mentarity and impact of NANOGrav data on gravita-tional leptogenesis: cosmic strings,” [arXiv:2009.13452[hep-ph]].[13] J. Liu, R. G. Cai and Z. K. Guo, “Large anisotropies ofthe stochastic gravitational wave background from cos-mic domain walls,” [arXiv:2010.03225 [astro-ph.CO]].[14] A. Paul, U. Mukhopadhyay and D. Majumdar, “Gravi-tational Wave Signatures from Domain Wall and StrongFirst-Order Phase Transitions in a Two Complex Scalarextension of the Standard Model,” [arXiv:2010.03439[hep-ph]].[15] M. Lewicki and V. Vaskonen, “Gravitational wavesfrom 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 DarkPhase Transition: Connecting NANOGrav Pulsar Tim-ing Data and Hubble Tension,” [arXiv:2009.09754 [astro-ph.CO]].[17] A. Addazi, Y. F. Cai, Q. Gan, A. Marciano and K. Zeng,“NANOGrav results and Dark First Order Phase Tran- sitions,” [arXiv:2009.10327 [hep-ph]].[18] N. Kitajima, J. Soda and Y. Urakawa, “Nano-Hzgravitational wave signature from axion dark matter,”[arXiv:2010.10990 [astro-ph.CO]].[19] N. Ramberg and L. Visinelli, “The QCD Axion andGravitational 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 TimingArrays,” [arXiv:2012.09857 [astro-ph.CO]].[21] M. Gorghetto, E. Hardy and H. Nicolaescu, “Ob-serving 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 phasetransition in the early universe,” Phys. Rev. D ,L041302 (2021) [arXiv:2009.14174 [astro-ph.CO]].[23] K. T. Abe, Y. Tada and I. Ueda, “Induced gravitationalwaves as a cosmological probe of the sound speed dur-ing the QCD phase transition,” [arXiv:2010.06193 [astro-ph.CO]].[24] J. M. Quashnock, A. Loeb and D. N. Spergel, “MagneticField Generation During the Cosmological QCD PhaseTransition,” Astrophys. J. Lett. , L49 (1989).[25] G. Sigl, A. V. Olinto and K. Jedamzik, “Primor-dial magnetic fields from cosmological first orderphase transitions,” Phys. Rev. D , 4582 (1997)[arXiv:astro-ph/9610201 [astro-ph]].[26] M. M. Forbes and A. R. Zhitnitsky, “Primordial galac-tic magnetic fields from domain walls at the QCDphase transition,” Phys. Rev. Lett. , 5268 (2000)[arXiv:hep-ph/0004051 [hep-ph]].[27] D. Boyanovsky, H. J. de Vega and M. Simionato, “Largescale magnetogenesis from a nonequilibrium phase tran-sition in the radiation dominated era,” Phys. Rev. D ,123505 (2003) [arXiv:hep-ph/0211022 [hep-ph]].[28] L. S. Kisslinger, S. Walawalkar and M. B. John-son, “Basic treatment of QCD phase transition bub-ble nucleation,” Phys. Rev. D , 065017 (2005)[arXiv:hep-ph/0503157 [hep-ph]].[29] P. V. Buividovich, M. N. Chernodub, E. V. Luschevskayaand M. I. Polikarpov, “Numerical evidence of chiral mag-netic effect in lattice gauge theory,” Phys. Rev. D ,054503 (2009) [arXiv:0907.0494 [hep-lat]].[30] F. R. Urban and A. R. Zhitnitsky, “Large-Scale MagneticFields, Dark Energy and QCD,” Phys. Rev. D , 043524(2010) [arXiv:0912.3248 [astro-ph.CO]].[31] D. H. Rischke, “The Quark gluon plasma in equilib-rium,” Prog. Part. Nucl. Phys. , 197 (2004) [arXiv:nucl-th/0305030].[32] J. Kiskis, R. Narayanan and H. Neuberger, “Does thecrossover from perturbative to nonperturbative physicsin QCD become a phase transition at infinite N?,” Phys.Lett. B , 65 (2003) [arXiv:hep-lat/0308033 [hep-lat]].[33] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Sz-abo, “The Order of the quantum chromodynamics transi-tion predicted by the standard model of particle physics,”Nature , 675 (2006) [arXiv:hep-lat/0611014].[34] J. P. Blaizot, E. Iancu and Y. Mehtar-Tani, “Medium-induced QCD cascade: democratic branching and waveturbulence,” Phys. Rev. Lett. , 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. , 021301 (2018)[arXiv:1708.07614 [astro-ph.CO]].[36] E. Witten, “Cosmic Separation of Phases,” Phys. Rev. D , 272 (1984).[37] J. H. Applegate and C. J. Hogan, “Relics of CosmicQuark Condensation,” Phys. Rev. D , 3037 (1985).[38] S. Midorikawa, “Bubble Collisions in the CosmologicalQuark - Hadron Phase Transition,” Phys. Lett. B ,107 (1985).[39] M. Hindmarsh, “Axions and the QCD phase transition,”Phys. Rev. D , 1130 (1992).[40] W. N. Cottingham, D. Kalafatis and R. Vinh Mau, “Bub-ble nucleation rates in first order phase transitions,”Phys. Rev. B , 6788(1993).[41] J. Ignatius, K. Kajantie, H. Kurki-Suonio andM. Laine, “The growth of bubbles in cosmologicalphase transitions,” Phys. Rev. D , 3854 (1994)[arXiv:astro-ph/9309059 [astro-ph]].[42] J. Ignatius, “Early stages of growth of QCD and elec-troweak bubbles,” [arXiv:hep-ph/9708383 [hep-ph]].[43] A. Strumia and N. Tetradis, “Bubble nucleation rates forcosmological phase transitions,” JHEP , 023 (1999)[arXiv:hep-ph/9904357 [hep-ph]].[44] A. Strumia and N. Tetradis, “A Consistent calculation ofbubble nucleation rates,” Nucl. Phys. B , 719 (1999)[arXiv:hep-ph/9806453 [hep-ph]].[45] D. J. Schwarz, ”The first second of the universe,”Annalen Phys. , 220 (2003) [arXiv:astro-ph/0303574[astro-ph]].[46] A. Tawfik, “The Hubble parameter in the early universewith viscous QCD matter and finite cosmological con-stant,” Annalen Phys. , 423 (2011) [arXiv:1102.2626[gr-qc]].[47] A. Tawfik and T. Harko, “Quark-Hadron Phase Tran-sitions in Viscous Early Universe,” Phys. Rev. D ,084032 (2012) [arXiv:1108.5697 [astro-ph.CO]].[48] M. Maggiore, ”Gravitational Waves. Vol. 1: Theory andExperiments” (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 theEarly Universe,” Galaxies , 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. , 5 (2019) [arXiv:1811.08826[astro-ph.HE]].[52] A. Sesana, A. Vecchio and C. N. Colacino, “The stochas-tic gravitational-wave background from massive blackhole binary systems: implications for observations withPulsar Timing Arrays,” Mon. Not. Roy. Astron. Soc. , 192 (2008) [arXiv:0804.4476 [astro-ph]].[53] S. A. Sanidas, R. A. Battye and B. W. Stappers, “Con-straints on cosmic string tension imposed by the limiton the stochastic gravitational wave background fromthe European Pulsar Timing Array,” Phys. Rev. D ,122003 (2012) [arXiv:1201.2419 [astro-ph.CO]].[54] C. Cutler, S. Burke-Spolaor, M. Vallisneri, J. Lazio andW. Majid, “The Gravitational-Wave Discovery Space ofPulsar Timing Arrays,” Phys. Rev. D , 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 observa-tion,” Phys. Lett. B , 392 (2018) [arXiv:1709.02434[astro-ph.CO]].[56] P. Campeti, E. Komatsu, D. Poletti and C. Bacci-galupi, “Measuring the spectrum of primordial gravita-tional waves with CMB, PTA and Laser Interferometers,”JCAP , 012 (2021) [arXiv:2007.04241 [astro-ph.CO]].[57] L. P. Grishchuk, “Amplification of gravitational waves inan istropic universe,” Sov. Phys. JETP , 409 (1975).[58] V. A. Rubakov, M. V. Sazhin and A. V. Veryaskin,“Graviton Creation in the Inflationary Universe and theGrand Unification Scale,” Phys. Lett. B , 189 (1982).[59] A. A. Starobinsky, “Spectrum of relict gravitational ra-diation and the early state of the universe,” JETP Lett. , 682(1979).[60] C. Caprini, R. Durrer and X. Siemens, “Detection ofgravitational waves from the QCD phase transition withpulsar timing arrays,” Phys. Rev. D , 063511 (2010)[arXiv:1007.1218 [astro-ph.CO]].[61] S. Capozziello, M. Khodadi and G. Lambiase, “Thequark chemical potential of QCD phase transition andthe stochastic background of gravitational waves,” Phys.Lett. B , 626 (2019) [arXiv:1808.06188 [gr-qc]].[62] A. Kobakhidze, C. Lagger, A. Manning and J. Yue,“Gravitational waves from a supercooled electroweakphase transition and their detection with pulsar timingarrays,” Eur. Phys. J. C , 570 (2017) [arXiv:1703.06552[hep-ph]].[63] M. A. Ruderman and P. G. Sutherland, “Theory of pul-sars: Polar caps, sparks, and coherent microwave radia-tion,” Astrophys. J. , 51 (1975).[64] T. Gold, ”Rotating Neutron Stars as the Origin of thePulsating Radio Sources”, Nature, , 731 (1968).[65] M. Maggiore, ”Gravitational Waves. Vol. 2: ”Astro-physics 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 mitigat-ing 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 internationalpulsar timing array project: using pulsars as a gravita-tional wave detector,” Class. Quant. Grav. , 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 andFundamental Physics,” [arXiv:1908.05356 [astro-ph.IM]].[70] A. Brandenburg, K. Enqvist and P. Olesen, “Large scalemagnetic fields from hydromagnetic turbulence in thevery early universe,” Phys. Rev. D , 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 hydro-magnetic turbulence from the electroweak phase transi-tion,” Phys. Rev. D , 123528 (2017). [arXiv:1711.03804[astro-ph.CO]].[72] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky and S. Mandal, “The timestep constraint insolving the gravitational wave equations sourced by hy-dromagnetic turbulence,” Geophys. Astrophys. Fluid Dy-namics , 130 (2020) [arXiv:1807.05479 [physics.flu-dyn]].[73] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kah-niashvili and A. Kosowsky, “Numerical Simulations ofGravitational Waves from Early-Universe Turbulence,”Phys. Rev. D , 083512 (2020) [arXiv:1903.08585[astro-ph.CO]].[74] G. K. Batchelor and I. Proudman, “The large-scale struc-ture of homogeneous turbulence,” Phil. Trans. Roy. Soc.Lond. A, , 369 (1956).[75] P. G. Saffman, “Note on decay of homogeneous turbu-lence,” Phys. Fluids , 1349 (1967).[76] A. Brandenburg, T. Kahniashvili, S. Mandal, A. R. Pol,A. G. Tevzadze and T. Vachaspati, “The dynamo effectin decaying helical turbulence,” Phys. Rev. Fluids. ,024608 (2019) [arXiv:1710.01628 [physics.flu-dyn]].[77] D. N. Hosking and A. A. Schekochihin, “Reconnection-controlled decay of magnetohydrodynamic turbulenceand 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. Heine-mann, B. Dintrans, D. Mitra, S. Candelaresi, J. War-necke, P. J. K¨apyl¨a, A. Schreiber, P. Chatterjee, M. J.K¨apyl¨a, X.-Y. Li, J. Kr¨uger, J. R. Aarnes, G. R. Sar-son, 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¨a, and C. Qian, “ThePencil Code, a modular MPI code for partial differen-tial equations and particles: multipurpose and multiuser-maintained,” Journal of Open Source Software , 2807(2021). [arXiv:2009.08231 [astro-ph.IM]].[79] A. Brandenburg, Y. He, T. Kahniashvili, M. Rhein-hardt, and J. Schober, “Gravitational waves fromthe chiral magnetic effect,” Astrophys. J., in press,arXiv:2101.08178 (2021).[80] U. Frisch, A. Pouquet, J. L´eorat, and A. Mazure, “Pos-sibility of an inverse cascade of magnetic helicity in hy-drodynamic turbulence,” J. Fluid Mech. , 769 (1975).[81] A. Brandenburg and T. Kahniashvili, “Classes of hydro-dynamic and magnetohydrodynamic turbulent decay,”Phys. Rev. Lett. , 055102 (2017). [arXiv:1607.01360[physics.flu-dyn]].[82] J. J. Blanco-Pillado, K. D. Olum and J. M. Wachter,“Comparison of cosmic string and superstring models toNANOGrav 12.5-year results,” [arXiv:2102.08194 [astro-ph.CO]].[83] A. Brandenburg and S. Boldyrev, “The turbulent stressspectrum in the inertial and subinertial ranges,” As-trophys. J. , 80 (2020). [arXiv:1912.07499 [astro-ph.CO]].[84] A. Brandenburg, E. Clarke, Y. He, T. Kahni-ashvili, Datasets for “Can we observe the QCDphase transition-generated gravitational wavesthrough pulsar timing arrays?” (v2021.02.24), https://doi.org/10.5281/zenodo.4560423 ; see also