Circular Polarization of Gravitational Waves from Early-Universe Helical Turbulence

We perform direct numerical simulations to compute the net circular polarization of gravitational waves from helical (chiral) turbulent sources in the early universe for a variety of initial conditions, including driven (stationary) and decaying turbulence. We investigate the resulting gravitational wave signal assuming different turbulence geneses such as magnetically or kinetically driven cases. Under realistic physical conditions in the early universe we compute numerically for the first time the total (integrated over all wavenumbers) polarization degree of the gravitational waves. We find that the spectral polarization degree strongly depends on the initial conditions. The peak of the spectral polarization degree occurs at twice the typical wavenumber of the source, as expected, and for fully helical decaying turbulence, it reaches its maximum (100%) only at the peak. We determine the temporal evolution of the turbulent sources as well as the resulting gravitational waves, showing that the dominant contribution to their spectral energy density happens shortly after the source activation. Only through an artificially prolonged decay of the turbulence can further increase of the gravitational wave amplitude be achieved. We estimate the detection prospects for the net polarization arguing that its detection contains clean information (including the generation mechanisms, time and strength) about the sources of possible parity violations in the early universe.

We perform direct numerical simulations to compute the net circular polarization of gravitational waves from helical (chiral) turbulent sources in the early universe for a variety of initial conditions, including driven (stationary) and decaying turbulence. We investigate the resulting gravitational wave signal assuming different turbulence geneses such as magnetically or kinetically driven cases. Under realistic physical conditions in the early universe we compute numerically for the first time the total (integrated over all wavenumbers) polarization degree of the gravitational waves. We find that the spectral polarization degree strongly depends on the initial conditions. The peak of the spectral polarization degree occurs at twice the typical wavenumber of the source, as expected, and for fully helical decaying turbulence, it reaches its maximum (100%) only at the peak. We determine the temporal evolution of the turbulent sources as well as the resulting gravitational waves, showing that the dominant contribution to their spectral energy density happens shortly after the source activation. Only through an artificially prolonged decay of the turbulence can further increase of the gravitational wave amplitude be achieved. We estimate the detection prospects for the net polarization arguing that its detection contains clean information (including the generation mechanisms, time and strength) about the sources of possible parity violations in the early universe.
PACS numbers: 98.70.Vc, 98.80.-k A remarkable possible source of stochastic gravitational waves (GWs) from the early universe is turbulence in the primordial plasma induced either from cosmological first-order (electroweak or QCD) phase transitions [1][2][3], or from the primordial magnetic fields that are coupled to the cosmological plasma [4][5][6][7]. The GW signal is potentially detectable by the Laser Interferometer Space Antenna (LISA) in the case of strong enough turbulent motions present at the electroweak scale (assuming that the total energy in turbulence is up to 1-10% of the total thermal energy at the moment of generation) [8][9][10][11][12]. Interestingly, if the GWs are produced during the QCD phase transitions, they are potentially detectable through Pulsar Timing Arrays (PTAs) and this has been proposed to explain the reported NANOGrav signal [13]. Since GWs propagate almost freely from the moment of generation until today 1 , the detection of GWs sourced by primordial turbulence will open a new window to understand physical processes in the very early stages of the evolution of the universe (at the femtoseconds timescale); see Ref. [17] and references therein. Moreover, several theoretical extensions of the standard model (SM) of particle physics and cosmology (which is insufficient to explain the matter-anti-matter asymmetry in the universe) imply parity symmetry violation at the electroweak energy scale being possibly manifested through helical (chiral) turbulent motions and/or magnetic fields [18,19]. As expected, such parity-violating turbulent sources will produce circularly polarized GWs [20][21][22][23][24][25][26] analogously to the GWs produced via Chern-Simons coupling [27], and, if detected, the GW polarization can be a clean measure of the deviations from the SM and will provide us with a better understanding of the nature of parity symmetry and its violation.
The detection of circular polarization of the stochastic GW background is a challenging task [28,29], and the planar interferometers cannot measure net polarization in the case of isotropic backgrounds [30,31]. However, the dipolar anisotropy induced by our proper motion with respect to the cosmic reference frame makes it possible to measure the net circular polarization of the stochastic GW background [32,33], and recently it has been shown that the net polarization of GWs could be detected with a signal-to-noise ratio of order one by LISA if the strength of the signal achieves h 2 0 Ω GW ∼ 10 −11 (with Ω GW the fraction between the GW energy density E GW and the critical density today E cr = 3H 2 0 /(8πG), with H 0 = 100h 0 kms −1 Mpc −1 the Hubble parameter today and G is the gravitational constant) [34]. These findings make it extremely important to properly compute all characteristics (such as the amplitude, the spectral shape, and the polarization degree) of the GW signal from primordial helical (chiral) sources.
When computing the GW signal from early-universe turbulent sources, previous studies (with the exception of Ref. [35]) assumed stationary hydrodynamic turbulence and a turbulence duration set by a fraction of the Hubble time at the moment of generation (H −1 ⋆ ), making it possible to use the simplified GW equation with a discarded term that describes the expansion of the universe (∼ H = a −1 da/dt ph , where t ph denotes the physical time, and a is the scale factor). Moreover, the magnetic field and primordial plasma coupling, and correspondingly, the turbulence decay have been neglected, making it impossible to study the GW source dynamics and the temporal dependence of the GW amplitude and spectral characteristics beyond the dilution due to the universe expansion, i.e., E GW ∝ a −4 .
Recently we have performed for the first time direct numerical simulations of magnetohydrodynamic (MHD) turbulence in the early universe accounting for the expansion of the universe using the Pencil Code [36], and numerically computed the resulting stochastic GW background and relic magnetic fields [35]. To scale out expansion effects, we use appropriately scaled comoving variables and conformal time. The full set of MHD equations is then similar to the usual MHD equations [4], and the GW equation is written for the scaled strains and the comoving total (magnetic and kinetic) tracelesstransverse stress tensor T TT ij . These simulations allowed us to conclude that the proper inclusion of coupling effects results in the extension of the GW spectrum at lower frequencies due to the power transfer at large scales, until the causal horizon determined by the comoving Hubble frequency, f H = 1.65 · 10 −5 Hz (g ⋆ /100) 1/6 (T ⋆ /100GeV) with g ⋆ and T ⋆ being the relativistic degrees of freedom and the temperature at the moment of the GW source activation. In terms of the GW energy density per logarithmic frequency interval as a fraction of critical density today, h 2 0 Ω GW (f ), the spectrum in the frequency range and N determines the number of turbulent eddies per linear Hubble length scale) 3 is ∝ f as opposed to f 3 for the causal low-frequency (f < f H ) tail obtained 2 The factor "2" is due to the quadratic nature of the turbulent source. 3 The parameter N is determined by the physical stirring wave number of the source, k phys 0 , via H⋆/N = 2π/k phys 0 , and to the comoving peak wave number through k 0 = k phys 0 (a⋆/a 0 ). Here a⋆ and a 0 are scale factors corresponding to the moment of generation and today, and a⋆/a 0 = 8.0 × 10 −16 (100/g⋆) 1/3 (100GeV/T⋆).
analytically [11], while in terms of the comoving dimensionless strain amplitude h c , conventionally written as , we observe the scaling h c ∝ f −1/2 in the frequencies range f H < f < f S . At this point we distinguish three parts of Ω GW (f ): the low-frequency region below causal horizon f < f H , the intermediate region f H < f < f S , and the high-frequency region f > f S . However, due to computational limitations, we are unable to reproduce the entire spectrum in our numerical simulations. At high frequencies, f > f S , numerical simulations agree well with the analytical estimate; see Ref. [35] for more details: in particular, for Kolmogorov turbulence with a slope −5/3, the high frequency tail h c (f ) (Ω GW (f )) scales [25] showed that the high-frequency scaling of the GW spectrum depends strongly on the assumptions about the turbulence and the modeling of the time-decorrelation function. Since the spectrum of the stochastic GW background determines the detectability of GWs in a given experiment, realistic turbulence simulations are essential for establishing the sensitivity of upcoming GW experiments to early-universe physics [37].
Interestingly, the spectral form of the GW spectrum at low frequencies is independent of the initial conditions and the turbulence model while the amplitude of the signal strongly depends on the model chosen. The universal form of the GW spectrum does not allow us to discriminate between helical and nonhelical sources and thus limits our ability to determine the parity violation in the early universe. This leads us to the present study with its main focus on GW circular polarization estimates and the question whether the detection of polarization can help in the identification of distinct properties of the source.
As it was shown in Ref. [35], the GW spectrum becomes stationary shortly after the driving of the source ends (i.e., when the free decay stage of the source starts), while the energy density of the source is still present. To demonstrate this, we drive magnetic fields with an electromotive force consisting of plane waves that are deltacorrelated in time. In Fig. 1, we show the temporal evolution of the source and the growth of the GW energy density for the driven (1 < t < 1.1) and decaying stages (t > 1.1), where the driving decreases linearly for a time τ = 0.1-2, although τ > 0.5 may be unrealistic. During the statistically stationary stage, the GW energy density growth rate is proportional to the duration of turbulence, as was estimated through analytical modeling of Ref. [10]. In reality, the driving stage is short compared to the Hubble time-scale, and consists of the few largest eddy turnover times.
Moreover, the GW generation can be split into three cases (see Fig. 2 for kinetically and magnetically driven turbulence): non-helical, partially helical (the fractional helicity is described by the ratio σ < 1 of actual to max- imally allowed helicity), and fully helical. Surprisingly, the kinetically driven turbulence is more efficient in producing GW energy; see Figs. 2(c) and (f). However in this case the presence of kinetic helicity does not affect the source amplitude -contrary to the magnetically driven case where the amplitude of the source increases substantially with increasing fractional helicity; see Figs. 2(a) and (d). We also present the kinetic and magnetic energy dissipation rates, ǫ K and ǫ M , respectively. The dissipation rates remain almost unchanged during the stationary stage as we can expect. In addition, we see that they are almost unaffected by the presence of helicity. One may have expected a correlation between ǫ K (or ǫ M ) and Ω GW (f ), but in the magnetically driven case, larger values of σ produce even slightly less dissipation at early times. Nevertheless, Ω GW (f ) clearly increases with σ.
To estimate the polarization degree, we follow the standard procedure described in Refs. [37,38]. Using the usual circular polarization basis tensors, we decompose the Fourier transform of the GW strains into two states -right-(h + ) and left-handed (h − ) circularly polarized GWs, and the polarization degree is defined P characterize the GW amplitude and polarization, respectively.
The GW polarization degree depends, as expected, on the fractional helicity of the source. On the other hand, the evolution of helical sources is determined by the initial fractional helicity: the coupling between the partially helical magnetic field and the plasma motions leads to a reconfiguration of the magnetic field at large scales through free decay, resulting in the growth of the fractional helicity due to the increase of the correlation length and the corresponding decrease of the magnetic energy until the fully helical stage is developed and inverse cascading starts [39]. However, for weakly helical sources, a substantial time-period is needed to reach a fully helical configuration. On the other hand, our simulations show that the dominant contribution to the GW signal occurs shortly after the source has reached its maximum: the subsequent decay of the magnetic field causes a decline of the turbulent driving of GWs. This decline is further enhanced by the expansion of the universe, although this effect is small if the decay time of the turbulence is short compared with the Hubble time; see Fig. 1. We see that, even with a substantially extended decay phase of the turbulence of twice the Hubble time, the final GW production is enhanced by only a factor of about four. Thus, the GW polarization will retain information about the initial fractional helicity of the source.
As we have highlighted above, previous works [20,22] considered stationary turbulence with two different models for helical turbulence realization: (i) Kolmogorovlike helical turbulence with two different spectral slopes, −5/3 and −11/3 for the spectral energy and helicity densities, E M (k) and H M (k), respectively, and (ii) helicitytransferring turbulence (if helicity transfer and smallscale helicity dissipation dominates [40]) with the spectral indices −7/3 and −10/3 for E M (k) and H M (k) [41] 4 . The former case seems most suitable for describing the usual nonhelical turbulence experiencing forward cascading [42]. The difference in using these spectral shapes is determined by the effect of helicity on the energy dissipation length. For highly helical turbulence, the helicity dissipation length is larger, so the two region description [20] might be justified using Kolmogorov-like turbulence 4 Note that Ref. [20] used the different convention for the spectral indices referring to the power spectra for the symmetric (P S (k) ∝ E M (k)/k 2 ) and the antisymmetric parts at large wavenumbers and approximating the low wave number tail with helicity transfer turbulence [43]. Following this description, Ref. [21] and later also Ref. [26] discussed two stages -first the fractionally helical one and later the fully helical one with inverse transfer to describe more precisely the GW generation by helical MHD turbulence. However, these estimates suffer due to (i) the assumption of stationary turbulence; (ii) neglecting the decay and correspondingly temporal dynamical effects (the GW spectrum becomes stationary shortly after source activation). We show the polarization degree spectra in Fig. 3 for continuous pumping of kinetic or magnetic energy and helicity at intermediate scales. We see a substantial difference between the kinetic and magnetic initial sources at the low frequency tail due to the inverse cascade for the magnetic sources that is absent in the kinetically driven case: more precisely the energy density spectrum is unchanged after the decay stage starts, while the transfer of magnetic helicity to the large scales results in the increasing of the polarization degree. These results confirm that the polarization degree is scale dependent: ∝ k −0.5 at large wavenumbers, which is shallower than the k −1 expected for Kolmogorov-like helical turbulence with different spectral indices for the magnetic spectral energy density (n S = −5/3) and the spectral helicity density (n H = −11/3). In our simulations, the actual indices are a bit smaller [38], which also explains the shallower slope in the polarization degree. The departure from the theoretical predictions is due to the assumption of a scaleindependent time-decorrelation function for magnetic energy and helicity densities. Interestingly, the spectral shape of the polarization degree is independent of the actual indices for energy and helicity, but depends on the difference between them [20]. Obviously, in a realistic case, the proper consideration of time-decorrelation and its dependence on wave numbers is required. In fact, even in the simplified description, different forms of the time-decorrelation function for different models of turbulence (including both compressible and incompressible cases) and its scale-dependence leads to different scaling of the GW spectrum at high frequencies [25].
In this Letter we present the first numerical simulations of the circular polarization degree of GWs generated through parity violating turbulent sources in the early universe. We present our results for GWs generated at electroweak energy scales, but the formalism is not limited to the specific moment of GW generation, and can be adjusted to primordial turbulence sources at any time after inflation and before recombination epochs. We have confirmed that the GW signal reaches its maximal strength faster then the turbulence decays. We have also shown that the slope of the low frequency tail of the GW spectrum is independent of the nature of the turbulent source (i.e., the nature of initial driving, the presence of helicity, etc). This restricts the discrimina-tion between helical and non-helical sources as well between kinetic and magnetic drivings if the signal will be detected. On the other hand, the polarization spectrum not only retains information about the initial characteristic frequency (as well as the GW signal does) and the strength of parity violation of the source, but also manifests the dependence on the driving mechanisms. We have shown that the previously used assumption of stationary turbulence does not predict the GW polarization spectrum for realistic turbulence (the scaling of the spectrum at low and high frequencies, the peak position, etc), and thus might result in inadequate estimates of the detection prospects. Fortunately, numerical simulations have now become an affordable tool to address these and other questions of relic GW generation and give a more complete picture for the detection prospects by LISA (as for electroweak phase transitions) and by PTAs (as for QCD phase transitions), and/or any future planned missions, including atomic interferometry [44].
Data availability-The source code used for the simulations of this study, the Pencil Code, is freely available from Ref. [36]. The simulation setups and the corresponding data are freely available from Ref. [45].
We thank Arthur Kosowsky and Andrii Neronov for useful discussions. Support through the Swedish Research Council, grant 2019-04234, and Shota Rustaveli GNSF (grants FR/18-1462 and FR/19-8306) 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.

Supplemental Material 1 Polarization Degree
To compute the polarization degree we follow the standard procedure: choosing the coordinate system so that unit vectorê 3 points in the GW propagation direction (k = k/k with k = |k| is the unit vector) and using the usual circular polarization basis tensors e ± ij = −(e 1 ± ie 2 ) i × (e 1 ± ie 2 ) j / √ 2, we decompose the Fourier transform of the GW strains h ij (k) = d 3 x e ik·x h ij (x) into two states -right-(h + ) and left-handed (h − ) circularly polarized GWs The GW circular polarization degree is given by [20] 5 where H(k) and H(k) characterize the GW amplitude and polarization (chirality) defined through the Gaussiandistributed GWs wave number-space two-point function (for simplicity of notations we omit the time-dependence): , with 4M ijlm (k) ≡ P il P jm + P im P jl − P ij P lm , and 8A ijlm (k) ≡k q (P jm ǫ ilq + P il ǫ jmq + P im ǫ jlq + P jl ǫ imq ) are tensors, P ij (k) ≡ δ ij −k ikj is the projector operator, δ ij and ǫ ijl are the Kronecker delta and the fully antisymmetric tensor, respectively.

Energy and helicity spectra
In Fig. 4 we show the numerator and denominator of the degree of polarization, H(k) and H(k), respectively. The underlying kinetic and magnetic energy and helicity spectra are shown in Fig. 5.  (a,b) kinetically and (c,d) magnetically forced cases with σ = 0 (black) 0.1 (blue), 0.3 (green), 0.5 (orange), and 1 (red). 5 As an alternative we can use the decomposition using the linear polarization basis tensors e + ij (k) = e 1 i e 1 j − e 2 i e 2 j and e × ij (k) = e 1 i e 2 j + e 2 i e 1 j as h ij (k) = h + (k)e + ij (k) + h × (k)e × ij (k) where h + and h × are gauge independent components corresponding to two polarization modes. In this case the polarization degree will be equal to P  Figure 4 shows that the GW polarization is larger at larger length scales. To see whether this could be related to inverse cascading in the magnetic case, we now show in Fig. 5 the corresponding energy and helicity spectra. They show clear inverse cascading of the magnetic energy and helicity spectra in the magnetically driven case. Inverse cascading is clearly absent in the kinetic energy and helicity spectra in the kinetically driven case.