1 Introduction

Gravitational wave (GW) once generated, propagates almost unhindered through the space-time. This property makes GW a very powerful probe of the source which produces it as well as the medium through which it propagates (see [1,2,3] and references therein). From the cosmological point of view, the most interesting gravitational radiation is that of the stochastic gravitational wave (SGW) background. Such gravitational radiations are produced by events in the early stages of the Universe and hence, may decipher the physics of those epochs. Several attempts have already been made in this regard and various sources of SGW have been considered. List of SGW source includes quantum fluctuations during inflation [4,5,6,7], bubble wall collision during phase transition [8,9,10,11,12,13], cosmological magnetic fields [14,15,16,17] and turbulence in the plasma [16, 18, 19].

In the early universe, before the electroweak phase transition, many interesting phenomena have taken place. For instance, it has been shown by several authors [20,21,22,23,24,25] that in presence of asymmetry in the left-handed and the right-handed particles in the early Universe, there will be an instability which leads to the generation of turbulence in the plasma as well as (hyper-charge) magnetic fields. In reference [26], it has been shown that these magnetic fields can be generated even in absence of net chiral charge but due to the gravitational anomaly. The magnetic field generated via this mechanism are helical in nature. However, helicity (\({\mathscr {H}}_B =\frac{1}{V} \int d^3 x\, \, \mathbf{A} \cdot \mathbf{B}\)) of these magnetic fields are not completely conserved due to the fact that large but finite conductivity gives a slight time variation of the helicity density. In presence of chiral imbalance, the chiral charge conservation equation, which is valid at temperature \(T>T_R \approx 80\) TeV [20], is given as: \(\partial _\eta (\varDelta \mu + \frac{\alpha '}{\pi } {\mathscr {H}}_B) = 0\) where \(\varDelta \mu =\mu _R - \mu _L\), \(\alpha '\) and \(\eta \) are the asymmetry in the chiral chemical potentials and \(U(1)_Y\) fine structure constant and conformal time respectively. At the onset of the instability, \({\mathscr {H}}_B \approx 0\) and at subsequent time, helicity will grow at the expense of chiral chemical potential. In this regime, \(\varDelta \mu \) can be regarded as constant. On the other hand, at temperature \(T< 80\) TeV the above conservation equation is not valid as chiral flipping \(\varGamma _{f}\) rate is non-vanishing and the hence, the conservation equation is given by [22]

$$\begin{aligned} \frac{d}{d \eta }\left( \varDelta \mu + \frac{\alpha '}{\pi } {\mathscr {H}}_B\right) =-\varGamma _f \varDelta \mu . \end{aligned}$$
(1)

Here \( \varGamma _f=\frac{T_R}{M_*}T\) and \(M_* = \left( {90\over 8\pi ^3g_{\mathrm{eff}}}\right) ^{1/2}M_{pl}\,\), where \(g_{\mathrm{eff}}\) and \(M_{pl}=1/\sqrt{G}\) are relativistic degree of freedom and the Planck mass respectively [20]. In this regime, non-linearity sets in and the magnetic fields are generated. The generated magnetic fields show inverse cascade behaviour, where magnetic energy is transfered from small scale to large-scale. In reference [27], it is shown that the currently observed baryon asymmetry (\(\eta _B \sim 10^{-10}\)) can be generated if the magnetic fields produced above electroweak scale undergoes the inverse cascade and the strength of the magnetic field is of the order \((10^{-14}-10^{-12})\) G at the galactic scale. The generated magnetic fields induce a anisotropic stress so that their energy density \(B^2/8\pi \) must be a small perturbation, in order to preserve the isotropy of the Friedmann-Robertson-Walker background. This condition allows us to use the linear perturbation theory. In this perturbation scheme, peculiar velocity and magnetic fields are considered to be first order in perturbations. At sufficiently large length scale, the effect of the fluid on the evolution of magnetic fields can be neglected. However, at small length scales the interaction between the fluid and the magnetic field become very crucial. At the intermediate length scale, plasma undergoes Alfven oscillations and on a very small scale (viscous scale), these fields undergo exponential damping due to the shear viscosity [28]. Thus, the large-scale magnetic fields are important for the physics at the cosmic scale.

We have already mentioned that the seed magnetic field can be generated even in absence of net chiral charge but due to gravitational anomaly [26] and the magnetic field thus generated can produce instability in the plasma. These magnetic fields contribute a anisotropic stress to the energy-momentum tensor and hence can act as a source for the generation of the GWs. The underlying physics of GW generation is completely different from previously considered scenarios. Therefore, it is important to investigate the generation and evolution of GW in this context. In this article, we compute the metric tensor perturbation due to the chiral magnetic field. Since chiral magnetic field, which sources the tensor perturbations, has a unique spectrum, the GWs generated is expected to have a unique signature in its spectrum as well. Moreover, we compute the amplitude and frequency of the GW and show its dependences on the model parameters. Consequently, any detection of SGW in future measurements like eLISA will constrain or rule out such theoretical constructs.

This paper is organized as follows: in Sect. 2 we outline the generation and evolution of magnetic field due to gravitational anomaly and chiral imbalance. We discuss the generation of SGW in Sect. 3. We present our results in Sect. 4 and finally conclude in (5). Throughout this work, we have used \(\hslash =c=k_B=1\) unit. We have also considered Friedman-Robertson-Walker metric for expanding background space-time

$$\begin{aligned} ds^2=a^2(\eta )\left( -d\eta ^2+\delta _{ij}\,dx^i\,dx^j\right) \, , \end{aligned}$$
(2)

where scale factor \(a(\eta )\) have dimension of length, whereas conformal time \(\eta \) and conformal coordinate \(x^i\) are dimensionless quantities. In the radiation dominated epoch \(a=1/T\), we can define conformal time \(\eta = M_*/T\). Unless stated otherwise, we will work in terms of comoving variables defined as,

$$\begin{aligned} B_c = a^2(\eta )B(\eta ),~~ \mu _c = a(\eta )\mu ,~~ k_c =a\, k, \end{aligned}$$
(3)

where, B, \(\mu \) and k represents the physical magnetic fields, chemical potential, and the wave number respectively. It is clear from the convention used here that all the comoving quantities are dimensionless. In terms of these comoving variables, the evolution equations of fluid and electromagnetic fields are form invariant [29,30,31]. Therefore, we will work with the above defined comoving quantities and omit the subscript “c” in our further discussion.

2 Gravitational anomaly and magnetic fields in the early universe

Although the origin of large-scale magnetic field is still an unsolved issue in cosmology, several attempts have been made to address the issue. It has been discussed in the literature that there are processes in the early universe, much above electroweak scale, which can lead to more number of right-handed particle than the left-handed ones and remains in thermal equilibrium via its coupling with the hypercharge gauge bosons [32, 33]. Furthermore, if the plasma has rotational flow or external gauge field present, there could be a current in the direction parallel to the vorticity due to rotational flow or parallel to the external field. The current parallel to the vorticity is known as chiral vortical current and the phenomenon is called as chiral vortical effects (CVE) [34,35,36,37,38]. Similarly, the current parallel to the external magnetic field is known as chiral magnetic current and the phenomenon is called as chiral magnetic effect (CME) [39,40,41,42,43]. CVE and CME are characterized by the transport coefficients \(\xi \) and \(\xi ^{(B)}\) respectively. The form of these coefficients can be obtained by demanding the consistency with the second law of thermodynamics (\(\partial _\mu \,s^\mu \ge 0\), with \(s^\mu \) being the entropy density). Thus, in presence of chiral imbalance and gravitational anomaly, which arises due to the coupling of spin with gravity [44], the coefficients for each right and left particle have the following form [37, 38, 41]

$$\begin{aligned} \xi _i= & {} C\,\mu _i^2\,\left[ 1-\frac{2\,n_i\,\mu _i}{3\,(\rho + p)}\right] \, + \frac{D\,T^2}{2}\left[ 1-\frac{2\,n_i\,\mu _i}{(\rho + p)}\right] \, , \nonumber \\\end{aligned}$$
(4)
$$\begin{aligned} \xi _i^{(B)}= & {} C\,\mu _i\,\left[ 1-\frac{n_i\,\mu _i}{2\,(\rho + p)}\right] \, - \frac{D}{2}\left[ \frac{n_i\,T^2}{(\rho + p)}\right] \, . \end{aligned}$$
(5)

In above equations ‘i’ stands for each species of the chiral plasma. The constants C and D are related to those of the chiral anomaly and mixed gauge-gravitational anomaly and are given as \(C=\pm 1/4\pi ^2\) and \(D=\pm 1/12\) for right and left-handed chiral particles respectively. The variables n, \(\rho \) and p are respectively the number density, energy density and pressure density.

Using the effective Lagrangian for the standard model, one can derive the generalized Maxwell’s Eq. [45],

$$\begin{aligned} \mathbf {\nabla }\times {\mathbf {B}}= {\mathbf {j}}\, , \end{aligned}$$
(6)

where \({\mathbf {j}}\) is defined as \({\mathbf {j}} = {{\mathbf {j}}}_{\mathrm{v}} + {{\mathbf {j}}}_5\) with \({{\mathbf {j}}}_{\mathrm{v}}\) being the vector current and \({{\mathbf {j}}}_5 \) is the axial current. Vector and axial currents respectively takes the following form:

$$\begin{aligned} j^\mu _{\mathrm{v}}= & {} n_{\mathrm{v}}\,u^\mu + \sigma E^\mu + \xi _\mathrm{v}\,\omega ^\mu + \xi _{\mathrm{v}}^{(B)}\,B^\mu \, , \end{aligned}$$
(7)
$$\begin{aligned} j_5^\mu= & {} n_5\,u^\mu + \xi _{\mathrm{5}}\,\omega ^\mu + \xi _{5}^{(B)}B^\mu \, . \end{aligned}$$
(8)

In the above equations, any quantity \(x_{\mathrm{v, (5)}}\) denotes the sum (difference) of the quantities pertaining to right and left handed particles. Also \(E^\mu =u_{\nu } F^{\mu \nu }\), \(B^\mu =1/2 \varepsilon ^{\mu \nu \sigma \delta } u_\nu F_{\sigma \delta }\), and \(\omega ^\mu = 1/2 \varepsilon ^{\mu \nu \sigma \delta } u_\nu \partial _\sigma u_\delta \) are the electric, magnetic and the vorticity four vectors respectively. We have ignored the displacement current in Eq. (6). Taking \(u^\mu = (1, \mathbf{v})\) and using Eqs. (7)– (8), one can show that

$$\begin{aligned} j^{0}= & {} n = n_{\mathrm{v}} + n_5 \end{aligned}$$
(9)
$$\begin{aligned} {\mathbf {j}}= & {} n{\mathbf {v}} + {\sigma (\mathbf {E}} +{\mathbf {v}}\times {\mathbf {B}}) + \xi \, {\mathbf {\omega }} + \xi ^{(B)}\,{\mathbf {B}} \, , \end{aligned}$$
(10)

with \(\xi = \xi _{\mathrm{v}} + \xi _5\) and \(\xi ^{(B)} = \xi ^{(B)}_{\mathrm{v}} + \xi ^{(B)}_5\). Assuming the velocity field to be divergence free field, i.e. \(\mathbf {\nabla }\cdot \mathbf {v} = 0\) and taking curl of Eq. (6) along with the expression for current from Eq. (10) we obtain,

$$\begin{aligned} \frac{\partial {\mathbf {B}}}{\partial \eta }= & {} \frac{n }{\sigma }{\varvec{\omega }} +\frac{1}{\sigma }\nabla ^2{\mathbf {B}}+ \mathbf {\nabla }\times ({\mathbf {v}}\times {\mathbf {B}}) \nonumber \\&+ \frac{\xi }{\sigma }\mathbf {\nabla }\times {\varvec{\omega }} + \frac{\xi ^{(B)}}{\sigma }\mathbf {\nabla }\times {\mathbf {B}}. \end{aligned}$$
(11)

In our previous work [26], we discussed that the seed magnetic field (for which \({\mathbf {B}}\) in the right hand side of Eq. (11) is zero) can be generated even if \(n = 0\). The \(T^2\) term in \(\xi _i\) [see Eq. (4)], which arises due to the gravitational anomaly, acts as a source for the generation of seed field. On the other hand, presence of finite chiral imbalance such that \(\mu /T \ll 1\), \(T^2\) term in \(\xi \) still acts as source of seed magnetic field but non-zero \(\xi ^{(B)}\) triggers instability in the system. This result is in agreement with the previous studies where it was shown that in the presence of a chiral imbalance in the plasma, much above the Electroweak scale (\(T>100~\mathrm{GeV}\)), there can be instability known as chiral plasma instability [46].

The production and evolution of the magnetic field can be seen through the evolution equation given in Eq. (11). In order to do so, we decompose the divergence-free vector fields, e.g. magnetic field, in the orthonormal helicity basis, \(\varepsilon ^{\pm }_i\), defined as

$$\begin{aligned} \varvec{ \varepsilon }^\pm (\mathbf{k}) = \frac{-i}{\sqrt{2}}\left[ \mathbf{e}_1(\mathbf{k})\,\pm \,i\, \mathbf{e}_2(\mathbf{k})\right] \exp (i \, \mathbf{k} \cdot \mathbf{x}), \end{aligned}$$
(12)

where \((\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3 = \hat{\mathbf{k}})\) form a right-handed orthonormal basis with \(\mathbf{e}_2 = {\hat{\mathbf{k}}} \times \mathbf{e}_1\). We choose \(\mathbf{e_1} \) to remain invariant under the transformation \(\mathbf{k}\rightarrow -\mathbf{k}\) while \(\mathbf{e_2} \) flip its sign. In this basis, the magnetic field can be decomposed as

$$\begin{aligned} B_i(\eta , \mathbf{k})= B^+(\eta , \mathbf{k})\,\varepsilon _i^+(\mathbf{k})\, + \, B^-(\eta , \mathbf{k})\, \varepsilon _i^-(\mathbf{k})\, . \end{aligned}$$
(13)

In this basis, the evolution equation for the \(|B^{\pm }|^2\) can be obtained from Eq. (11)

$$\begin{aligned} \frac{\partial |{\tilde{B}}^\pm |^2}{\partial \eta }= & {} \frac{2}{\sigma }\left( -k^2\pm \xi ^{(B)} k\right) |{\tilde{B}}^\pm |^2 + \frac{2}{\sigma ^2}\left( \pm nk + \xi k^2\right) ^2 \nonumber \\&\quad \times |\mathrm{v}^\pm |^2{{\mathscr {F}}}_s \end{aligned}$$
(14)

where \({{\mathscr {F}}}_s=\eta -\eta _0\) for \(\eta -\eta _0\le 2\pi /(k\mathrm{v})\) and zero for \(\eta -\eta _0\ge 2\pi /(k\mathrm{v})\). Also magnitude of the wave vector \(\mathbf{k}\) is represented as k, i.e. \(|\mathbf{k}|=k\).  From Eq. (14), it is clear that, when first term is dominant over the second term, magnetic modes will grow exponentially with time as

$$\begin{aligned} |B^\pm |^2 = |B_0|^2 \exp \left( \frac{2\eta }{\sigma }k_\mathrm{ins}^2\right) \exp \left( -\frac{2\eta }{\sigma }(k\mp k_{\mathrm{ins}})^2\right) , \end{aligned}$$

where \(k_{\mathrm{ins}}= \alpha ' \varDelta \mu /\pi \) is the value of the k at which magnetic modes have maximum growth rate. The exponential growth of the magnetic modes is true only in the linear regime. In this regime, chiral chemical potential remain constant and magnetic fields are generated at the cost of chiral imbalance. However, when generated magnetic field is sufficiently large, the non-linear effects become prominent. In this case, we need to consider the evolution of the chemical potential [20, 23], which is given by Eq. (1).

At temperature \(T> 80\) TeV, perturbative processes that lead to the flipping of the chirality are small compared to the expansion rate of the Universe and hence we can consider the chiral symmetry to be an exact symmetry of the theory. Therefore, a chemical potential \(\mu _{R, L}\) for each species can be defined. However, at \(T< 80\) TeV, the chiral symmetry is not exact due to the chiral flipping of the right-handed particles to the left-handed particles and vice versa. As a consequence, the number densities of these particle are not conserved and therefore, we can not define the chemical potentials for right and left handed particles [22, 47]. In order to obtain the velocity profile, we used the scaling symmetry [48, 49] rather than solving the Navier-Stokes equation. In our earlier work [26], we have obtained \(|\mathrm{v}^+|=|\mathrm{v}^-| = \pi k^{-3/2}\mathrm{v}\). Here \(\mathrm{v}\) is given by

$$\begin{aligned} \mathrm{v}(k, \eta )= \mathrm{v}_i \left( \frac{k}{k_i(\eta _0)}\right) ^{(1+m)/2}\left( \frac{\eta }{\eta _0}\right) ^{(m-2)/6}. \end{aligned}$$
(15)

In above equation, m is a positive integer and \(\mathrm{v}_i\), \(k_i\) are arbitrary function encoding the boundary conditions. We also showed that the scaling law allows more power in the magnetic field. For \(n = 0\), \(E_B\propto k^7\) at larger length scale whereas for \(n\ne 0\), \(E_B\propto k^5\) instead of \(k^7\) [23]. However, in the both scenarios \(E_B\) is more than that of the case without considering the scaling symmetry [23]. This aspect is important for this analysis as more power in the helical magnetic field can generate larger anisotropic stress.

The two point correlation function for the helical magnetic fields is given as:

$$\begin{aligned} \langle B_i({{\varvec{k}})}B^*_j({{\varvec{k}}'})\rangle= & {} {(2\pi )^3 \over 2} \delta ({{\varvec{k}}-{\varvec{k}}'}) \left[ P_{ij}S(k)\right. \nonumber \\&\left. +i\,\epsilon _{ijl}{\hat{k}}_l{\mathscr {A}}(k)\right] ~~ \end{aligned}$$
(16)

where S(k) is the symmetric and \({\mathscr {A}}(k)\) is the helical part of the magnetic field power spectrum. \(P_{ij} = \delta _{ij} - {\hat{k}}_i\,{\hat{k}}_j\) is the transverse plane projector which satisfies: \(P_{ij}{\hat{k}}_j = 0 \), \(P_{ij}P_{jk} = P_{ik}\) with \({\hat{k}}_i = k_i/k\) and \(\epsilon _{ijk}\) is the totally antisymmetric tensor. Using Eq. (16) and the reality condition \(B^{\pm *}({\varvec{k}})=-B({\varvec{k}})^{\pm }\), one can show that

$$\begin{aligned}&\langle B^+({\varvec{k}})B^{+*}({\varvec{k}'}) + B^-({\varvec{k}})B^{-*}({\varvec{k}'})\rangle \nonumber \\&\quad = -(2\,\pi )^3 S(k)\delta ({\varvec{k} -\varvec{k}'}) \end{aligned}$$
(17)
$$\begin{aligned}&\langle B^+({\varvec{k}})B^{+*}({\varvec{k}'}) - B^-({\varvec{k}})B^{-*}({\varvec{k}'})\rangle \nonumber \\&\quad = (2\,\pi )^3\,{\mathscr {A}}(k)\delta ({\varvec{k} - \varvec{k}'}). \end{aligned}$$
(18)

Note that, \({\mathscr {A}}({{\varvec{k}}})\) represents the difference in the power of left-handed and right-handed magnetic fields, however a maximally helical magnetic fields configuration can be achieved when \({\mathscr {A}}({{\varvec{k}}})=S({{\varvec{k}}})\).

2.1 Anisotropic stress

Tensor component of metric perturbation, which results in the gravitational waves, are sourced by the transverse-traceless part of the stress-energy tensor. In this work, we assume that the anisotropic stress is generated by the magnetic stress-energy tensor which is given by

$$\begin{aligned} T_{ij}= \frac{1}{a^2}\left( B_i\,B_j-\frac{1}{2}\delta _{ij} \, B^2\right) \, , \end{aligned}$$
(19)

Note that the spatial indices are raised, lowered and contracted by the Kronecker delta such that \(B^2 = \delta ^{ij}B_iB_j\). The magnetic field component \(B_i\) then coincide with the comoving magnetic field which in our notation is \(B_c = a^2\, B\). In Fourier space, the stress-energy tensor for the magnetic field take the following form

$$\begin{aligned} T_{ij}({\varvec{k}})= & {} \frac{a^{-2}}{2(2\pi )^4}\int d^3q \left[ B_i(\varvec{q})B^*_j({\varvec{q}}-{\varvec{k}}) \right. \nonumber \\&\left. - \frac{1}{2}B_l({\varvec{q}})B_l^*({\varvec{q}}-{\varvec{k}})\delta _{ij} \right] . \end{aligned}$$
(20)

We are interested in the generation of GW and hence, we need to extract the transverse traceless component of the stress energy tensor given in Eq. (20). This can be done by using the projection operator: \({\mathscr {P}}_{iljm}({{\varvec{k}}})=[P_{il}({{\varvec{k}}})\, {P}_{ jm}({{\varvec{k}}})-\frac{1}{2} \,P_{ij}({{\varvec{k}}})\,P_{lm}({{\varvec{k}}})]\) which leads to [50]

$$\begin{aligned} \varPi _{ij}({\varvec{k}}) = {\mathscr {P}}_{iljm}({{\varvec{k}}})\, T_{lm}(k) . \end{aligned}$$
(21)

At this stage, we will evaluate the two-point correlation function of the energy-momentum tensor which will be used in the later part of the calculation. The two point correlation of the stress-energy tensor takes the form

$$\begin{aligned} \begin{aligned} \langle T_{ij}({\varvec{k}})T_{lm}^*({\varvec{k}}') \rangle&= {1 \over (2\pi )^6}{1\over (4\pi )^2 a^4}\int d^3 p\int d^3 q \\&\quad \times \bigg [\langle B_i(\varvec{p})\,B_j^*(\varvec{p}-{\varvec{k}})\,B_l^*(\varvec{q}) B_m(\varvec{q} -{\varvec{k}}') \rangle \\&\quad + (..)\delta _{ij} + (..)\, \delta _{lm} + (..)\delta _{ij} \delta _{lm}\bigg ] \end{aligned} \end{aligned}$$
(22)

It was shown in Ref. [51] that only first term in the angular bracket will have a non-vanishing contribution in the two-point correlation function of the anisotropic stress \(\langle \varPi _{ij}\varPi ^*_{lm}\rangle \). Therefore, we will not evaluate other terms. Moreover, around the chiral instability, the magnetic field profile is Gaussian and the major contribution to the anisotropic stress come from this regime only. Therefore, we can safely assume that the magnetic fields are Gaussian and hence four-point correlator in the integrand can be expressed, using Wick’s theorem, in terms of two-point correlators as

$$\begin{aligned} \begin{aligned}&\langle B_i({\varvec{k}}_i)B_j({\varvec{k}}_j)B_l({\varvec{k}}_l)B_m({\varvec{k}}_m)\rangle \\&\quad =\langle B_i({\varvec{k}}_i) B_j({\varvec{k}}_j)\rangle \,\langle B_l({\varvec{k}}_l) B_m({\varvec{k}}_m)\rangle \\&\qquad + \langle B_i({\varvec{k}}_i) B_l({\varvec{k}}_l)\rangle \,\langle B_j({\varvec{k}}_j) B_m({\varvec{k}}_m)\rangle \\&\qquad + \langle B_i({\varvec{k}}_i) B_m({\varvec{k}}_m)\rangle \,\langle B_j({\varvec{k}}_j) B_l({\varvec{k}}_l)\rangle \,. \end{aligned} \end{aligned}$$
(23)

After a bit of lengthy but straightforward calculation, we obtain the two point correlations of the energy momentum tensor and the transverse-traceless part of the energy momentum tensor is

$$\begin{aligned} \begin{aligned}&\langle \mathrm{T}_{\mathrm{ij}}(\varvec{k})\mathrm{T}_{\mathrm{lm}}^*(\varvec{k}') \rangle \\\quad&={1\over 4 (4\pi )^2 a^4}\delta ({\varvec{k}}-{\varvec{k}}')\int d^3 p \bigg [ S(p)S(k-p) \\&\quad \left\{ P_{il}(\hat{\varvec{p}})\, P_{jm}\widehat{(\varvec{k} - \varvec{p})} + P_{im}(\hat{\varvec{p}})\, P_{jl}\widehat{(\varvec{k} - \varvec{p})}\right\} \\&\qquad - {\mathscr {H}}(p)\,{\mathscr {H}}(k-p)\big \{\epsilon _{ila}\,\epsilon _{jmb }\,\hat{\varvec{p}}_a\,\widehat{({\varvec{k}} - \varvec{p})}_b \\&\qquad + \epsilon _{imc}\,\epsilon _{jld }\,\hat{\varvec{p}}_c\,\widehat{({\varvec{k}}- \varvec{p})}_d\big \} \!+ i\, {\mathscr {H}}(p)S(k-p) \\&\qquad \big \{\epsilon _{ila}\,P_{jm}(\widehat{\varvec{k} - \varvec{p}})\,\hat{\varvec{p}}_a\, + \epsilon _{imc}\,P_{jl}(\widehat{\varvec{k} - \varvec{p}})\,\hat{\varvec{p}}_c\,\big \} \\&\qquad + i\, {\mathscr {H}}(k-p)S(p)\big \{\epsilon _{jmb}\,P_{il}(\hat{\varvec{p}})\, (\widehat{\varvec{k} - \varvec{p}})_b\\&\qquad + \epsilon _{jld}P_{im}(\hat{\varvec{p}})\,(\widehat{\varvec{k} - \varvec{p}})_d\,\big \} \bigg ] \end{aligned} \end{aligned}$$
(24)

and

$$\begin{aligned} \langle \varPi _{\mathrm{ab}}({{\varvec{k}}})\,\varPi _{\mathrm{cd}}^*({{\varvec{k}}}') \rangle = {\mathscr {P}}_{\mathrm{aibj}}({{\varvec{k}}})\,{\mathscr {P}}_{\mathrm{cldm}}({{\varvec{k}}}')\langle \mathrm{T}_{\mathrm{ij}}({{\varvec{k}}})\,\mathrm{T}_{\mathrm{lm}}^*({{\varvec{k}}}') \rangle . \end{aligned}$$
(25)

Above equation can also be written in terms of a most general isotropic transverse traceless fourth rank tensor which obeys \(\mathcal{P}_{abcd}={\mathscr {P}}_{bacd}={\mathscr {P}}_{abdc}={\mathscr {P}}_{cdab}\) as [28, 52]

$$\begin{aligned} \langle \varPi _{\mathrm{ab}}({\varvec{k}})\varPi _{\mathrm{cd}}^*({\varvec{k}}') \rangle \!=\! \frac{1}{4a^4}[{\mathscr {M}}_{abcd}f(k)+i{\mathscr {A}}_{abcd} g(k)]\delta ({\varvec{k}}-{\varvec{k}}'), \end{aligned}$$
(26)

with a definition of \({\mathscr {M}}_{abcd}\) and \({\mathscr {A}}_{abcd}\) as

$$\begin{aligned} {\mathscr {M}}_{abcd}= & {} P_{ac}P_{bd}+P_{ad}P_{bc}-P_{ab}P_{cd} \,\nonumber \\ {\mathscr {A}}_{abcd}= & {} \frac{1}{2}\hat{{\varvec{k}}}_e (P_{bd}\epsilon _{ace}+P_{ac}\epsilon _{bde}+P_{ad}\epsilon _{bce} +P_{bc}\epsilon _{ade})\nonumber \\ \end{aligned}$$
(27)

which follows following properties:

$$\begin{aligned}&{\mathscr {M}}_{abcd} = {\mathscr {M}}_{cdab}= {\mathscr {M}}_{abdc}={\mathscr {M}}_{bacd}\nonumber \,\\&{\mathscr {A}}_{abcd} = {\mathscr {A}}_{cdab}= -{\mathscr {A}}_{abdc}=-{\mathscr {A}}_{bacd}\nonumber \,\\&{\mathscr {M}}_{abab} = 4\nonumber \,\\&{\mathscr {M}}_{aacd} = {\mathscr {M}}_{abcc} = 0\nonumber \,\\&P_{ea}{\mathscr {M}}_{abcd} = {\mathscr {M}}_{ebcd} \nonumber \, \\&P_{ea}{\mathscr {A}}_{abcd} = {\mathscr {A}}_{ebcd} \nonumber \,\\&{\mathscr {M}}_{abcd}\,{\mathscr {M}}_{abcd} = {\mathscr {M}}_{abcd}\,{\mathscr {M}}_{abcd} =8\nonumber \, \\&{\mathscr {A}}_{abcd}\,{\mathscr {M}}_{abcd} = 0 \nonumber \,\\&{\mathscr {A}}_{abab} = {\mathscr {A}}_{aacd} = {\mathscr {A}}_{abcc}=0. \end{aligned}$$
(28)

The functions f(k) and g(k) are defined as follows:

$$\begin{aligned} \delta ({{\varvec{k}}}-{{\varvec{k}}}')\, f(k)= & {} \frac{1}{2}{\mathscr {M}}_{abcd}\langle T_{ab}({{\varvec{k}}})\,T_{cd}^*({{\varvec{k}}}')\rangle \nonumber \, \\ \delta ({{\varvec{k}}}-{{\varvec{k}}}')\, g(k)= & {} -\frac{i}{2}{\mathscr {A}}_{abcd}\langle T_{ab}({{\varvec{k}}})\,T_{cd}^*({{\varvec{k}}}')\rangle . \end{aligned}$$
(29)

We point out that the functions f(k) and g(k) also depends on time because of the time dependence of the magnetic fields. The integral form of these two functions are

$$\begin{aligned} f(k)= & {} \frac{1}{4}\frac{1}{(4\pi )^2} \int d^3p\,\big [(1+\gamma ^2)(1+\beta ^2)\,S(p)S(k-p) \nonumber \\&+\, 4\,\gamma \,\beta \,{\mathscr {H}}(p){\mathscr {H}}(k-p) \, \big ] \end{aligned}$$
(30)
$$\begin{aligned} g(k)= & {} \frac{1}{2}\frac{1}{(4\pi )^2}\int d^3p\, \left[ (1+\gamma ^2)\,\beta \,S(p)\, {\mathscr {H}}(k-p)\right] .\nonumber \\ \end{aligned}$$
(31)

where \(p=|\mathbf{p}|\), \((k-p)=|{{\varvec{k}}}-\mathbf{p}|\), \(\gamma = \hat{{{\varvec{k}}}}\cdot \hat{\mathbf{p}}\) and \(\beta = \hat{{{\varvec{k}}}}\cdot (\widehat{{{\varvec{k}}}-\mathbf{p}})= (k-p\gamma )/\sqrt{k^2+p^2-2\gamma p k}\).

3 Gravitational waves from chiral magnetic fields

We have seen that the chiral magnetic field generated at very high temperature can produce anisotropic stress which leads to tensor perturbation in the metric. To linear order, the small tensor perturbation in the FLRW background can be written as:

$$\begin{aligned} ds^2= a^2(\eta )[-d\eta ^2+(\delta _{ij}+2 h_{ij})dx^i\, dx^j], \end{aligned}$$
(32)

where the tensor perturbation satisfies the following conditions \(h_i^i=0\) and \(\partial ^i \, h^j_i=0\). In this gauge, these tensor perturbations describe the GW whose evolution equation can be obtained by solving Einstein’s equation which, to the linear order in \(h_{ij}\), is given as:

$$\begin{aligned} h_{ij}'' + 2\,H\, h_{ij}' \, +k^2\, h_{ij}= 8\,\pi \, G \, \varPi _{ij}, \end{aligned}$$
(33)

where prime denotes the derivative with respect to the conformal time \(\eta \) and \(H =\frac{1}{a(\eta )} \frac{\partial a(\eta )}{\partial \eta }\). The time dependence in the right hand side of the Eq. (33) comes from the fact that the magnetic field is frozen in the plasma. Therefore, \(\varPi _{ij}(k,\eta )\) is a coherent source, in the sense that each mode undergoes the same time evolution. Assuming that the tensor perturbations has a Gaussian distribution function, the statistical properties can be described by the two-point correlation function given as,

$$\begin{aligned} \begin{aligned} \langle h_{ij}^{*'}({\varvec{k}}, \eta )\, h_{lm}'({\varvec{k}}', \eta )\rangle&= \frac{1}{4} \delta ^3({\varvec{k}}-{\varvec{k}}')\big [{\mathscr {M}}_{ijlm} S_{GW}(k,\eta ) \\&\quad + i\, {\mathscr {A}}_{ijlm}\,{\mathscr {H}}_{GW}(k,\eta ) \big ] \, , \end{aligned} \end{aligned}$$
(34)

where \(S_{GW}\) and \({\mathscr {A}}_{GW}\) characterizes the amplitude and polarization of the GWs. With the above definition, we can write,

$$\begin{aligned} \delta ({{\varvec{k}}}-{{\varvec{k}}}')\,S_{GW}= & {} \frac{1}{(2\pi )^3}{\mathscr {M}}_{ijlm}\langle h_{ij}'({{\varvec{k}}})\, h_{lm}^{*'}({{\varvec{k}}}')\rangle \, \end{aligned}$$
(35)
$$\begin{aligned} \delta ({{\varvec{k}}}-{{\varvec{k}}}')\,{\mathscr {H}}_{GW}= & {} \frac{1}{(2\pi )^3}{\mathscr {A}}_{ijlm}\langle h_{ij}'({{\varvec{k}}})\, h_{lm}^{*'}({{\varvec{k}}}')\rangle . \end{aligned}$$
(36)

We now choose a coordinate system, for which unit vectors are \(\hat{\mathbf{e}}_1\), \(\hat{\mathbf{e}}_2\), \(\hat{\mathbf{e}}_3\), in which GW propagates in the \(\hat{\mathbf{e}}_3\) direction. Further, we introduce

$$\begin{aligned} e_{ij}^\pm =-\sqrt{\frac{3}{2}} \,(\varepsilon _i^\pm \times \varepsilon _j^\pm ) \end{aligned}$$
(37)

which forms basis for a tensor perturbations and satisfy the following properties: \(\delta _{ij}\,e^\pm _{ij}\,=\,0\), \(\hat{k}_i\,e^\pm _{ij}\,=\,0\) and \(e^\pm _{ij}\,e^\mp _{ij}\,=3/2\) [53, 54]. The right-handed and left handed circularly polarized state of the GWs are represented by \(+\) and − sign respectively. In this basis, polarization tensor and modes of the GWs can be written as follows:

$$\begin{aligned} \varPi _{ij}({{\varvec{k}}})= & {} e_{ij}^+\varPi ^+({{\varvec{k}}})+e_{ij}^-\varPi ^-({{\varvec{k}}}) \, , \end{aligned}$$
(38)
$$\begin{aligned} h_{ij}({{\varvec{k}}}, \eta )= & {} h^+({{\varvec{k}}}, \eta )e^+_{ij}+h^-({{\varvec{k}}}, \eta )e^-_{ij} \, . \end{aligned}$$
(39)

On using equations in (38), Eq. (29) can be expressed as:

$$\begin{aligned} \delta ({\varvec{k}}-{\varvec{k}}')f(k)= & {} \frac{3}{2}\langle \varPi ^+({\varvec{k}})\varPi ^{+*}({\varvec{k}}') \nonumber \\&+\varPi ^-({\varvec{k}}) \varPi ^{-*}({\varvec{k}'}) \rangle \, \end{aligned}$$
(40)
$$\begin{aligned} \delta ({\varvec{k}}-{\varvec{k}}')\,g(k)= & {} -\frac{3}{2}\langle \varPi ^+({\varvec{k}})\varPi ^{+*}({\varvec{k}}') \nonumber \\&-\varPi ^-({\varvec{k}}) \varPi ^{-*}({\varvec{k}'}) \rangle \, . \end{aligned}$$
(41)

Adding and subtracting Eqs. (40) and (41) we obtain,

$$\begin{aligned} f(k) + g(k)= & {} 3\langle \varPi ^-({{\varvec{k}}})\varPi ^{-*}({{\varvec{k}}})\rangle \, \approx 3 \langle |\varPi ^-|^2\rangle \, \end{aligned}$$
(42)
$$\begin{aligned} f(k) - g(k)= & {} 3\langle \varPi ^+({{\varvec{k}}})\varPi ^{+*}({{\varvec{k}}})\rangle \approx 3 \langle |\varPi ^+|^2\rangle \, . \end{aligned}$$
(43)

Similarly, we can write Eq. (36) as:

$$\begin{aligned} \delta ({\varvec{k}}-{\varvec{k}}')S_{GW}(k,\eta )= & {} \frac{3}{2} \langle h^{+}({\varvec{k}},\eta )h^{+*}({\varvec{k}}',\eta )\\&+ h^{-}({\varvec{k}},\eta )h^{-*}({\varvec{k}}',\eta ) \rangle \\ \delta ({\varvec{k}}-{\varvec{k}}'){\mathscr {H}}_{GW}(k,\eta )= & {} -\frac{3}{2} \langle h^{+}({\varvec{k}},\eta )h^{+*}({\varvec{k}}',\eta )\\&- h^{-}({\varvec{k}},\eta )h^{-*}({\varvec{k}}',\eta ) \rangle \end{aligned}$$

Therefore, components of the GWs evolve as

$$\begin{aligned} h^{\pm ''}({{\varvec{k}}}, \eta )+ 2\,{a'\over a}\, h^{\pm '}({{\varvec{k}}}, \eta )+k^2 h^{\pm }({{\varvec{k}}}, \eta )=8\pi G \varPi ^{\pm }({{\varvec{k}}})\, , \end{aligned}$$
(44)

here \(\varPi _{ij}\) is the mean square root value of the transverse traceless part of the energy momentum tensor. In terms of dimensionless variable \(x=k\eta \), the above equation reduces to

$$\begin{aligned} h^{\pm ''} + 2 \frac{\alpha }{x}\, h^{\pm '} + h^{\pm } = \frac{s^\pm }{k^2}\, , \end{aligned}$$
(45)

where \(s^\pm (k,\eta ) = \left( {8\pi G\over a^2}\right) \sqrt{\frac{f(k)\,\mp \,g(k)}{3}}\) and the parameters \(\alpha =1\) and \(\alpha =2\) indicates the radiation dominated and the matter dominated epoch respectively. In the radiation dominated epoch, the homogeneous solution of the Eq. (45) are the spherical Bessel function \(j_0\) and \(y_0\). In our case, magnetic field is generated at \(\eta _{in}\) in the radiation dominated epoch due to chiral instability leading to anisotropic stress which in turn generates the gravitational waves. Thus, the general solution of Eq. (45) can be given as,

$$\begin{aligned} h^\pm (x) = c_1^\pm (x)\,j_0(x)\,+\,c_2^\pm (x)\,y_0(x) \end{aligned}$$
(46)

where \(c_1(x)\) and \(c_2(x)\) are undetermined coefficients which is given as

$$\begin{aligned} c_1^\pm (x)= & {} - \int _{x_{in}}^x dx' \frac{s^\pm (x')}{\mathrm{w}(x')\, k^2}\,y_0(x') \, \end{aligned}$$
(47)
$$\begin{aligned} c_2^\pm (x)= & {} \int _{x_{in}}^x dx' \frac{s^\pm (x')}{\mathrm{w}(x')\, k^2}\, j_0(x') \end{aligned}$$
(48)

where \(\mathrm{w}(x) = j_0(x) \, y_0'(x) - y_0(x) \, j_0'(x) = \frac{1}{x^2}\). We have calculated \(c_1^\pm (x)\) and \(c_2^\pm (x)\) using equations Eqs. (47) and (48) under the limits of \(x>1\). In this limit, the second term with \(y_0\) diverges, therefore, first term dominates over second one. In this case, in the radiation dominated epoch, the two polarizations of the tensor perturbations can be written as:

$$\begin{aligned} h^+(x)= & {} c_1^+(1)j_0(x) \nonumber \\= & {} -{90 \over \pi ^2 g_{\mathrm{eff}}} \sqrt{\frac{f-g}{3}}\,j_0(x)\,\mathrm{log}(x_{in}) \end{aligned}$$
(49)
$$\begin{aligned} h^-(x)= & {} c_1^-(1)j_0(x) \,\nonumber \\= & {} -{90 \over \pi ^2 g_{\mathrm{eff}}} \sqrt{\frac{f+g}{3}}\,j_0(x)\,\mathrm{log}(x_{in})\, , \end{aligned}$$
(50)

here \(x=1\) in \(c_1(1)\) signifies the value at the time of horizon crossing. After horizon crossing, these gravitational waves propagate without any hindrance. However, their energy and polarization stretched by the scale factor, similar to the magnetic radiation energy. The time derivative of the of the Eqs. (49) and (50) is

$$\begin{aligned} h^{+'}(x)= & {} -{90 \over \pi ^2 g_{\mathrm{eff}}} \sqrt{\frac{f-g}{3}}\,j_0'(x)\,\mathrm{log}(x_{in}), \end{aligned}$$
(51)
$$\begin{aligned} h^{-'}(x)= & {} -{90 \over \pi ^2 g_{\mathrm{eff}}} \sqrt{\frac{f+g}{3}}\,j_0'(x)\,\mathrm{log}(x_{in})\, , \end{aligned}$$
(52)

In real space, the energy density of the gravitational waves is defined as

$$\begin{aligned} \rho _{GW} = {1\over 16\,\pi \, G\, a^2}\, \langle h'_{ij}h'_{ij}\rangle . \end{aligned}$$
(53)

Note that a factor of \(a^2\) in the denominator comes from the fact that \(h'\) is the derivative with respect to conformal time. In Fourier space, the energy density of the gravitational wave is given as

$$\begin{aligned} S_{GW}(k) = \int \frac{dk}{k}\frac{d\,S_{GW}}{d\,\mathrm{log}\, k} \, \end{aligned}$$
(54)

with

$$\begin{aligned} \frac{d\,S_{GW}(k)}{d\, \text {log}k}= & {} \frac{k^3 }{2\,M_*^4\,a^2\,(2\pi )^6 G}(|h^{+'}|^2+|h^{-'}|^2) \, . \end{aligned}$$
(55)

With this definition, we can define power spectrum evaluated at the time of generation as

$$\begin{aligned} \frac{d\varOmega _{GW,s}}{d\mathrm{log}k}= & {} \frac{1}{(\rho _{c,s}/M_*^4)}\frac{dS_{GW,s}}{d\mathrm{log}k}\nonumber \\= & {} {16\pi k^3\over 3(2\pi )^6 a_s^2} \left( {90\over \pi ^2g_\mathrm{eff}}\right) ^2{f(k)\over H_s^2}\left[ j_0'(x)\mathrm{log}(x_{in})\right] ^2 . \nonumber \\ \end{aligned}$$
(56)

where \(\rho _{c,s}\) is the critical density of the universe at the time of generation of GW. Once gravitational waves are produced, they are decoupled from the rest of the Universe. This implies that the energy density of the gravitational waves will fall as \(a^{-4}\) and frequency redshifts as \(a^{-1}\). Hence, the power spectrum at today’s epoch can be given as

$$\begin{aligned} \frac{d\varOmega _{{GW,0}}}{d\mathrm{log}k}\equiv \frac{d\varOmega _{GW,s}}{d\mathrm{log}k}\,\left( {a_s\over a_0}\right) ^4 \left( {\rho _{c,s}\over \rho _{c,0}}\right) \, . \end{aligned}$$
(57)

Assuming that the Universe has expanded adiabatically which implies that the entropy per comoving volume is conserved leads to

$$\begin{aligned} {a_s\over a_0} = \left( {g_{\mathrm{eff,0}}\over g_\mathrm{eff,s}}\right) ^{1/3}\left( {T_0\over T_s}\right) \, , \end{aligned}$$
(58)

where we have used \(g_{\mathrm{eff}}\) for the effective degrees of freedom that contribute to the entropy density also. This is due to the fact that effective degrees of freedom that contribute to the energy and entropy densities are same at very high temperature. Therefore, Eq. (57), using Eq.  (58) reads as

$$\begin{aligned} \frac{d\varOmega _{{GW,0}}}{d\mathrm{log}k}= & {} \left( {g_{\mathrm{eff,0}}\over g_{\mathrm{eff,s}}}\right) ^{4/3}\left( {T_0\over T_s}\right) ^4 \left( {H_s\over H_0}\right) ^2\,\frac{d\varOmega _{{GW,s}}}{d\mathrm{log}k}\, \nonumber \\= & {} {16\pi k^3\over 3(2\pi )^6 a_s^2}\times \left( {90\over \pi ^2g_\mathrm{eff}}\right) ^2\left( {g_{\mathrm{eff,0}}\over g_{\mathrm{eff,s}}}\right) ^{4/3} \nonumber \\&\left( {T_0\over T_s}\right) ^4 {f(k)\over H_0^2}\left[ j_0'(x)\mathrm{log}(x_{in})\right] ^2 \end{aligned}$$
(59)

In Figs. 1 and 2, we have shown the variation of GW spectrum with respect to k for different temperature at fix number density and for different number density at fix temperature.

Fig. 1
figure 1

Gravitational wave spectrum as a function of \(x=k\eta \). We have fixed \(n = 10^{-6}\) and varied temperature

Full size image
Fig. 2
figure 2

Gravitational wave spectrum as a function of \(x=k\eta \) at different number densities at temperature T\(\sim 10^8\) GeV. For a large number density, the effects at large length scale is more than at the small length scale

Full size image

4 Results and Discussion

Before discussing the results obtained, we would like to explain various important length scales useful for the magnetohydrodynamic discussion of the generation of GWs due to the chiral instability in presence of the external magnetic fields. Firstly, magnetic modes grow exponentially for \(k= k_{\mathrm{ins}}= \xi ^{(B)}/2 \approx g^2 (\frac{\mu }{T})T\) [25, 26]. Secondly, dissipation due to the finite resistivity of the plasma works at wave numbers \(k< k_{\mathrm{diss}}\sim \sqrt{\frac{g_{\mathrm{eff}}}{6}}gT\) [55]. Therefore, near instability, the wavenumber corresponding to the length scales of interest \(k \ll k_{\mathrm{diss}}\). In the present analysis, we have not considered any dissipation in the plasma and restricted ourself to the scales where there is maximally growing modes of the magnetic fields are available. Therefore, for k values larger than the instability scale, this analysis may not be reliable. In Figs. 1 and 2, we have found that GWs peak occurs at \(k_{\mathrm{ins}}\). It is evident that at higher values of k i.e. at small length scales, power increases after instability. This is related to the rise in the magnetic energy at large k. The rise in the magnetic energy is unphysical as we know that in turbulent system, energy accumulates at smallest scales. This effect in principle can be restricted by going to hyper diffusion scale (instead of \(\nabla ^2\) operator, one needs to introduce \(\nabla ^4\) operator) [56].

We would like to emphasize, in Figs. 1 and 2 that variable \(x=k\eta \) is a dimensionless quantity. In order to interpret the results, we convert x in frequency of the GW. Since, \(x=a\, k\eta = \frac{k}{T}\eta ~= \frac{2\pi \nu }{T}\eta \). From this we can get \(\nu \) in terms of x as \(\nu = \frac{T}{2\pi \eta } x\). Moreover, peak of the power spectrum of the GW occurs when growth of the instability is maximum which is given by \(\nu _{\mathrm{max}}\approx \frac{16}{9\pi }\left( \frac{\delta ^2}{\sigma /T}\right) T\) [25, 55]. Here \(\delta \) is defined as \(\delta = \alpha ( \mu _R-\mu _L) /T\). The red-shifted value of the frequency can be obtained using the relation: \(\nu _0=\frac{T_0}{T_*}\nu _{\mathrm{max}}\), where \(T_*\) is the temperature at which instability occurs. Hence, we can obtain the frequency at which maximum power is transferred from the magnetic field to GWs. The obtained formula of the frequency in simplified form is \(\nu \approx 10^{9} \, \delta ^2\) Hz, where we have used \(T_0=2.73\) K \(\approx 10^{11}\) Hz in our units and \(\sigma /T=100\). For temperature \(T\sim 10^6\) GeV, \(( \mu _R-\mu _L)/T \sim 10^{-3}\) [23] and thus, \(\delta ^2= 10^{-10}\) (with \(\alpha \approx 10^{-2}\)). Hence frequency where maximum power of GW occurs, is around \(10^{-1}\) Hz. Thus they may be detected in eLISA experiment [57]. Further, the strength of magnetic field changes when chiral charge density n change. Figure  2 shows the effect of n on GW spectrum. It is apparent that the \(k_{\mathrm{ins}}\) is not affected by the number density and hence, the peak does not shift. However, the power in a particular k mode enhances with an increase in n. This happens due to the fact that for a larger value of n, magnetic field strength is higher at larger k [26].

5 Conclusion

In the present work, we have extended our earlier works on the generation of primordial magnetic fields in a chiral plasma [25, 26] to the generation of GWs. This kind of source may exist much above electroweak scale. We have shown that the gravitational anomaly generates the seed magnetic field which evolves and create instability in the system. This instability acts as a source of anisotropic stress which leads to the production of gravitational waves. The production and evolution of the magnetic field has been studied using Eq. (11). In order to obtain the velocity profile, we have used scaling properties [48, 49] rather than solving the Navier-Stokes equation. This scaling property results in more power in the magnetic field at smaller k as compared to that of the case without scaling symmetries (see [26]). We have calculated power spectrum of the produced GWs and shown that the spectrum has a distinct peak at \(k_{\mathrm{ins}}\) and hence correspond to the dominant frequency of GW. The GW generated at high temperature \(T \ge 10^{6}\) GeV via aforementioned method is potentially detectable in eLISA.

In this work, we have considered massless electrons much above electroweak scale and discussed the production of gravitational waves due to chiral instabilities in presence of Abelian fields belonging to \(U(1)_Y\) group. However, a similar situation can arise in the case of Quark-Gluon Plasma (QGP) at \(T\gtrsim 100\) MeV where quarks are not confined and interact with gluons which may result in instabilities. Thus, GW can be produced in QGP as well.

To conclude, the study of relic GWs can open the door to explore energy scales beyond our current accessibility and give insight into exotic physics.