# Dynamical stability and superconductivity of (Ce,La)H9 under high pressure

Prutthipong Tsuppayakorn-aek<EMAIL_ADDRESS>Extreme Condition
Physics Research Laboratory and Center of Excellence in Physics of Energy
Materials (CE:PEM), Department of Physics, Faculty of Science, Chulalongkorn
University, Bangkok 10330, Thailand  Thailand Center of Excellence in Physics
, Ministry of Higher Education, Science, Research and Innovation, 328 Si
Ayutthaya Road, Bangkok 10400, Thailand  Udomsilp Pinsook  Department of
Physics, Faculty of Science, Chulalongkorn University, 10330 Bangkok, Thailand
Thiti Bovornratanaraks  Extreme Condition Physics Research Laboratory and
Center of Excellence in Physics of Energy Materials (CE:PEM), Department of
Physics, Faculty of Science, Chulalongkorn University, Bangkok 10330, Thailand
Thailand Center of Excellence in Physics , Ministry of Higher Education,
Science, Research and Innovation, 328 Si Ayutthaya Road, Bangkok 10400,
Thailand

###### Abstract

Recent experiments have shown that CeH9 and (Ce,La)H9 can be synthesized under
high pressure between 90-170GPa and become a superconductor with a high value
of superconducting critical temperature ($T_{c}$) between 100-200K. In this
work, we performed a theoretical study of a (Ce,La)H9 compound where the Ce:La
ratio is equal to 1. We used the density functional theory and the ab initio
molecular dynamics (AIMD) method. From the phonon dispersion, there exist some
unstable modes around the K-point phonons. Then, we performed AIMD simulation
at around 203K and found that the compound becomes stable. The superconducting
spectral function can be calculated. We found that $\lambda$ is as high as
$3.0$ at 200GPa. By using Allen-Dynes-modified McMillan equation in the strong
coupling regime, we found that $T_{c}=87K$ at 200GPa.

##  I Introduction

In 1968, Ashcroft suggested that metallic hydrogen can become a conventional
superconductor with a high value of critical temperature [1]. He also pointed
out that metal hydrides could become a superconductor [2,3]. Since then, the
search for metallic hydrogen and hydride superconductors have been extensively
explored by both theory and experiments [4-20]. Metal superhydrides have
recently attracted much more of the scientific attention because they can be
synthesized and characterized under high pressure. Furthermore, some of them
have been found to be a superconductor with a close-to-room-temperature
superconducting critical temperature ($T_{c}$). For instance, the Tc of LaH10
has been measured to be 250K under high pressure [21].

In 2019, cerium superhydride CeH9 has been synthesized under high pressure
between 80 and 100GPa by using the laser-heated diamond anvil cell combined
with synchrotron X-ray diffraction [14]. It was found that CeH9 is a hexagonal
structure with a space group of P63/mmc. Its $T_{c}$ has been theoretically
estimated to be 105 - 117K [14]. Recently, the (La,Ce)Hx compounds have been
successfully synthesized at pressure lower than 130GPa by Chen, et al. [22].
They discovered that the compound resembles (La,Ce)H9 with $T_{c}$ equal to
176K at 100GPa.

At the current stage, the trend has moved from searching for the
superconductivity of a binary compound to searching for the superconductivity
of a ternary hydride compound. This is because the mixing of different atomic
species in the making of ternary compounds offers a slight and perhaps
controllable modification to the electronic and dynamical properties of the
hydride compounds, which sometimes tends to enhance $T_{c}$ [12, 17, 20,
23-26]. For example, Mg can be used as a substitution for Ca in the CaH6
compound [17]. It is found that Mg0.5Ca0.5H6 compound is stable at 200GPa and
become a superconductor with $T_{c}\approx$ 288K.

In this work, we performed the simulation of the Ce/La substituted of
(Ce,La)H9 at 200GPa, as both CeH9 and LaH9 share a common hexagonal structure
with a space group of P63/mmc [14, 27, 28]. The hexagonal structure is also
suggested by the experiment for the La-Ce-H compounds [22]. We calculated the
electron-phonon interaction and evaluated the Eliasberg spectral function by
the density functional theory (DFT). The $T_{c}$ can then be approximated. We
also preformed the ab initio molecular dynamics (AIMD) to examine the
dynamical stability of (Ce,La)H9.

##  II Computational details

In this work, we performed the DFT calculations with the generalized gradient
approximation of the Perdew–Burke–Ernzerhof (GGA-PBE) functional [29] for the
exchange-correlation functional, and the projector augmented wave (PAW) method
[30], as implemented in the quantum espresso (QE) [31,32]. For the calculation
of the electron-phonon interaction, a plane-wave energy cutoff of 60Ry was
used. The electron-phonon coupling (EPC) matrix elements were computed in the
first Brillouin-zone (BZ) on 4$\times$4$\times$2 q-meshes using individual EPC
matrices obtained with 24$\times$24$\times$16 k–point meshes. The Allen-Dynes-
modified McMillan equation [33] was used for the estimating Tc.

For the electronic structure, we used a plane-wave basis set with the cutoff
energy of 600eV and the first BZ with 12$\times$12$\times$8 k-point meshes.
The pseudocore radii of Ce, La, and H are 2.57Bohr, 2.80Bohr, and 1.1Bohr,
which are small enough to ensure that no core overlapping will occur at
200GPa. For ab initio molecular dynamics (AIMD), the simulation was performed
by using the $NPT$ ensemble with 160 atoms per supercell. The AIMD simulation
was carried out at 203K and 200GPa. Both the electronic structure and AIMD
were calculated by using the Vienna ab initio simulation package (VASP) [34].

##  III Results and discussion

As both CeH9 and LaH9 share a common hexagonal structure with a space group of
P63/mmc [14, 27, 28], as shown in Fig. 1(a). The unit cell of the P63/mmc
structure contains two metal sites. We modelled the (Ce,La)H9 compound by
inserting an La atom into one of the metal site, and a Ce atom into the
others. The Ce:La ratio is equal to 1. The symmetry is reduced from the space
group of P63/mmc to P${\bar{6}}$m2, as shown in Fig. 1(b).

Figure 1: (a) The hexagonal structure of CeH9, where the yellow spheres
represent the Ce atoms and the pink spheres represent the H atoms. (b) The
hexagonal structure of (Ce,La)H9 , where the yellow sphere represents the Ce
atom, the green sphere represents the La atom, and the pink spheres represent
the H atoms.

The electronic structures are shown in Fig. 2(a). Under high pressure, the
(Ce,La)H9 compound becomes metallic. The band dispersions exhibit the weaving
of the up-running bands from below EF and down-running bands above EF, similar
to those of CeH9 [14], CeH10 [35], and LaH10 [36]. The corresponding projected
density of states (PDOS) of (Ce,La)H9 is also shown in Fig. 2(a). There are
several places near the Fermi level where the electronic structures exhibit
flat dispersions, such as around the H-point and M-point. These flat
dispersions lead to the so-called van Hove singularity (vHs). The vHs is
marked by an arrow in Fig. 2(a). The Fermi surface topology (FST) is shown in
Fig. 2(b). The FST exhibits a number of parallel surfaces, which should lead
to Fermi surface nesting and enhance superconductivity.

Figure 2: Electronic structures of (Ce,La)H9 at 200 GPa. (a) The band
dispersions (left) and the projected density of states (PDOS). (b) The Fermi
surface topology (FST)

As Tc can be related to the number of electronic states at the Fermi surface,
it is useful to consider the local density of states (LDOS) [37,38], which is
defined as

$N(E,\mathbf{r})=\sum_{n}\int\frac{d^{3}k}{(2\pi)^{3}}\delta(E-\epsilon_{nk})|\psi_{nk}(\mathbf{r})|^{2},$
(1)

where $\psi$($\mathbf{k}$)($\mathbf{r}$) and $\epsilon$nk are the Kohn-Sham
eigenfunctions and eigenvalues of the system. The LDOS allows us to evaluate
the number of the electronic states at the Fermi level, which can be related
to electron pairings in the superconducting state. The LDOS indicates that the
Fermi surface of the (Ce,La)H9 compound contains 0.91 states/eV.

Next, we calculated the isotropic Eliashberg spectral function, which is
written as

$\alpha^{2}F(\omega)=\frac{1}{N_{F}}\sum_{\nu,k,k^{\prime}}|g^{\nu}_{kk^{\prime}}|^{2}\delta(\epsilon_{k})\delta(\epsilon_{k}^{\prime})\delta(\omega-\omega_{q\nu}),$
(2)

where NF is the density of states at the Fermi level, k and k′ are the
electronic states, and q are the phonon modes. Nk (Nq) is the total number of
k (q) points in the simulation cell, $\epsilon$k is the energy eigenvalue of
the Kohn-Sham state with respect to the Fermi level, and
g${}^{\nu}_{kk^{\prime}}$ is the electron-phonon matrix elements for the
scattering between the electronic states k and k′ via a phonon with wave
vector $q=k^{\prime}-k$, and $\omega$qν are the frequencies of phonon modes q
and phonon band $\nu$. The calculated Eliashberg spectral function
$\alpha$2$F(\omega)$ is shown in the right panel of Fig. 3.

The corresponding phonon dispersions are also shown in the left panel of Fig.
3. It is obvious that the acoustic modes are well-seperated from the optical
modes. The optical modes are in the frequency range between 360cm-1 and
1814cm-1. The acoustic modes are in the frequency range between 0-281cm-1, and
contain several imaginary frequencies around the K-point phonon modes, as
shown in the left panel of Fig. 3. These modes indicate that the (Ce,La)H9
compound with the hexagonal structure would indeed be an unstable structure.
We will reexamine this point in more detail later.

Figure 3: (Left) Phonon dispersion curves of (Ce,La)H9 at 200GPa. (Right) The
Eliashberg spectral function $\alpha$2$F(\omega)$ and the corresponding
accumulated $\lambda(\omega)$ at 200GPa.

From $\alpha$2$F(\omega)$, we can evaluate the average coupling strength
$\lambda$ as a function of $\omega$ by

$\lambda(\omega)=2\int_{0}^{\omega}d\omega^{\prime}\frac{\alpha^{2}F(\omega^{\prime})}{\omega^{\prime}},$
(3)

This is the so-called accumulated $\lambda(\omega)$. The result of
$\lambda(\omega)$ is also shown in the right panel of Fig. 3. The total
$\lambda$ can be obtained by integrating over all available frequencies. We
found that for (Ce,La)H9, $\lambda$ is as high as 3.0. Thus, the electron-
phonon interaction in (Ce,La)H9 enters the strong coupling regime. As a
consequence, we chose the Allen-Dynes-modified McMillan equation with strong
coupling correction for describing the superconducting critical temperature,
as follows;

$T_{c}=f_{1}f_{2}\frac{\omega_{log}}{1.2}\exp\Big{[}-\frac{1.04(1+\lambda)}{\lambda-\mu^{*}(1+0.62\lambda)}\Big{]},$
(4)

We chose $\mu^{*}=0.1$. We found that f1 and f2 are 1.22 and 1.20,
respectively. The $\omega_{log}$ is 327K. We found that Tc is 87K. Despite of
the large value of $\lambda$ and NF, the $\omega_{log}$ is very low. This is
because the bandwidth of $\alpha$2$F(\omega)$ is very large [15]. As a result,
Tc becomes not as high as we might expect.

At this stage, we discuss the structure stability. From the phonon dispersion
in the left panel of Fig. 3, there exist some unstable phonon modes at 0K. In
order to take the effects of temperature into the account, we performed the
AIMD simulation, with the NPT ensemble. The simulation cell was construct from
a 2$\times$2$\times$2 hexagonal supercell with 160 atoms. We chose temperature
of 203K and pressure of 200GPa. The energy as a function of time is shown in
Fig. 4 (a). The time step is 1 fs. We performed the AIMD integration to 1,000
time step. After the simulation is in equilibrium, we took several snapshots,
as shown in Fig. 4 (b)-(e). The hexagonal structure remains stable during the
simulation time. The AIMD simulation suggests that the effects of temperature
would help stabilize the hexagonal structure of the (Ce,La)H9 compound.

Figure 4: (a) Energy as a function of time of (Ce,La)H9 from AIMD simulations
(b-d) the snapshots of the AIMD simulations at different time steps.

##  IV Conclusion

We studied the dynamical stability and superconductivity of the (Ce,La)H9
compound. The symmetry is reduced from the space group of P63/mmc, as of CeH9
and LaH9, to $P{\bar{6}}m2$. The phonon dispersions showed some unstable
phonon modes around the K-point phonons. However, the Elaisberg spectral
function can still be evaluated. We found that $\lambda$ is as high as 3.0 and
$N(E_{F})$ is as large as 0.91 states/eV. Despite of high values of these
favorable parameters, the $\omega_{log}$ appears to be very low. This results
in a low $T_{c}$ = 87K.

##  V ACKNOWLEDGMENTS

We acknowledge the Computational Materials Physics project SLRI, Thailand for
providing support for computational resources. We thank the Chalawan HPC
cluster, which is operated and maintained by the National Astronomical
Research Institute of Thailand (NARIT) under the Ministry of Science and
Technology of the Royal Thai government. This research project was supported
by the Second Century Fund (C2F), Chulalongkorn University. This project is
funded by National Research Council of Thailand (NRCT): (NRCT5-RSA63001-04).
This research is partially funded by Chulalongkorn University; Grant for
Research.

##  VI References

[1] N. W. Ashcroft, Metallic hydrogen: A high-temperature superconductor?,
Phys. Rev. Lett. 21, 1748 (1968).

[2] N. W. Ashcroft, Bridgman’s high-pressure atomic destructibility and its
growing legacy of ordered states, J. Condens. Matter Phys. 16, S945 (2004).

[3] N. W. Ashcroft, Hydrogen dominant metallic alloys: High temperature
superconductors?, Phys. Rev. Lett. 92, 187002 (2004).

[4] E. Zurek, R. Hoffmann, N. W. Ashcroft, A. R. Oganov, and A. O. Lyakhov, A
little bit of lithium does a lot for hydrogen, Proc. Natl. Acad. Sci. U.S.A.
106, 17640 (2009).

[5] T. A. Strobel, A. F. Goncharov, C. T. Seagle, Z. Liu, M. Somayazulu, V. V.
Struzhkin, and R. J. Hemley, High pressure study of silane to 150 GPa, Phys.
Rev. B 83, 144102 (2011).

[6] A. Shamp and E. Zurek. Superconducting high pressure phases composed of
hydrogen and iodine, J. Phys. Chem, 6, 4067 (2015).

[7] I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y. Li, H.
Liu, Y. Zhang, Ya. Ma, and F. Mauri, High pressure hydrogen sulfide from first
principles: A strongly anharmonic phonon-mediated superconductor, Phys. Rev.
Lett. 114, 157004 (2015).

[8] A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I.
Shylin, Conventional superconductivity at 203 kelvin at high pressures in the
sulfur hydride system, Nature, 525, 73 (2015).

[9] C. M. Pépin, G. Geneste, A. Dewaele, M. Mezouar, and P. Loubeyre.
Synthesis of FeH5: A layered structure with atomic hydrogen slabs, Science,
357, 382 (2017).

[10] B. Liu, W. Cui, J. Shi, L. Zhu, J. Chen, S. Lin, R. Su, J. Ma, Ka. Yang,
M. Xu, J. Hao, A. P. Durajski, J. Qi, Y. Li, and Y. Li, Effect of covalent
bonding on the superconducting critical temperature of the H-S-Se system,
Phys. Rev. B 98, 174101 (2018).

[11] C. Heil, S. Cataldo, G. B. Bachelet, and L. Boeri. Superconductivity in
sodalite-like ytttrium hydride clathrates, Phys. Rev. B 99, 220502 (2019).

[12] Y. Sun, J. Lv, Y. Xie, H. Liu, and Y. Ma, Route to a superconducting
phase above room temperature in electron-doped hydride compounds under high
pressure, Phys. Rev. Lett. 123, 097001 (2019).

[13] X. Liang, A. Bergara, L. Wang, B. Wen, Z. Zhao, X. F. Zhou, J. He, G.
Gao, and Y. Tian. Potential high-Tc superconductivity in CaYH12 under
pressure, Phys. Rev. B 99, 100505 (2019).

[14] N. P. Salke, M. M. D. Esfahani, Y. Zhang, I. A. Kruglov, J. Zhou, Y.
Wang, E. Greenberg, V. B. Prakapenka, J. Liu, A. R. Oganov, and J. F. Lin,
Synthesis of clathrate cerium superhydride CeH9 at 80-100 GPa with atomic
hydrogen sublattice, Nat. Commun, 10, 4453 (2019).

[15] U. Pinsook, In search for near-room-temperature superconducting critical
temperature of metal superhydrides under high pressure: A review. J. Met.
Mater. Miner, 30(2) (2020).

[16] D. V. Semenok, A. G. Kvashnin, A. G. Ivanova, V. Svitlyk, V. Yu.
Fominski, A. V. Sadakov, O. A. Sobolevskiy, V. M. Pudalov, I. A. Troyan, and
A. R. Oganov. Superconductivity at 161 K in thorium hydride ThH10: Synthesis
and properties, Mater. Today, 33, 36 (2020).

[17] W. Sukmas, P. Tsuppayakorn-aek, U. Pinsook, and T. Bovornratanaraks,
Near-room-temperature superconductivity of Mg/Ca substituted metal hexahydride
under pressure. J. Alloys Compd, 849, 156434 (2020).

[18] P. Tsuppayakorn-aek, U. Pinsook, W. Luo, R. Ahuja, and T.
Bovornratanaraks, Superconductivity of superhydride CeH10 under high
pressure.Mater. Res. Express, 7, 086001 (2020).

[19] H. Xie, Y. Yao, X. Feng, D. Duan, H. Song, Z. Zhang, S. Jiang, S. A. T.
Redfern, V. Z. Kresin, C. J. Pickard, and T. Cui, Hydrogen pentagraphenelike
structure stabilized by hafnium: A high-temperature conventional
superconductor, Phys. Rev. Lett. 125, 217001 (2020).

[20] P. Tsuppayakorn-aek, P. Phansuke, P. Kaewtubtim, R. Ahuja, and T.
Bovornratanaraks, Enthalpy stabilization of superconductivity in an alloying
S-P-H system: First-principles cluster expansion study under high pressure.
Comput. Mater. Sci, 190, 110282 (2021).

[21] A. P. Drozdov, P. P. Kong1, V. S. Minkov, S. P. Besedin, M. A.
Kuzovnikov, S. Mozaffari, L. Balicas, F. F. Balakirev, D. E. Graf, V. B.
Prakapenka, E. Greenberg, D. A. Knyazev, M. Tkacz, and M. I. eremets,
Superconductivity at 250 K in lanthanum hydride under high pressures. Nature,
569, 528 (2019).

[22] W. Chen, X. Huang, D. V. Semenok, S Chen, K. Zhang, A. R. Oganov, and T.
Cui, Enhancement of superconducting critical temperature realized in La-Ce-H
system at moderate pressures, arXiv:2203.14353 (2022).

[23] F. Fan, D.A. Papaconstantopoulos, M.J. Mehl, and B.M. Klein, High-
temperature superconductivity at high pressures for H3SixP1-x, H3PxS1-x, and
H3ClxS1-x, J. Phys. Chem. Solids, 99, 105 (2016).

[24] Y. Ge, F. Zhang, and Y. Yao, First-principles demonstration of
superconductivity at 280 K in hydrogen sulfide with low phosphorus
substitution, Phys. Rev. B 93, 224513 (2016).

[25] A. Nakanishi, T. Ishikawa, and K. Shimizu, First-principles study on
superconductivity of P- and Cl-doped H3S, J. Phys. Soc. Jpn. 87, 124711
(2018).

[26] A. P. Durajski and R. Szczęśniak. Gradual reduction of the
superconducting transition temperature of H3S by partial replacing sulfur with
phosphorus, Physica C Supercond, 554, 38 (2018).

[27] I. A. Kruglov, D. V. Semenok, H. Song, R. Szczęśniak, I. A. Wrona, R.
Akashi, M. M. D. Esfahani, D. Duan, T. Cui, A. G. Kvashnin, and A. R. Oganov,
Superconductivity of LaH10 and LaH16 polyhydrides, Phys. Rev. B 101, 024508
(2020).

[28] A. M. Shipley, M. J. Hutcheon, M. S. Johnson, R. J. Needs, and C. J.
Pickard, Stability and superconductivity of lanthanum and yttrium
decahydrides, Phys. Rev. B, 101, 224511 (2020).

[29] J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient
approximation made simple. Phys. Rev. Lett. 77, 3865 (1996).

[30] P. E. Blöchl. Projector augmented-wave method. Phys. Rev. B 50, 17953
(1994).

[31] S. Baroni, S. D. Gironcoli, A. D. Corso, and P. Giannozzi, Phonons and
related crystal properties from density-functional perturbation theory. Rev.
Mod. Phys, 73, 515 (2001).

[32] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D.
Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. D.
Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A.
Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S.
Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero,
A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Quantum
espresso: a modular and open-source software project for quantum simulations
of materials. J. Condens. Matter Phys, 21, 395502, (2009).

[33] P. B. Allen and R. C. Dynes, Transition temperature of strong-coupled
superconductors reanalyzed, Phys. Rev. B, 12, 905 (1975).

[34] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio
total-energy calculations using a plane-wavebasis set, Phys. Rev. B 54, 11169
(1996).

[35] L. Liu, C. Wang, S. Yi, K. W. Kim, J. Kim, and Jun-Hyung Cho. Microscopic
mechanism of room-temperature superconductivity in compressed LaH10, Phys.
Rev. B 99, 140501 (2019).

[36] C. Heil, G. B. Bachelet, and L. Boeri, Absence of superconductivity in
iron polyhydrides at high pressures. Phys. Rev. B 97, 214510 (2018).

[37] S. Cataldo, W. Linden, and L. Boeri, Phase diagram and superconductivity
of calcium borohyrides at extreme pressures, Phys. Rev. B 102, 014516 (2020).


# Circular polarization of gravitational waves from early-Universe helical
turbulence

Tina Kahniashvili<EMAIL_ADDRESS>McWilliams Center for Cosmology and
Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA
Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi,
Georgia  Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia
Department of Physics, Laurentian University, Sudbury, ON P3E 2C, Canada  Axel
Brandenburg<EMAIL_ADDRESS>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  Grigol Gogoberidze
<EMAIL_ADDRESS>Faculty of Natural Sciences and Medicine,
Ilia State University, 0194 Tbilisi, Georgia  Sayan Mandal
<EMAIL_ADDRESS>Physics and Astronomy Department, Stony Brook
University, Stony Brook, New York 11794, USA  Faculty of Natural Sciences and
Medicine, Ilia State University, 0194 Tbilisi, Georgia  Alberto Roper Pol
<EMAIL_ADDRESS>Université de Paris, CNRS, Astrophysique et Cosmologie,
Paris, F-75013, France  Faculty of Natural Sciences and Medicine, Ilia State
University, 0194 Tbilisi, Georgia

###### Abstract

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 turbulent geneses such as
magnetically or kinetically driven cases. Under realistic physical conditions
in the early Universe we compute numerically the wave number-dependent
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 of nearly 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 activation of the source. 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:

98.70.Vc, 98.80.-k

††preprint: NORDITA-2020-102

##  I Introduction

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 Witten:1984rs
; Hogan:1986qda ; Kamionkowski:1993fg , or from the primordial magnetic fields
that are coupled to the cosmological plasma Brandenburg:1996fc ;
Christensson:2000sp ; Kahniashvili:2010gp ; Brandenburg:2017rnt . 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 the turbulence is up to
1–10% of the total thermal energy at the moment of generation) Kosowsky:2001xp
; Caprini:2006jb ; Gogoberidze:2007an ; Kahniashvili:2008pe ; Caprini:2009yp .
Since GWs propagate almost freely from the moment of generation until
today111We discard the GW damping due to neutrino free streaming Durrer:1997ta
; Weinberg:2003ur  or from anisotropic stresses Deryagin:1986qq ., 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. Caprini:2019egz  and
references therein. Moreover, several theoretical extensions of the standard
model (SM) of particle physics and cosmology (which is insufficient to explain
the matter-antimatter asymmetry in the Universe) imply parity symmetry222 The
baryon asymmetry of the Universe can be explained through spontaneous lepton
number symmetry breaking at a cosmological phase transition Cohen:1990it ;
Cohen:1990py . violation at the electroweak energy scale being possibly
manifested through helical (chiral) turbulent motions and/or magnetic fields
Long:2013tha ; Dorsch:2016nrg . As expected, such parity-violating turbulent
sources will produce circularly polarized GWs Kahniashvili:2005qi ;
Kahniashvili:2008er ; Kisslinger:2015hua ; Alexander:2018fjp ; Anand:2018mgf ;
Niksa:2018ofa ; Ellis:2020uid  analogously to the GWs produced via Chern-
Simons coupling Alexander:2004us ; Lyth:2005jf . Furthermore, chiral
inflationary GWs might be responsible (through the gravitational anomaly) for
gravitational leptogenesis which, in turn, manifests itself in the neutrino
sector by being successful in either the Dirac or Majorana neutrino mass
scenario Adshead:2017znw . Moreover, the gravitational anomaly effect and,
correspondingly, a successful electroweak baryogenesis scenario appear viable
in the case when the GWs are generated through helical magnetic fields
Abedi:2018top . Additionally, the primordial magnetic fields that are
responsible for the circularly polarized GWs and, correspondingly, for the
baryon asymmetry, can serve as seeds for large-scale magnetic fields in the
Universe; see Ref. Vachaspati:2020blt  for a review and references therein.
Under this consideration the scheme is as follows: helical magnetic fields act
as seeds for large-scale magnetic fields, produce circularly polarized GWs,
and enable successful lepto- and baryogenesis through the gravitational
anomaly. Obviously, 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. In fact, the consideration of
the different forms of “driving” is dictated by different possible sources:
for scenarios in which the magnetic field is the primary (dominant) source of
turbulence (e.g., inflationary or phase transition-generated magnetic fields)
we deal with “magnetically driven” turbulence, while turbulence with primary
fluid motions (e.g., sound waves) is “kinetically driven”, and could then
result in magnetic field generation through dynamo action. This can lead to a
growth of the magnetic energy by orders of magnitude even though the kinetic
energy of the turbulence decays Brandenburg:2017rnt . However, this is not
considered in the present paper.

Each of these sources might be parity symmetry violating ones. As we show
below, the polarization spectra are different for different forms of driving,
and thus the measurement of these spectra can lead to an identification of the
source and the energy-scale at which parity violation happened. Interestingly,
our results may help us understand the origin of the cosmic seed magnetic
field, i.e., it may help discriminate between astrophysical and cosmological
magnetogenesis scenarios. Perhaps most importantly, the investigation of the
GW background in combination with its polarization might convincingly
establish that observations demand a cosmological magnetic field that cannot
be generated by a mechanism operating within the confines of the SM.

The detection of circular polarization of the stochastic GW background is a
challenging task Seto:2006hf ; Smith:2016jqs ; Masui:2017fzw , and the planar
interferometers cannot measure net polarization in the case of isotropic
backgrounds Seto:2007tn ; Seto:2008sr . 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 Seto:2006dz ; Crowder:2012ik , 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_{0}^{2}\Omega_{\rm
GW}\sim 10^{-11}$ (with $\Omega_{\rm GW}$ the fraction between the GW energy
density and the critical density today ${\mathcal{E}}_{\rm
cr}=3H_{0}^{2}/(8\pi G)$, with $H_{0}=100h_{0}\,{\rm km}\,{\rm s}^{-1}{\rm
Mpc}^{-1}$ the Hubble parameter today and $G$ is the gravitational constant)
Domcke:2019zls . 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. Pol:2019yex ) assumed stationary
hydrodynamic turbulence and a turbulence duration set by a fraction of the
Hubble time at the moment of generation ($H_{\star}^{-1}$), making it possible
to use the simplified GW equation with a discarded term that describes the
expansion of the Universe ($\sim H=a^{-1}{da/dt_{\rm ph}}$, where $t_{\rm 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 in those analytic studies, 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
expansion of the Universe. These shortcomings will be avoided by numerically
simulating nonstationary turbulence with all terms included, allowing for a
full coupling between magnetic fields and plasma motions. We apply
hydrodynamic and magnetic forcing terms that have a realistic representation
of the time dependent generation of turbulence.

Figure 1: Evolution of magnetic energy (left) and growth of GW energy density
(right) for simulations where the driving is turned off at $t=1.1$ (black
dotted line), or the strength of the driving is reduced linearly in time over
the duration $\tau=0.2$ (green), $0.5$ (blue), $1$ (red), or $2$ (black). Time
is in units of the Hubble time at the moment of source activation. The
magnetic and GW energy densities are in units of the radiation energy density.

##  II Numerical Modeling

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 PC , and numerically computed
the resulting stochastic GW background and relic magnetic fields Pol:2019yex .
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 (Brandenburg:1996fc, ), and the GW equation is written for the
scaled strains and the comoving total (magnetic and kinetic) traceless-
transverse stress tensor $T_{ij}^{\rm TT}$. 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\cdot 10^{-5}{\rm
Hz}\,(g_{\star}/100)^{1/6}(T_{\star}/100{\rm GeV})$ with $g_{\star}$ and
$T_{\star}$ 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_{0}^{2}\Omega_{\rm GW}(f)$, the spectrum in the frequency range $f_{\rm
H}\leq f\leq f_{\rm S}$ (where the typical frequency of the source is $f_{\rm
S}=2Nf_{\rm H}$333The factor “2” is due to the quadratic nature of the
turbulent source., and $N$ determines the number of turbulent eddies per
linear Hubble length scale)444The parameter $N$ is determined by the physical
stirring wave number of the source, $k_{0}^{\rm phys}$, via
$H_{\star}/N=2\pi/k_{0}^{\rm phys}$, and to the comoving peak wave number
through $k_{0}=k_{0}^{\rm phys}(a_{\star}/a_{0})$. Here $a_{\star}$ and
$a_{0}$ are scale factors corresponding to the moment of generation and today,
and $a_{\star}/a_{0}=8.0\times 10^{-16}(100/g_{\star})^{1/3}(100\,{\rm
GeV}/T_{\star})$. is $\propto f$, as opposed to $f^{3}$ for the causal low-
frequency ($f<f_{H}$) tail obtained analytically Kahniashvili:2008pe , while
in terms of the comoving dimensionless strain amplitude $h_{c}$,
conventionally written as $h_{c}(f)=1.263\cdot
10^{-18}(f/1{\rm\,Hz})^{-1}\left[h_{0}\Omega_{\rm GW}(f)\right]^{1/2}$
Romano:2016dpx , we observe the scaling $h_{c}\propto f^{-1/2}$ in the
frequency range $f_{H}<f<f_{S}$. At this point we distinguish three parts of
$\Omega_{\rm GW}(f)$: the low-frequency region below the 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_{\rm S}$, numerical simulations agree well with the
analytical estimate. In particular, for Kolmogorov turbulence with a slope
$-5/3$, the high frequency tail $h_{c}(f)$ ($\Omega_{\rm GW}(f)$) scales as
$\propto f^{-7/3}$ ($\propto f^{-8/3}$); see Ref. Pol:2019yex  for more
details. Earlier work Niksa:2018ofa  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
Romano:2016dpx .

Figure 2: Evolution of (a) ${\cal E}_{\rm K}$, (b) $\epsilon_{\rm K}$, and (c)
${\cal E}_{\rm GW}$ for kinetically driven cases with $\sigma=0$ (Run K0,
black), $0.5$ (Run K05, blue), and $1$ (Run K1, red), and of (d) ${\cal
E}_{\rm M}$, (e) $\epsilon_{\rm M}$, and (f) ${\cal E}_{\rm GW}$ for
magnetically driven cases with $\sigma=0$ (Run M0, black), $0.3$ (Run M03,
blue), and $1$ (Run M1, red).

##  III Results

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.

Table 1:  Characteristic parameters of the runs of Fig. 1.

Run ${\cal P}$ ${\cal E}_{\rm M}^{\max}$ ${\cal E}_{\rm K}^{\max}$ ${\cal
E}_{\rm GW}^{\rm sat}$ ${h}_{\rm rms}^{\rm sat}$ $\lambda=0.1$ $-13.5$
$0.0367$ $0.0134$ $6.3\times 10^{-8}$ $2.5\times 10^{-6}$ $\lambda=0.2$
$-9.37$ $0.0368$ $0.0135$ $7.6\times 10^{-8}$ $2.8\times 10^{-6}$
$\lambda=0.5$ $-3.32$ $0.0372$ $0.0137$ $1.1\times 10^{-7}$ $3.4\times
10^{-6}$ $\lambda=1$ $-1.43$ $0.0378$ $0.0140$ $1.6\times 10^{-7}$ $4.1\times
10^{-6}$ $\lambda=2$ $-0.63$ $0.0381$ $0.0141$ $2.4\times 10^{-7}$ $5.1\times
10^{-6}$

###  III.1 Decay with Decreased Driving

We first address the temporal evolution of the GW spectrum. As it was shown in
Ref. Pol:2019yex , 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, ${\cal
F}_{i}=(\delta_{ij}+\sigma\epsilon_{ijl}\hat{k}_{l}){\cal F}_{j}^{(0)}$,
consisting of plane waves that are delta-correlated in time. Here,
$-1\leq\sigma\leq 1$ quantifies the fractional helicity, and ${\cal
F}_{j}^{(0)}$ is a nonhelical plane-wave forcing term. Plasma motions are
selfconsistently driven by the Lorentz force. Purely hydrodynamic motions are
driven by a ponderomotive force analogous to ${\cal F}_{i}$.

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 duration $\tau=0.1$–$2$, although
$\tau>0.5$ may be unrealistic. As in Ref. Pol:2019yex , we use $1152^{3}$
meshpoints for the runs in Fig. 1 and put $\sigma=0$; see Table
LABEL:Tsummary_ramp for a summary relevant quantities. 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. Gogoberidze:2007an . In reality, the driving stage
is short compared to the Hubble time-scale, and consists of the few largest
eddy turnover times.

In Table LABEL:Tsummary_ramp, we have quoted the values of ${\cal E}_{\rm
GW}^{\rm sat}$ and ${h}_{\rm rms}^{\rm sat}$ obtained at the end of the
simulation at $t=t_{\rm end}$. We emphasize that ${\cal E}_{\rm GW}$ is the
comoving GW energy density normalized by the critical energy density, which is
the same as the radiation energy density during the simulation, and ${h}_{\rm
rms}$ corresponds to the scaled strain. To compute the relic observable
$h_{0}^{2}{\Omega}_{\rm GW}$ at the present time, we have to multiply ${\cal
E}_{\rm GW}^{\rm sat}$ by a factor $(H_{\ast}/H_{0})^{2}a_{0}^{-4}$; see Refs.
Pol:2019yex ; Pol:2018pao  for details, although there we used the symbol
${\Omega}_{\rm GW}$ also for the latter. For the aforementioned fiducial
parameters of $g_{*}=100$ and $T_{*}=100\,{\rm GeV}$, this factor is
$1.64\times 10^{-5}$. The largest value of ${\cal E}_{\rm GW}^{\rm sat}$
quoted in Table LABEL:Tsummary_ramp is $2.4\times 10^{-7}$ and corresponds
therefore to $h_{0}^{2}{\Omega}_{\rm GW}=4\times 10^{-12}$. Likewise, the
values of ${h}_{\rm rms}^{\rm sat}$ in Table LABEL:Tsummary_ramp have to be
multiplied by $a_{0}^{-1}=8.0\times 10^{-16}$ to obtain the observable $h_{c}$
at the present time. Again, the largest value of ${h}_{\rm rms}^{\rm
sat}=5.1\times 10^{-6}$ corresponds therefore to the observable $h_{c}=4\times
10^{-21}$.

Table 2:  Characteristic parameters of the runs of Fig. 1.

Run ${\cal P}$ ${\cal E}_{\rm M}^{\max}$ ${\cal E}_{\rm K}^{\max}$ ${\cal
E}_{\rm GW}^{\rm sat}$ ${h}_{\rm rms}^{\rm sat}$ K0 $0.01$ $0$ $0.0080$
$4.5\times 10^{-9}$ $7.0\times 10^{-7}$ K01 $0.31$ $0$ $0.0082$ $4.4\times
10^{-9}$ $6.9\times 10^{-7}$ K03 $0.73$ $0$ $0.0085$ $4.8\times 10^{-9}$
$6.4\times 10^{-7}$ K05 $0.88$ $0$ $0.0091$ $5.3\times 10^{-9}$ $5.7\times
10^{-7}$ K1 $0.95$ $0$ $0.0114$ $6.4\times 10^{-9}$ $5.5\times 10^{-7}$ M0
$-0.01$ $0.0059$ $0.0020$ $1.6\times 10^{-9}$ $4.1\times 10^{-7}$ M01 $0.57$
$0.0078$ $0.0021$ $2.0\times 10^{-9}$ $4.6\times 10^{-7}$ M03 $0.86$ $0.0119$
$0.0024$ $3.5\times 10^{-9}$ $6.9\times 10^{-7}$ M05 $0.94$ $0.0148$ $0.0026$
$4.5\times 10^{-9}$ $8.6\times 10^{-7}$ M1 $0.97$ $0.0168$ $0.0025$ $5.3\times
10^{-9}$ $9.0\times 10^{-7}$

###  III.2 Approach to a Stationary State

The GW generation can be split into three cases (see Fig. 2 for kinetically
and magnetically driven turbulence): non-helical, partially helical
($|\sigma|<1$), and fully helical. Here, we used $512^{3}$ meshpoints and
arranged the forcing amplitude such that ${\cal E}_{\rm K}$ and ${\cal E}_{\rm
M}$ are around $10^{-2}$. The viscosity and magnetic diffusivity have equal
values and are chosen such that the dissipative subrange is resolved; see
Table LABEL:Tsummary_pol for a summary and Ref. DATA  for further details.
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, $\epsilon_{\rm
K}$ and $\epsilon_{\rm 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 $\epsilon_{\rm K}$ (or $\epsilon_{\rm M}$) and
${\cal E}_{\rm GW}$, but in the magnetically driven case, larger values of
$\sigma$ produce even slightly less dissipation at early times. Nevertheless,
${\cal E}_{\rm GW}$ clearly increases with $\sigma$.

###  III.3 Polarization Degree

To estimate the polarization degree, we follow the procedure described in Ref.
Pol:2018pao . We use the usual circular polarization basis tensors
$e^{\pm}_{ij}=-({\bf e}_{1}\pm i{\bf e}_{2})_{i}\times({\bf e}_{1}\pm i{\bf
e}_{2})_{j}/\sqrt{2}$, we decompose the Fourier transform of the GW strains
$h_{ij}({\bf k})=\int d^{3}x\,e^{i{\bf k}\cdot{\bf x}}h_{ij}({\bf x})$ into
two states – right- ($h_{+}$) and left-handed ($h_{-}$) circularly polarized
GWs $h_{ij}=h_{+}e^{+}_{ij}+h_{-}e^{-}_{ij}$. The GW circular polarization
degree is given by Kahniashvili:2005qi 555As an alternative we can use the
decomposition using the linear polarization basis tensors
$e^{+}_{ij}({\bf\hat{k}})=e_{i}^{1}e_{j}^{1}-e_{i}^{2}e_{j}^{2}$ and
$e^{\times}_{ij}({\bf\hat{k}})=e_{i}^{1}e_{j}^{2}+e_{i}^{2}e_{j}^{1}$ as
$h_{ij}({\bf k})=h_{+}({\bf k})e_{ij}^{+}({\bf\hat{k}})+h_{\times}({\bf
k})e_{ij}^{\times}({\bf\hat{k}})$ where $h_{+}$ and $h_{\times}$ are gauge
independent components corresponding to two polarization modes. In this case
the polarization degree will be equal to ${\mathcal{P}}(k)={\langle
h^{\star}_{\times}({\mathbf{k}})h_{+}({\mathbf{k}^{\prime}})-h^{\star}_{+}({\mathbf{k}})h_{\times}({\mathbf{k}^{\prime}})\rangle}/{\langle
h^{\star}_{+}({\mathbf{k}})h_{+}({\mathbf{k}^{\prime}})+h^{\star}_{\times}({\mathbf{k}})h_{\times}({\mathbf{k}^{\prime}})\rangle}$

${\mathcal{P}}(k)=\frac{\langle
h^{\star}_{+}({\mathbf{k}})h_{+}({\mathbf{k}^{\prime}})-h^{\star}_{-}({\mathbf{k}})h_{-}({\mathbf{k}^{\prime}})\rangle}{\langle
h^{\star}_{+}({\mathbf{k}})h_{+}({\mathbf{k}^{\prime}})+h^{\star}_{-}({\mathbf{k}})h_{-}({\mathbf{k}^{\prime}})\rangle}=\frac{{\mathcal{H}}(k)}{H(k)},~{}$
(1)

where $H({k})$ and ${\mathcal{H}}({k})$ characterize the GW amplitude and
polarization (chirality) defined through the Gaussian-distributed GWs wave
number-space two-point function (for simplicity of notations we omit the time-
dependence):

$\frac{\langle
h^{\star}_{ij}({\mathbf{k}})h_{lm}({\mathbf{k}^{\prime}})\rangle}{(2\pi)^{3}}=\delta^{(3)}({\bf
k}-{\bf
k^{\prime}})\left[{\mathcal{M}}_{ijlm}H(k)+i{\mathcal{A}}_{ijlm}{\mathcal{H}}(k)\right],$

where $4{\mathcal{M}}_{ijlm}({\mathbf{\hat{k}}})\equiv
P_{il}P_{jm}+P_{im}P_{jl}-P_{ij}P_{lm}$, and
$8{\mathcal{A}}_{ijlm}({\mathbf{\hat{k}}})\equiv{\hat{\bf
k}}_{q}(P_{jm}\epsilon_{ilq}+P_{il}\epsilon_{jmq}+P_{im}\epsilon_{jlq}+P_{jl}\epsilon_{imq})$
are tensors,
$P_{ij}({\bf\hat{k}})\equiv\delta_{ij}-{\hat{k}}_{i}{\hat{k}}_{j}$ is the
projection operator, $\delta_{ij}$ and $\epsilon_{ijl}$ are the Kronecker
delta and the fully antisymmetric tensor, 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 Tevzadze:2012kk . 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.

Figure 3: Degree of circular polarization for (a) kinetically and (b)
magnetically forced cases with $\sigma=0$ (black) $0.1$ (blue), $0.3$ (green),
$0.5$ (orange), and $1$ (red). Approximate error bars based on the temporal
fluctuations and statistical spread for different random seeds of the forcing
are shown as solid black lines for $\sigma=0$ and as dotted lines otherwise.
The wave number is in units of the comoving Hubble frequency.  Figure 4:
$H(k)$ (left) and $|{\cal H}(k)|$ (right) for (a,b) kinetically and (c,d)
magnetically forced cases with $\sigma=0$ (Runs K0 and M0, black), $0.1$ (Runs
K01 and M01, blue), $0.3$ (Runs K03 and M03, green), $0.5$ (Runs K05 and M05,
orange), and $1$ (Runs K1 and M1, red).  Figure 5: (a) $2E_{K}(k)$ and (b)
$|H_{K}(k)/k|$ for kinetically forced cases, and (c) $2E_{M}(k)$, and (d)
$|kH_{M}(k)|$ for magnetically forced cases with $\sigma=0$ (Runs K0 and M0,
black), $0.1$ (Runs K01 and M01, blue), $0.3$ (Runs K03 and M03, green), $0.5$
(Runs K05 and M05, orange), and $1$ (Runs K1 and M1, red).

As we have highlighted above, previous works Kahniashvili:2005qi ;
Kisslinger:2015hua  considered stationary turbulence with two different models
for helical turbulence realization: (i) Kolmogorov-like helical turbulence
with two different spectral slopes, ${-5/3}$ and $-11/3$ for the spectral
energy and helicity densities, $E_{\rm M}(k)$ and $H_{\rm M}(k)$,
respectively, and (ii) helicity-transferring turbulence (if helicity transfer
and small-scale helicity dissipation dominates Kraichanan:1973 ) with the
spectral indices $-7/3$ and $-10/3$ for $E_{\rm M}(k)$ and $H_{\rm M}(k)$
Moiseev:1996 666Note that Ref. Kahniashvili:2005qi  used a different
convention for the spectral indices referring to the power spectra of the
symmetric $[P_{S}(k)\propto E_{\rm M}(k)/k^{2}]$ and the antisymmetric parts
$[P_{A}(k)\propto H_{\rm M}(k)/k]$.. The former case seems most suitable for
describing the usual nonhelical turbulence experiencing forward cascading
Borue:1997 . 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 Kahniashvili:2005qi  might be justified using Kolmogorov-like
turbulence at large wavenumbers and approximating the low wave number tail
with helicity transfer turbulence Ditlevsen:2001 . Following this description,
Ref. Kahniashvili:2008er  and later also Ref. Ellis:2020uid  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).

###  III.4 Dependence on the Driving

We now investigate the GW polarization spectra dependence on the nature of
driving. We show the polarization degree spectra in Fig. 3 for continuous
pumping of kinetic or magnetic energy and helicity at intermediate scales. The
parameters of these simulations are the same as those used for Fig. 2; see
also Ref. DATA . 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: $\propto 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,
which also explains the shallower slope in the polarization degree.

###  III.5 Comparison of the Spectra

In Fig. 4 we show the numerator and denominator of the degree of polarization,
$H(k)$ and ${\cal H}(k)$, respectively. The underlying kinetic and magnetic
energy and helicity spectra are shown in Fig. 5. The figure 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.

The departure from the theoretical predictions is due to the assumption of a
scale-independent 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 Kahniashvili:2005qi . 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
Niksa:2018ofa .

##  IV Conclusions

In this paper we present numerical simulations of the circular polarization
degree of GWs generated through parity violating (helical) 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 discrimination
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 in 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. In
particular, we have shown that the polarization spectra increase at low
frequencies due to the inverse cascade, reflecting redistribution of the
helical structures, while the energy density of the GWs sustains the scaling
that was established soon after the turbulent source activation. In fact,
previous works Kahniashvili:2005qi ; Ellis:2020uid  predicted a completely
different picture for both helical Kolmogorov turbulence and helicity-
transferring turbulence. More precisely, no previous work addressed the
dependence of the polarization spectra on the temporal characteristics of the
source, nor was the helicity transfer-induced increase of polarization spectra
at low frequencies detected. There was never a self-consistent description of
GW polarization nor an attempt to disentangle the nature of the underlying
source (kinetic vs magnetic driving). Based on our results, we argue that the
inverse cascade time scale determines the height of the polarization spectra
at low frequencies (a second peak) in the MHD case – even for low fractional
helicity. The second (smaller) peak is located at lower frequencies.
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 Caprini:2010xv ;
Neronov:2020qrl ), and/or any future planned missions, including atomic
interferometry Dimopoulos:2007cj .

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

###### Acknowledgements.

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.

## References

  * (1) M. Kamionkowski, A. Kosowsky and M. S. Turner, “Gravitational radiation from first order phase transitions,” Phys. Rev. D 49, 2837 (1994).
  * (2) E. Witten, “Cosmic Separation of Phases,” Phys. Rev. D 30, 272 (1984).
  * (3) C. J. Hogan, “Gravitational radiation from cosmological phase transitions,” MNRAS 218, 629 (1986).
  * (4) 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).
  * (5) M. Christensson, M. Hindmarsh and A. Brandenburg, “Inverse cascade in decaying 3-D magnetohydrodynamic turbulence,” Phys. Rev. E 64, 056405 (2001).
  * (6) T. Kahniashvili, A. Brandenburg, A. G. Tevzadze and B. Ratra, “Numerical simulations of the decay of primordial magnetic turbulence,” Phys. Rev. D 81, 123002 (2010).
  * (7) A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze and T. Vachaspati, “The dynamo effect in decaying helical turbulence,” Phys. Rev. Fluids. 4, 024608 (2019).
  * (8) A. Kosowsky, A. Mack and T. Kahniashvili, “Gravitational radiation from cosmological turbulence,” Phys. Rev. D 66, 024030 (2002).
  * (9) C. Caprini and R. Durrer, “Gravitational waves from stochastic relativistic sources: Primordial turbulence and magnetic fields,” Phys. Rev. D 74, 063521 (2006).
  * (10) G. Gogoberidze, T. Kahniashvili and A. Kosowsky, “The Spectrum of Gravitational Radiation from Primordial Turbulence,” Phys. Rev. D 76, 083002 (2007).
  * (11) T. Kahniashvili, L. Campanelli, G. Gogoberidze, Y. Maravin and B. Ratra, “Gravitational Radiation from Primordial Helical Inverse Cascade MHD Turbulence,” Phys. Rev. D 78, 123006 (2008).
  * (12) C. Caprini, R. Durrer and G. Servant, “The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition,” JCAP 0912, 024 (2009).
  * (13) R. Durrer and T. Kahniashvili, “CMB anisotropies caused by gravitational waves: A Parameter study,” Helv. Phys. Acta 71, 445 (1998).
  * (14) S. Weinberg, “Damping of tensor modes in cosmology,” Phys. Rev. D 69, 023503 (2004).
  * (15) D. V. Deryagin, D. Y. Grigoriev, V. A. Rubakov and M. V. Sazhin, “Possible Anisotropic Phases in the Early Universe and Gravitational Wave Background,” Mod. Phys. Lett. A 1, 593 (1986).
  * (16) C. Caprini, et al. “Detecting gravitational waves from cosmological phase transitions with LISA: an update,” JCAP 03, 024 (2020).
  * (17) A. G. Cohen, D. B. Kaplan and A. E. Nelson, “Baryogenesis at the weak phase transition,” Nucl. Phys. B 349, 727 (1991).
  * (18) A. G. Cohen, D. B. Kaplan and A. E. Nelson, “WEAK SCALE BARYOGENESIS,” Phys. Lett. B 245, 561 (1990).
  * (19) A. J. Long, E. Sabancilar and T. Vachaspati, “Leptogenesis and Primordial Magnetic Fields,” JCAP 1402, 036 (2014).
  * (20) G. C. Dorsch, S. J. Huber, T. Konstandin and J. M. No, “A Second Higgs Doublet in the Early Universe: Baryogenesis and Gravitational Waves,” JCAP 1705, 052 (2017).
  * (21) T. Kahniashvili, G. Gogoberidze and B. Ratra, “Polarized cosmological gravitational waves from primordial helical turbulence,” Phys. Rev. Lett. 95, 151301 (2005).
  * (22) T. Kahniashvili, G. Gogoberidze and B. Ratra, “Gravitational Radiation from Primordial Helical MHD Turbulence,” Phys. Rev. Lett. 100, 231301 (2008).
  * (23) L. Kisslinger and T. Kahniashvili, “Polarized Gravitational Waves from Cosmological Phase Transitions,” Phys. Rev. D 92, 043006 (2015).
  * (24) S. Alexander, E. McDonough and D. N. Spergel, “Chiral Gravitational Waves and Baryon Superfluid Dark Matter,” JCAP 1805, 003 (2018).
  * (25) S. Anand, J. R. Bhatt and A. K. Pandey, “Chiral plasma instability and primordial Gravitational wave,” Eur. Phys. J. C 79, 119 (2019).
  * (26) P. Niksa, M. Schlederer and G. Sigl, “Gravitational Waves produced by Compressible MHD Turbulence from Cosmological Phase Transitions,” Class. Quant. Grav. 35, 144001 (2018).
  * (27) J. Ellis, M. Fairbairn, M. Lewicki, V. Vaskonen and A. Wickens, “Detecting circular polarisation in the stochastic gravitational-wave background from a first-order cosmological phase transition,” JCAP 10, 032 (2020).
  * (28) S. H. S. Alexander, M. E. Peskin and M. M. Sheikh-Jabbari, “Leptogenesis from gravity waves in models of inflation,” Phys. Rev. Lett. 96, 081301 (2006).
  * (29) D. H. Lyth, C. Quimbay and Y. Rodriguez, “Leptogenesis and tensor polarisation from a gravitational Chern-Simons term,” JHEP 03, 016 (2005).
  * (30) P. Adshead, A. J. Long and E. I. Sfakianakis, “Gravitational Leptogenesis, Reheating, and Models of Neutrino Mass,” Phys. Rev. D 97, 043511 (2018).
  * (31) H. Abedi, M. Ahmadvand and S. S. Gousheh, “Electroweak baryogenesis via chiral gravitational waves,” Phys. Lett. B 786, 35 (2018).
  * (32) T. Vachaspati, “Progress on Cosmological Magnetic Fields,” [arXiv:2010.10525 [astro-ph.CO]].
  * (33) N. Seto, “Prospects for direct detection of circular polarization of gravitational-wave background,” Phys. Rev. Lett. 97, 151101 (2006).
  * (34) T. L. Smith and R. Caldwell, “Sensitivity to a Frequency-Dependent Circular Polarization in an Isotropic Stochastic Gravitational Wave Background,” Phys. Rev. D 95, 044036 (2017).
  * (35) K. W. Masui, U. L. Pen and N. Turok, “Two- and Three-Dimensional Probes of Parity in Primordial Gravity Waves,” Phys. Rev. Lett. 118, 221301 (2017).
  * (36) N. Seto and A. Taruya, “Measuring a Parity Violation Signature in the Early Universe via Ground-based Laser Interferometers,” Phys. Rev. Lett. 99, 121101 (2007).
  * (37) N. Seto and A. Taruya, “Polarization analysis of gravitational-wave backgrounds from the correlation signals of ground-based interferometers: Measuring a circular-polarization mode,” Phys. Rev. D 77, 103001 (2008).
  * (38) N. Seto, “Quest for circular polarization of gravitational wave background and orbits of laser interferometers in space,” Phys. Rev. D 75, 061302 (2007).
  * (39) S. G. Crowder, R. Namba, V. Mandic, S. Mukohyama and M. Peloso, “Measurement of Parity Violation in the Early Universe using Gravitational-wave Detectors,” Phys. Lett. B 726, 66 (2013).
  * (40) V. Domcke, et al. J. Garcia-Bellido, M. Peloso, M. Pieroni, A. Ricciardone, L. Sorbo and G. Tasinato, “Measuring the net circular polarization of the stochastic gravitational wave background with interferometers,” JCAP 05, 028 (2020).
  * (41) 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).
  * (42) The Pencil Code, https://github.com/pencil-code, DOI:10.5281/zenodo.2315093
  * (43) J. D. Romano and N. J. Cornish, “Detection methods for stochastic gravitational-wave backgrounds: a unified treatment,” Living Rev. Rel. 20, 2 (2017).
  * (44) 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 Dyn. 114, 130 (2020).
  * (45) T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal, and A. Roper Pol, Datasets for “Circular Polarization of Gravitational Waves from Early-Universe Helical Turbulence” (v2020.11.07), https://doi.org/10.5281/zenodo.4256906; see also http://www.nordita.org/~brandenb/projects/CircPol/ for easier access.
  * (46) A. G. Tevzadze, L. Kisslinger, A. Brandenburg and T. Kahniashvili, “Magnetic Fields from QCD Phase Transitions,” Astrophys. J. 759, 54 (2012).
  * (47) R. H. Kraichnan, ”Helical turbulence and absolute equilibrium”, J. Fluid Mech., 59, 745 (1973).
  * (48) S. S. Moiseev and O. Chkhetiani, ”Helical scaling in turbulence”, JETP 83, 192, (1996).
  * (49) V. Borue and S. A. Orszag, ”Spectra in helical three-dimensional homogeneous isotropic turbulence” Phys. Rev. E 55, 7005 (1997).
  * (50) P. D. Ditlevsen and P. Giuliani, ”Cascades in helical turbulence”, Phys. Rev. E 63, 036304 (2001).
  * (51) 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).
  * (52) 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]].
  * (53) S. Dimopoulos, P. W. Graham, J. M. Hogan, M. A. Kasevich and S. Rajendran, “Gravitational Wave Detection with Atom Interferometry,” Phys. Lett. B 678, 37 (2009).