---
abstract: 'Given that gravity waves are the future probes of the universe, it is important to test various physical effects of these on matter. This article studies a perturbation of the field of a massless scalar particle caused by a gravity wave. We also discuss the gravitational waves’ effect on a chiral fermion. Our results have physical implications for early universe cosmology and also for ground based observations.'
author:
- 'Siddhartha E. Morales G.$^{1, }$ $^2$'
- Arundhati Dasgupta$^1$
title: Scalar and Fermion field interactions with a gravity wave
---

i

ł Ł c

Introduction
============

With the discovery of gravity waves at the LIGO detector in 2016, a new era of observational physics has started [@ligo; @ligo2]. The gravity wave will soon become a major probe of the universe and it’s dynamics. In this article, we discuss the question: how does a gravity wave affect the behaviour of a scalar wave? This scalar wave could be a matter wave in a condensed matter system, or could represent the inflaton of the early universe [@cosmo]. Gravity waves are important in the dynamics of the early universe, though they are expected to be non-linear and at a much higher amplitude than that detected on earth recently [@ligo]. In our calculations, we test whether a linear gravity wave can interact with a scalar field to give non-trivial dynamics. The relevance for early universe cosmology has to be explored further. Whereas we have done a search for papers which report a similar scalar field gravity wave interaction [@singleton], our calculations are unique. In particular, we seem to have obtained a ‘non-trivial’ interaction of the massless scalar wave with the gravity wave, which we discuss and try to obtain new physical implications. The scalar field responds to the passing gravitational wave, and is perturbed to an oscillatory reaction evident energetically, The simple superposition of two massless waves in the same direction usually gives us a resultant mode of frequency $\omega\pm\omega_0$. If we take a gravitational wave in the $z$ direction, then in the perturbation equations, the inhomogeneous terms are proportional to $k_{0x}^2 - k_{0y}^2$ and $k_{0x} k_{0y}$ (where $k_{0x}, k_{0y}$ are the $x-y$ components of the scalar field’s direction of propagation). Thus, at least one of the transverse directions of the scalar wave vector has to be non-zero for the gravity wave in the $+z$ direction to have any non-trivial effect on the scalar wave. The spherical propagator used to obtain the scalar wave perturbation is also crucial for the results, if the scalar waves are obtained using the planar wave kernel as we show in the massive scalar example, the plane waves are perturbed as ordinary waves. Our solutions, are obtained using Kirchoff’s law and Duhamel’s principle, and the separation of variables method is not used to solve the Laplacian. The new result is that the perturbations have a new mode, obtained using the spherical propagator for the solution to the inhomogeneous wave equation. The new mode can be interpreted as massless waves of frequency $\omega$ with the dispersion $\omega^2= (\textbf {k}+ \textbf {k}_0)^2$, where $\textbf k$ is the wave vector of the gravity wave, and $\textbf k_0$ is the wave vector of the scalar/fermion wave. We check for resonant modes, in the solutions for the scalar field as has been found in ‘driven systems’ and find none.

Similarly we investigate the effect of the gravitational wave on a chiral fermion (neutrino), and obtain the first order effect on the particle. Similar to the scalar particle, our solutions have different behaviour not seen previously in [@fermion]. We comment on this, and how this might affect gravitational interferometers on earth. As the neutrino density is modified, we can expect corrections to the effective number of neutrino species for the cosmic neutrino background. The effective number of neutrino species is around the physically expected 3.0 $\pm$ corrections, the fluctuations have been attributed to non-thermal interactions in the early universe [@neut1; @neut2]. Here we estimate a correction due to the gravity wave interaction, and the perturbation effect’s contribution to the background cosmic neutrino density. The perturbation of the density is very small to contribute to the effective degrees of freedom, however the neutrino interaction with the cross polarization has a resonant mode. We comment on how the resonant run away mode can contribute to the background density of neutrinos.

These perturbations might also affect the neutrino flux for the multi-messenger astronomy observations from neutron star collisions. As of now neutrinos are yet to be detected from the gravity wave events. More importantly, the fluctuations in the neutrino flux produced due to the gravity wave interaction can be used to detect gravity waves on earth. Particularly for fermions perturbed by the cross polarization of the gravity waves, there is a resonant mode which grows with time. The physical implications of this have to be investigated in a real experimental situation in consultation with existing neutrino experiments [@adg].

The paper is organized as follows: in the next section, we describe the Klein-Gordon equation and its solution in the gravity wave background. In the third section we describe the chiral fermion and its solution in the gravity wave background. In the fourth section, we interpret the solution, and provide some practical uses of the solution, and we conclude.

\[sec: Klein Gordon\] The Klein Gordon equation on the background of a gravity wave\
====================================================================================

The Klein-Gordon equation in the background of a curved metric $g_{\mu \nu}$ is $$\partial_\mu \left( \sqrt{\text{-}g} \  
 \partial^\mu \phi \right) =  \partial_{\nu} \left( \sqrt{\text{-}g} \ g^{\mu\nu} \partial_{\mu} \phi \right)= 0.$$

In this article we solve this equation in the background where $g_{\mu \nu}$ is due to a gravity wave. We take the example of a perturbation due to gravity wave of a flat geometry, which could be the situation on Earth, when the gravity wave arrives from a distant source.

\[sec: Gravity Wave\] Gravity Waves\
------------------------------------

The simplest kind of gravity waves emerge in the frame of the linearized gravity theory, obtained by approximating the metric tensor up to first order. In the linear theory approximation one can say that the metric is flat space-time with fluctuations $h_{\mu \nu}$ as $$g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu},$$ where $\lVert h_{\mu\nu} \rVert << 1$. Using this, it can be proved that the inverse of the metric can be written as $$g^{\mu\nu} = n^{\mu\nu} - h^{\mu\nu}$$ (if we want to keep everything up to first order).\
Introducing the “traceless" metric $$\bar{h}_{\mu\nu} = h_{\mu\nu} - \frac{1}{2}\eta_{\mu\nu}h^{\alpha}_\alpha$$ Einstein’s equation can be written as $${\mathop{}\!\mathbin\Box}\bar{h}_{\mu\nu} = -16 \pi T_{\mu\nu}.$$

(in units of $G=c=1$ or Newton’s constant and speed of light set to 1). In the absence of matter the last equation is just the three-dimensional wave equation, which has as solutions plane waves. $$\bar{h}_{\mu\nu} = A_{\mu\nu} \ \mathrm{e}^{i k_{\alpha}x^{\alpha}}$$ where $A_{\mu \nu}$ is the tensor polarization amplitudes. Moreover, using gauge freedom we can choose the so called Transverse-Traceless gauge (TT-gauge), which impose on $\bar{h}_{\mu\nu}$ the following conditions, $\bar{h}_{\mu}^\mu = 0$, then $\bar{h}_{\mu\nu} = h_{\mu\nu}$ and $\partial_\alpha h_{\mu\nu} = 0$. The metric is divergence-less $ h_{\mu\nu} u^\mu = 0$, for some 4-velocity $u$.

Moreover, since $h_{\mu\nu}$ is symmetric it has only $10- 4 -4 = 2$ independent components. We can set the coordinate axes so that the gravity wave is propagating in the $z$ direction. With these conditions and the TT-gauge, the linearized metric tensor can be expressed as $$\begin{aligned}
g_{\mu\nu} = \begin{pmatrix}
    -1 & 0 & 0 & 0 \\
    0 & 1+h_+ & h_\times & 0 \\
    0 & h_\times & 1-h_+ & 0 \\    
    0 & 0 & 0 & 1
  \end{pmatrix}    
 \,\,\,\,\,\, \ ; \ \ \ \ \ \ \ \ 
 g^{\mu\nu} = \begin{pmatrix}
    -1 & 0 & 0 & 0 \\
    0 & 1-h_+ & -h_\times & 0 \\
    0 & -h_\times & 1+h_+ & 0 \\    
    0 & 0 & 0 & 1
  \end{pmatrix}  \end{aligned}$$

(where $h_+ = A_+ \cos(\omega(z-t))$ and $h_{\times}=A_\times \cos(\omega(z-t)+ \delta)$, $\omega$ the frequency of the wave and $A_{+}, A_{\times}$ being the amplitudes of the two polarizations)

Let $g$ be the determinant of $g_{\mu\nu}$, we can see that due to the block diagonal form of the metric, the determinant is $$g = -(1 - h_{+}^2 - h_{\times}^2) \approx -1$$ which is approximately $1$ up to first order.

For numerical estimates; the binary white dwarf systems in our galaxy, a good approximation of the wave strain is [@flan]

$$h \approx 10^{-21} \left( \frac{M}{2M_\odot} \right)^{5/3} \left( \frac{1 \, \mbox{hour}}{P} \right)^{2/3} \left( \frac{1 \, \mbox{Kpsc}}{r} \right).$$

For binary neutron stars spiralling together in our galaxy, a good approximation of the wave strain is

$$h \approx 10^{-22} \left( \frac{M}{2.8 M_\odot} \right)^{5/3} \left( \frac{0.01 \, \mbox{sec}}{P} \right)^{2/3} \left( \frac{100 \, \mbox{Mpsc}}{r} \right)$$

where $P=2\pi/\Omega$, $\Omega$ is the orbital frequency and $r$ is the distance from the source. M is the total mass and $M_{\odot}$ is the mass of the sun. Similarly, quasi-normal ring down of stellar mass black hole mergers are expected to give waves of amplitudes [@ligo2] $$h \approx 10^{-21} \left(\frac{\nu}{0.25}\right) \left(\frac{M}{20 M_{\odot}}\right)\left(\frac{200 Mpsc}{r}\right)$$ ($\nu$ is the symmetric mass ratio [@ligo2]).

\[sec: Wave equation in 3d\] The interaction of the scalar wave and the gravity wave\
=====================================================================================

Let’s say that we have a field satisfying Klein-Gordon’s equation for a massless particle in a metric produced by a gravity wave. The equation for the scalar field in a curved space is given by $\partial_\mu \left( \sqrt{\text{-}g} \partial^\mu \phi \right) =  \partial_{\nu} \left( \sqrt{\text{-}g} g^{\mu\nu} \partial_{\mu} \phi \right)= 0 $, and remembering that the determinant for the metric of a gravity wave is equal to $-1$, the equation is just $$\partial_{\nu} \left(  g^{\mu\nu} \partial_{\mu} \phi \right)= 0.$$ Expanding all components in the equation and noting that $g_{\mu\nu}$ has only $6$ six non-zero components and that $h_{\mu\nu}$ is function only of $z$ and $t$, the equation is explicitly $$0 = -\partial_{tt} \phi +\partial_{xx} \phi +\partial_{yy} \phi +\partial_{zz} \phi - \left[ h_+ ( \partial_{xx} \phi - \partial_{yy} \phi ) + 2h_{\times} \partial_y \partial_x \phi \right]$$ or in a more compact form $${\mathop{}\!\mathbin\Box}\phi - \left[ h_+ ( \partial_{xx}  - \partial_{yy} ) + 2h_{\times} \partial_y \partial_x  \right] \phi=0$$

Due to the smallness of the strain of a gravity wave $(10^{-21}-10^{-22})$, we can make a Taylor expansion up to first order in $A_+$ and $A_\times$ as continuous variables, thus this result would be $$\phi \approx \phi_0 + A_+ \phi_+ + A_\times \phi_\times .$$

Substituting this expression into the Klein-Gordon with the metric of a gravity wave, and keeping terms up to first order gives $$\begin{aligned}
 {\mathop{}\!\mathbin\Box}\phi_0 - \left[ A_+ \cos(\omega(z-t)) (\partial_{xx} \phi_0 - \partial_{yy} \phi_0) + 2A_{\times}\cos(\omega(z-t) + \delta) \partial_x \partial_y \phi_0 \right] && \nonumber \\+ A_+ {\mathop{}\!\mathbin\Box}\phi_+ + A_\times{\mathop{}\!\mathbin\Box}\phi_\times&=&0.\end{aligned}$$

Since $A_+$ and $A_\times$ are arbitrary constants we get the following set of equations $$\begin{aligned}
{\mathop{}\!\mathbin\Box}\phi_0 &= 0, \ \ \ \mbox{(unperturbed wave equation) } \\
{\mathop{}\!\mathbin\Box}\phi_+ &= \cos(\omega(z-t)) (\partial_{xx} - \partial_{yy})\phi_0 \\
 {\mathop{}\!\mathbin\Box}\phi_{\times} &= 2\cos(\omega(z-t)+\delta) (\partial_{x} \partial_{y})\phi_0\end{aligned}$$

For sake of simplicity let’s suppose that $$\phi_0 (\textbf{x}, t)= \Re \left( A_0\mathrm{e}^{i(\textbf{k}_0 \cdot \textbf{x} - \omega_0 t)} \right)$$

(Notation: we use boldface to represent three vectors)

The perturbations are next solved for using techniques of inhomogeneous PDE’s.

The wave equation in 3+1 d
--------------------------

The wave equation in $\mathbb{R}^3$ is just the D’ Alambertian of a function equated to zero, ${\mathop{}\!\mathbin\Box}\phi = 0$. Further, let’s suppose that we have Cauchy (initial) conditions, $$\begin{aligned}
{\mathop{}\!\mathbin\Box}\phi &= 0 \\
\phi(\textbf{x},0) &= F(\textbf{x}) \\
\frac{\partial \phi}{\partial t} (\textbf{x},0) &= G(\textbf{x})\end{aligned}$$ where $\textbf{x} \ \in \mathbb{R}^3$, $t>0$, $ F  \in C^3(\mathbb{R}^3)$ and $G \ \in \ C^2(\mathbb{R}^3)$.

The solution for the wave equation with Cauchy conditions is given by the theorem [@kirch].

Suppose that $ F  \in C^3(\mathbb{R}^3)$ and $G \ \in \ C^2(\mathbb{R}^3)$. Then the solution of the initial value problem of the wave equation in three dimensions is given by $$\phi(\textbf{x},t) = \frac{1}{4\pi t} \int_{S(\textbf{x},t)} G(\textbf{x}') d \sigma_t + \frac{\partial}{\partial t} \left[ \frac{1}{4 \pi t} \int_{S(\textbf{x},t)} F(\textbf{x}') d \sigma_t \right]$$ and the solution is in $C^2$ for $\textbf{x} \in \mathbb{R}^3$ and $t \leq 0$.

The above formula is called Kirchoff’s formula.

The solution for the inhomogeneous wave equation is given by *Duhamel’s principle*, which says that if $\phi(\textbf{x},t)$ is the solution of the following initial value inhomogeneous wave equation, [@pde] $$\begin{aligned}
{\mathop{}\!\mathbin\Box}\phi = h(\textbf{x},t); \ \ \ \textbf{x} \in \ \mathbb{R}^3, \ \ \ t > 0, \\
\phi(\textbf{x},0) = 0, \ \ \ \ \frac{\partial \phi}{\partial t}(\textbf{x},0) = 0; \ \ \textbf{x} \in \ \mathbb{R}^3.\end{aligned}$$ Let $v(x,t;\tau)$ be the solution of the associated “pulse problem" $$\begin{aligned}
{\mathop{}\!\mathbin\Box}v = 0; \ \ \ \textbf{x} \in \mathbb{R}^3, \ \ \ t > \tau, \\
v(x,\tau;\tau) = 0, \ \ \ \frac{\partial}{\partial t} v(x,\tau:\tau) = - h(\textbf{x},\tau)\end{aligned}$$ Then the solution for the inhomogeneous wave equation is given by $$\phi(\textbf{x},t) = \int_{0}^t v(x,t:\tau) \ d \tau$$

Combining Duhamel’s principle and Kirchoff’s formula, the solution for the inhomogeneous wave equation in three dimensions can be written in compact form as $$\phi(\textbf{x},t) = -\frac{1}{4\pi} \int_{\bar{B}(\textbf{x},t)} \frac{h(\textbf{x}',t-r)}{r}d^3 x'$$ where $r = \lVert \textbf{x} - \textbf{x}' \rVert$ and $\bar{B}(\textbf{x},t)$ is the ball in $\mathbb{R}^3$ with center at $\textbf{x}$ and radius $t$.

The solution of $\phi_+$ and $\phi_{\times}$ would be $$\begin{aligned}
\phi_+ (\textbf{x},t) &= \frac{1}{4\pi} (k_{0x}^2 - k_{0y}^2) \int_{\bar{B}(\textbf{x},t)} d^3 \textbf{x}' \frac{\phi_0 (\textbf{x}',t')}{r} \cos(\omega(z'-t')) \\
\phi_\times (\textbf{x},t) &= \frac{1}{4\pi} (2k_{0x}k_{0y}) \int_{\bar{B}(\textbf{x},t)} d^3 \textbf{x}' \frac{\phi_0 (\textbf{x}',t')}{r} \cos(\omega(z'-t')+\delta)\end{aligned}$$ where $\textbf{k}_0 = (k_{0x},k_{0y},k_{0z}) $, $r = \lVert \textbf{x}' - \textbf{x} \rVert$ and $t' = t - r$.

Now, in order to find the perturbation explicitly, let’s solve the more general integral. $$I(\delta_0,\delta) = \int_{\bar{B}(\textbf{x},t)} d^3 \textbf{x}' \ \mathrm{e}^{i(\textbf{k}_0 \cdot \textbf{x}' - \omega_0 (t') + \delta_0 )}  \frac{\cos(\omega(z'-t') + \delta)}{r}$$ where $r = \lVert \textbf{x}' - \textbf{x} \rVert$ and $t'=t-r$. Since we are integrating over the ball centred at $\textbf{x}$ we can make a change of variable (reference system) and center the ball at the origin, the translation would be $$\begin{aligned}
\textbf{x}' &= \textbf{x} + \textbf{r}  \label{eqn:define}\\
d^3 \textbf{x}' &= r^2 \sin(\theta) d\phi d\theta dr \nonumber\end{aligned}$$

Moreover, let’s write the cosine term in its exponential form, then the integral is converted into $$\begin{aligned}
I &= \int_{\bar{B}(0,t)} \mathrm{e}^{i(\textbf{k}_0 \cdot (\textbf{x} + \textbf{r}) - \omega_0(t-r)+\delta_0)} \frac{\mathrm{e}^{i(\omega(z' - (t-r)) + \delta)}+\mathrm{e}^{-i(\omega(z' - (t-r)) + \delta)}}{2r} r^2\sin\theta d\phi d\theta dr \end{aligned}$$

This integral can be done using [@gradshteyn] and Mathematica, the result is: $$\mathrm{e}^{i(\textbf{k}_0 \cdot (\textbf{x}) - \omega_0(t)+\delta_0)} \times$$ $$\frac{\pi \mathrm{e}^{i(\omega(z - t) + \delta)}}{\lVert \textbf{k} + \textbf{k}_0\rVert} \frac{ \left(-\lVert \textbf{k} + \textbf{k}_0\rVert+e^{i t \left(\omega +\omega _0\right)} \left(\lVert \textbf{k} + \textbf{k}_0\rVert \cos (\lVert \textbf{k} + \textbf{k}_0\rVert t)-i \left(\omega +\omega _0\right) \sin
   (\lVert \textbf{k} + \textbf{k}_0\rVert t)\right)\right)}{\omega \omega_0 (1-\cos\beta)} + (\omega\rightarrow -\omega)$$ In the above $\beta$ is the angle between the wave vectors of the scalar wave and the gravity wave in three dimensions. From the above one can see that apart from the usual super-position of two waves where the frequencies add or subtract, there is a term of the form $\cos(|| \textbf{k} + \textbf{k}_0|| \ t)$. To analyze the results, we write the waveform in the following way: $$\frac{\pi}{ \omega\omega_0(1-\cos\beta)} e^{\left(i \left(\textbf{k}_0+\textbf{k}\right)\cdot\textbf{x}\right)} \left[-e^{-i(\omega+\omega_0) t} +   \left[ \cos(||\textbf{k}+\textbf{k}_0||  ~ t) - i \frac{(\omega+\omega_0)}{||\textbf{k}+\textbf{k}||}\sin(||\textbf{k}+\textbf{k}_0||~t)\right]\right] + (\omega\rightarrow -\omega)
 \label{eqn:masslesssc}$$

A plot of the above using MAPLE shows the behaviour of the solution in time. The first part comprises of an ordinary part of frequencies $\omega\pm\omega_0$, and the second part comprises of a wave of frequency $||\textbf{k}+\textbf{k}_0|| $, as in Figure (\[fig:Scalar\]).

The solution does not show resonance, as if we set $\beta=0$ (the angle in momentum space between the gravity wave and the scalar wave), which will make the amplitude singular, the solution itself is zero. Despite this, the solution shows the promise that matter waves, can be used to detect the passing of a gravity wave on earth. Note that the unperturbed wave does not show any amplitude oscillation, and is simply propagating with constant energy. The presence of the gravity wave causes interesting oscillation patterns as in Figure(\[fig:Scalar\]) albeit these perturbations will be very weak. The physical reason for the oscillating amplitude is of course the ‘energetic interaction’ of the scalar wave with the gravity wave in response to the stretching and straining of the space-time through which it propagates. To estimate the amplitude of the perturbation let’s say that the gravity wave has only $+$ polarization and that $\textbf{k}_0$ propagates only in the $x$ direction, then $$\begin{aligned}
\Delta \phi & \approx A_{+} \phi_{+} \sim h \phi_+ \\
& \sim  h \frac{A_0 k_{0}^2}{2\omega\omega_0} \sim hA_0 \left( \frac{\omega_0}{\omega} \right)\end{aligned}$$

Then, the order of the perturbation is proportional to $ h \left( \frac{\omega_0}{\omega} \right)$, where $h$ is the gravity wave strain. We already know that this strain is of the order of $10^{-21} - 10^{-22}$, and for $ \omega_0 = 160.23 Ghz$, the peak CMB frequency as an example; [^1] and $\omega \in (10^{-7} Hz -100 Hz)$, the amplitude of the perturbation would be in the range of $10^{-13} - 10^{-3}$. Even if it is weak compared to the scalar wave, which we take as order 1, it is certainly detectable, the frequency of oscillation would be however of the order of $\sqrt{\omega_0^2+\omega^2}\approx \omega_0$. If we take the ‘oscillatory’ phase of the inflaton, and use its mass to set a frequency scale then that is related to the Planck mass and $\omega_0\sim 10^{-6} M_{\rm pl}$.

\[sec: KG massive\] Klein-Gordon for a massive particle\
========================================================

Now, let’s consider the Klein-Gordon equation for a massive particle

$$-\frac{1}{\sqrt{\text{-}g}}\partial_\mu \left( \sqrt{\text{-}g} \partial^\mu \phi \right)+m^2 \phi =-  \partial_{\nu} \left( g^{\mu\nu} \partial_{\mu} \phi \right) +m^2 \phi = 0 ,$$

along with the metric of linear gravity, so $g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu}$ and $g \approx -1$. Expanding all the terms in the equation and using the same calculations as for the massless particle we get

$$0 = {\mathop{}\!\mathbin\Box}\phi -m^2\phi - \left[ h_+ ( \partial_{xx}  - \partial_{yy} ) + 2h_x \partial_y \partial_x  \right] \phi$$

Taking the first order approximation (with respect to the gravity wave strain), we can take $\phi \approx \phi_0 + A_+ \phi_+ + A_\times \phi_\times $. Substituting this approximation into the Klein-Gordon equation would give $$\begin{aligned}
{\mathop{}\!\mathbin\Box}\phi_0 -m^2(\phi_0 + A_+ \phi_+ + A_\times \phi_\times) - \left[ A_+ \cos(\omega(z-t)) (\partial_{xx} \phi_0 - \partial_{yy} \phi_0)+ 2 A_\times\cos(\omega(z-t) + \delta) \partial_x \partial_y \phi_0 \right] &&  \\ + A_+ {\mathop{}\!\mathbin\Box}\phi_+ + A_\times{\mathop{}\!\mathbin\Box}\phi_\times&=&0 \nonumber\end{aligned}$$

Since $A_+$ and $A_\times$ are arbitrary constants we get the following set of equations $$\begin{aligned}
{\mathop{}\!\mathbin\Box}\phi_0  - m^2 \phi_0 &= 0, \ \ \ \mbox{(unperturbed Klein-Gordon equation) } \\
{\mathop{}\!\mathbin\Box}\phi_+ - m^2 \phi_+ &= \cos(\omega(z-t)) (\partial_{xx} - \partial_{yy})\phi_0 \\
 {\mathop{}\!\mathbin\Box}\phi_\times  - m^2 \phi_\times &= 2  \cos(\omega(z-t)+\delta) (\partial_{x} \partial_{y})\phi_0\end{aligned}$$

The simplest solution for the unperturbed Klein-Gordon equation is just a plane wave with a modified dispersion relation: $ k_{\mu}k^\mu=-\omega^2 + \lVert \textbf{k}\rVert^2 =- m^2$. Let’s take for simplicity $\phi_0$ as a plane wave with one frequency, $\omega_0$. Then the solution for the first equation would be $$\phi_0 = A_0 \mathrm{e}^{i(\textbf{k}_0 \cdot \textbf{x} - {\omega}_0 t)}$$

In order to solve the inhomogeneous Klein-Gordon equation, we have to convert it into a momentum space through a Fourier transform. Then, if we have the equation $ {\mathop{}\!\mathbin\Box}\phi - m^2 \phi = f$ the solution is given through the propagator

$$G(\textbf{x}) = \frac{1}{\sqrt{(2\pi)^2}} \int_{\mathbb{R}^4} \frac{\mathrm{e}^{ik_\mu x^\mu}}{k_\mu k^\mu + m^2}$$

and the solution for $\phi$ would be $$\phi(x) = - \lim_{\epsilon \rightarrow 0} \int_{\mathbb{R}^4} \frac{d^4 k}{(2\pi)^4} \mathrm{e}^{i(k_\mu x^\mu)} \frac{\tilde{f}(\textbf{k})}{-k_\mu k^\mu - m^2 + i \epsilon}$$ Where $\tilde{f}$ is the Minkowski-Fourier transform of $f$. which is give by $$\tilde{f}(\textbf{k}) = \int_{\mathbb{R}^4} d^4 x \  \mathrm{e}^{-i(k_\mu x^\mu)} f(\textbf{x})$$

Now let’s calculate the integral for $f(\textbf{x}) =  \cos(\omega(z-t)) (\partial_{xx} - \partial_{yy})\phi_0 = A_0(k_{0y}^2 - k_{0x}^2)\cos(\omega(z-t)) \mathrm{e}^{ik_{0\mu}x^\mu}$. Then $$\begin{aligned}
\tilde{f}(\textbf{k}) &= A_0 (k_{0y}^2 - k_{0x}^2) \int_{\mathbb{R}^4} d^4 x \  \mathrm{e}^{-i(k_\mu x^\mu)} \cos(\omega(z-t)) \mathrm{e}^{i k_{0\mu}x^\mu} \\
&=  A_0  (k_{0y}^2 - k_{0x}^2) \int_{\mathbb{R}^4} d^4 x \ \mathrm{e}^{-i(k_\mu x^\mu)} \mathrm{e}^{ik_{0\mu}x^\mu} \left( \frac{\mathrm{e}^{i\omega(z-t)} +\mathrm{e}^{-i\omega(z-t)}}{2}\right) \\
&= A_0 (k_{0y}^2 - k_{0x}^2) \int_{\mathbb{R}^4} d^4 x \ \mathrm{e}^{-i(\textbf{k}\cdot \textbf{x} - \omega t)} \mathrm{e}^{i(\textbf{k}_0\cdot \textbf{x} - \omega_0 t)} 
\left( \frac{\mathrm{e}^{i\omega(z-t)} +\mathrm{e}^{-i\omega(z-t)}}{2}\right) \\
&=\frac{A_0}{2} (k_{0y}^2 - k_{0x}^2) (2\pi)^4 \left( \delta(\omega - (\omega_0 + \omega))\delta(k_{0x} - k_x)\delta(k_{0y} - k_y) \delta(k_{0z} + \omega- k_z) +   \right. \\
 &\left. \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \delta(\omega - (\omega_0 - \omega))\delta(k_{0x} - k_x)\delta(k_{0y} - k_y)\delta(k_{0z} - \omega - k_z)\right) \nonumber\end{aligned}$$ Thus the solution for the inhomogeneous Klein-Gordon equation would be $$\begin{aligned}
 &= -\frac{A_0}{2} (k_{0y}^2 - k_{0x}^2)  \lim_{\epsilon \rightarrow 0} \int_{\mathbb{R}^4} d^4 k \  \mathrm{e}^{i(k_\mu x^\mu)}  \cdot\\
&\cdot \frac{\left( \delta(\omega \mbox{-} (\omega_0 + \omega))\delta(k_{0x} \mbox{-} k_x)\delta(k_{0y} \mbox{-} k_y) \delta(k_{0z} + \omega \mbox{-} k_z) +  \delta(\omega \mbox{-} (\omega_0 \mbox{-} \omega))\delta(k_{0x} \mbox{-} k_x)\delta(k_{0y} \mbox{-} k_y)\delta(k_{0z} \mbox{-} \omega \mbox{-} k_z)\right)}{-k_\mu k^\mu \mbox{-} m^2 + i \epsilon} 
\nonumber \\
&=  \frac{A_0}{2} (k_{0x}^2 - k_{0y}^2)\left( \frac{e^{i((\textbf {k}_0 + \textbf{k}) \cdot \textbf {x}
- (\omega_0 + \omega) t)}}{(\omega_0 + \omega)^2 - \lVert\textbf{k}_0 + \textbf{k} \rVert^2 - m^2} + \frac{\mathrm{e}^{i((\textbf{k}_0 - \textbf{k}) \cdot \textbf{x} - (\omega_0 - \omega)t)}}{(\omega_0 - \omega_g)^2 - \lVert \textbf{k}_0 - \textbf{k} \rVert^2 - m^2} \right) \\
&=\frac{A_0}{2} (k_{0x}^2 - k_{0y}^2) \left( \frac{\mathrm{e}^{i((\textbf{k}_0 + \textbf{k}) \cdot \textbf{x} - (\omega_0 + \omega)t)}}{2 \omega (\omega_0- k_0 \cos \beta)} - \frac{\mathrm{e}^{i((\textbf{k}_0 - \textbf{k}) \cdot \textbf{x} - (\omega_0 - \omega)t)}}{2 \omega (\omega_0- k_0 \cos \beta)} \right) \end{aligned}$$

Where $\cos \beta = \frac{\textbf{k} \cdot \textbf{k}_0}{ \omega \lVert \textbf{k}_0\rVert}$ and $k_0 = \lVert \textbf{k}_0 \rVert$. Here the solution to the massless Klein Gordon equation clearly shows two resultant waves of the scalar response, one with frequency $\omega_0+\omega$, and another with frequency $\omega_0-\omega$, as expected from interaction of two different waves. Note that in the way that the perturbed scalar wave is calculated, only the z-sector of the wave gets a frequency modification due to the superposition of the scalar wave and the gravity wave in the z-direction. There are no interferences or amplitude oscillations as observed in the massless example. This is due to the fact that in the Fourier transformed space, the waves superpose as two plane waves, with no added mode components of the perturbed wave. If we use the more difficult calculation of obtaining the emergent wave using the spherical Green’s function in position space, we again recover the interference pattern as obtained in the massless example. The position space Green’s function in spherical coordinates is G\_m( [**r**]{}, [**r**’]{}, t-t’) =1[4]{} \[eqn:green\] The above propagator is taken from Wikipedia and verified explicitly. Using this again gives interference patterns which could be used to detect gravitational waves on earth. The perturbation would be obtained as \_[+]{}\~(k\_[0x]{}\^2 - k\_[0y]{}\^2) A\_0 12 (e\^[i (z’-t’)]{} + e\^[-i (z’-t’)]{}) e\^[-i \_0 t’ + i **k**\_0 **x**’]{} G\_m( [**r**]{}, [**r**’]{}, t-t’)  d\^4 x’ \[eqn:mass\] One the $t'$ integral is obtained the form acquires expression very similar to the massless case and one obtains (details given in Appendix 2). \_+&\~& A\_0 (k\_[0x]{}\^2-k\_[0y]{}\^2) This is very close to the massless perturbation solution with $\tilde\omega=\omega+\omega_0$. As evident, such oscillation patterns can be detected in an actual experiment, even for the massive scalar.

\[sec: Dirac + \] Field perturbations for massless spin $1/2$ particles\
========================================================================

In this section we investigate the effect of the gravity wave on a massless fermion. This is motivated from the fact that most particles are fermions, and in particular neutrinos are almost massless fermions. We investigate the effect of a gravity wave on a ‘neutrino’, assumed massless in the first approximation.

Lets consider again a gravity wave in the TT gauge, but with only one of its two polarizations, $+$ or $\times$. The Dirac equation is written as $$i\gamma ^{a}e_{a}^{\mu }D_{\mu }\psi -m\psi =0.$$ Where $\gamma^a$ are a matrix representation of the Clifford algebra and these matrices satisfy the equation $\{ \gamma^a,\gamma^b \} = 2\eta^{ab}$, $e_{a}^{\mu }$ (a,b are indices in the tangent space) is the so called tetrad and it satisfies the equation $$g_{\mu\nu} = e^{i}_{\mu} e^{j}_{\nu} \eta_{ij},$$ and $D_\mu$ is the covariant derivative taking into account the spin of the particles, it is defined as: $$D_{\mu }=\partial _{\mu }-{\frac  {i}{4}}\omega _{{\mu }}^{{ab}}\sigma _{{ab}},$$ where $\omega _{{\mu }}^{{ab}}$ is the spin connection symbol and $\sigma _{{ab}}={\frac  {i}{2}}\left[\gamma _{{a}},\gamma _{{b}}\right],$ is the commutator of the Dirac matrices. One way to define the spin connection symbols is through the Christoffel symbols and the derivative of the metric as follows $${\displaystyle \omega _{\mu }^{\ ab}=e_{\nu }^{\ a} \left(\Gamma _{\ \sigma \mu }^{\nu }e^{\sigma b}+\partial _{\mu }e^{\nu b} \right)=e_{\nu }^{\ a}\Gamma _{\ \sigma \mu }^{\nu }e^{\sigma b}-e^{\nu a}\partial _{\mu }e_{\nu }^{\ b}}$$

Now let’s consider the cases of $+$ and $\times $ polarization separately. Further, let’s take the Weyl (Chiral) representation of the Gamma matrices, these matrices will be $${\displaystyle \gamma ^{0}={\begin{pmatrix}0&-I_{2}\\-I_{2}&0\end{pmatrix}},\quad \gamma ^{k}={\begin{pmatrix}0&\sigma ^{k}\\-\sigma ^{k}&0\end{pmatrix}},\quad \gamma ^{5}={\begin{pmatrix}I_{2}&0\\0&-I_{2}\end{pmatrix}},}$$ Where $\sigma^k$ are the standard Pauli matrices.

In addition, to simplify the equations we are going to consider massless particles. Plugging every component into Dirac’s equation in curved space for massless particles, $ i\gamma ^{a}e_{a}^{\mu }D_{\mu }\Psi = 0$, and keeping everything up to first order with respect of the gravity wave strain ($h$) we obtained the following systems of equations.\
For $+$ polarization, let’s say $h_{+} = \Re (a_{+} \mathrm{e}^{i\omega(z-t)}) = A_{+} \cos(\omega(z-t)) $, then the tetrad can be taken as : $$e^a_{\mu}\equiv\left(\begin{array}{cccc}-1&0&0&0\\0&\sqrt{1+h_{+}}&0&0\\0&0&\sqrt{1-h_{+}}&0\\0&0&0&1\end{array}\right)$$ The spin connections can be calculated and they are found as $$\begin{aligned}
\omega^{01}_{x} & =&  -\frac{1}{2} \left( \frac{\partial_t h_{+}}{\sqrt{1+h_{+}}} \right)    \    \     \     \    \    \   \    \   \omega^{13}_x = \frac12 \left(\frac{\partial_z h_{+}}{\sqrt{1+h_{+}}} \  \right)   \\
 \omega^{02}_{y} & = & -\frac12 \left( \frac{\partial_t h_{+}}{\sqrt{1-h_{+}}} \right)   \   \   \   \   \   \    \     \  \omega^{23}_{y} = \frac12 \left( \frac{\partial_z h_+}{\sqrt{1-h_+}} \right) \end{aligned}$$

The equations are then solved and separated as: $$\begin{aligned}
\bar{\sigma}^{\mu} \partial_{\mu} \psi_L & = & \sigma^2 \frac{h_+}{2} \partial_y \psi_L - \sigma^1 \frac{h_+}{2} \partial_x \psi_L \\
\sigma^{\mu}\partial_{\mu} \psi_R &=& \sigma^1 \frac{h_+}{2} \partial_x \psi_R - \sigma^2\frac{h_+}{2} \partial_y \psi_R\end{aligned}$$ where $\bar\sigma^{\mu}= (1,-\pmb{\sigma})$ and $\sigma^{\mu}=(1,\pmb{\sigma})$, and $\psi_L,\psi_R$ are the two component chiral spinors such that $\psi\equiv(\psi_L,\psi_R)$. As evident from the above equations the fermion wave has to have at least wave vectors in the $x$ or $y$ directions to produce a perturbation as observed in the scalar sector too.

Now, we take a perturbation $\psi=\psi_0+ A_{+} \tilde{\psi}$, and $\psi_0$ must satisfy the Weyl equation in flat space. These solutions are given by the right and left handed spinors (positive and negative helicity ), $\psi_L$ and $\psi_R$ respectively. If we denote $k^{\mu}_0 = \{\omega_0, \textbf{k}_0 \}$ and $x^\mu = \{ t, \textbf{x} \}$, then the plane wave solutions for the Weyl equation in flat space are given by

$$\begin{aligned}
\psi_{0L}(t,\textbf{x}) &= \frac{1}{\sqrt{2\omega_0(2\pi)^3}} \mathrm{e}^{ik^\mu_0 x_\mu} u^{(+)}(\textbf{k}_0) \\
\psi_{0R}(t,\textbf{x}) &= \frac{1}{\sqrt{2\omega_0(2\pi)^3}} \mathrm{e}^{ik^\mu_0 x_\mu} u^{(-)}(\textbf{k}_0),\end{aligned}$$

where the spinors $u^{(+)}(\textbf{k}_0)$ and $u^{(-)}(\textbf{k}_0)$ are determined by the eigenvalue equations $$\begin{aligned}
\hat{\sigma} \cdot \textbf{k}_0 u^{(+)}(\textbf{k}_0) &= -\omega_0 u^{(+)}(\textbf{k}_0) \\
\hat{\sigma} \cdot \textbf{k}_0 u^{(-)}(\textbf{k}_0) &= \omega_0 u^{(-)}(\textbf{k}_0)\end{aligned}$$

For the next example, let’s consider a particle moving in the $x$ direction. In this way $k^{\mu}_0 = \{\omega_0,\omega_0,0,0 \}$, $u^{(+)}(\textbf{k}_0) = \frac{1}{\sqrt{2}} \left( \begin{smallmatrix} 
&1   \\ \text{-}& 1  \end{smallmatrix} \right)$ and $u^{(-)}(\textbf{k}_0) = \frac{1}{\sqrt{2}} \left( \begin{smallmatrix} 
&1   \\ &1  \end{smallmatrix} \right)$, therefore

$$\psi_0 = \frac{1}{\sqrt{4\omega_0(2\pi)^3}}\left( \begin{smallmatrix}
&1  \\
\text{-}&1 \\
&1  \\
&1
\end{smallmatrix} \right) \mathrm{e}^{i\omega_0 (x-t)}
\label{eqn:solution}$$

Then, the perturbed equations are

$$\begin{aligned}
{\mathop{}\!\mathbin\Box}\tilde{\psi}_{L1} &= \omega_0^2e^{i \omega_0 (x-t)}\cos\omega(z-t)  \\
{\mathop{}\!\mathbin\Box}\tilde{\psi}_{L2 }&= e^{i \omega_0 (x-t)} \left(\omega_0^2 \cos (\omega  (z-t)) - i\omega\omega_0 \sin\omega(z-t) \right)\\
{\mathop{}\!\mathbin\Box}\tilde{\psi}_{R1} &= e^{i \omega_0(x-t)}\left(\omega_0^2 \cos(\omega(z-t)) + i \omega\omega_0 \sin(\omega(z-t))\right)\\
{\mathop{}\!\mathbin\Box}\tilde{\psi}_{R2} &=   \omega^2_0 e^{i \omega_0 (x-t)} \cos\omega  (z-t)\end{aligned}$$

We write the sines and cosines in exponentials and get forms like:

$${\mathop{}\!\mathbin\Box}\tilde{\psi}_{L2 }= \frac{\omega_0}{2}e^{i \omega_0 (x-t)} \left((\omega_0-\omega) e^{i \omega(z-t)}+ (\omega_0 +\omega)e^{-i \omega(z-t)}\right)$$

Therefore, the integral solutions of the above equations are obtained in the same way as a scalar wave perturbation. e.g.

$$\begin{aligned}
\tilde\psi_{L1+\omega} &= & \frac{ \pi \omega_0}{\omega}e^{i\omega_0(x-t)} e^{i\omega(z-t)}  \left[ e^{i (\omega+\omega_0) t} \cos \left(\sqrt{\omega_0^2 + \omega^2} ~t \right) + e^{i (\omega+\omega_0) t} \frac{i (\omega+\omega_0)}{\sqrt{\omega^2+\omega_0^2}}\sin\left(\sqrt{\omega^2 +\omega_0^2}~t\right) -1\right] \\ 
\tilde\psi_{L2+\omega} &=& \frac{\omega_0-\omega}{\omega_0}\psi_{L1+\omega}\end{aligned}$$

The complete $\tilde\psi_{L1(2)}=\tilde\psi_{L1(2)+\omega}+ \tilde\psi_{L1(2)-\omega}$.

As we don’t observe fermions by themselves, we build a Fermion bilinear: $\psi^{\dag}\gamma^0 \gamma^{\mu}\psi$, of which the 0-th component represents a ‘probability’ density $\psi^{\dag}\psi$. This when calculated for let us say the left handed spinor gives a function which oscillates in time. This can be calculated as \^= (\_[0]{}\^+ A\_[+]{} \^)(\_0 + A\_[+]{} ) = \^\_0\_0 + A\_[+]{}\^\_0 + A\_[+]{} \_0\^ We compute the change in the Fermion density due to the perturbation as &= & A\_[+]{}

The unperturbed fermion has a constant probability density, but the perturbed fermion wave has an oscillatory phase. The profile of the bilinear appears in Figure 2. The amplitude of the oscillations are proportional to $A_+$ and thus would be of the order of $10^{-21}$, compared to the initial fermion flux.

If we take the cross polarization, then, the tetrad can be taken as $$e_{\mu a}= \left(\begin{array}{cccc}1&0&0&0\\0&\sqrt{\frac{1-h_{\times}}{2}}&\sqrt{\frac{1-h_{\times}}{2}}&0\\0 &\sqrt{\frac{1+h_{\times}}{2}}&\sqrt{\frac{1+h_{\times}}{2}}&0\\0&0&0&1\end{array}\right)$$

The spin connections are found as \^[12]{}\_[t]{}=-12 \_t h\_

\_[x]{}\^[01]{}= -1[2]{} \_t h\_           \_x\^[02]{}= -1[2]{} \_t h\_          \_[x]{}\^[13]{} = 1[2]{} \_z h\_         \_[x]{}\^[23]{} = 1[2]{} \_z h\_ \^[01]{}\_[y]{} = 1[2]{} \_t h\_               \^[02]{}\_[y]{}- 1[2]{} \_t h\_             \^[13]{}\_[y]{} =1[2]{} \_z h\_           \^[23]{}\_[y]{} =-1[2]{} \_z h\_ \_[z]{}\^[12]{}=-1[2]{} \_z h\_

The Dirac equation is separable for the Weyl components as previously for the $(\psi_L,\psi_R)$ components though with a redefinition of the coordinates $u=\frac{1}{\sqrt{2}}\left(y-x\right)$ and $v= \frac{1}{\sqrt{2}}\left(x+y\right)$ and additional terms due to the tilting of the frame in the tangent space; $$\begin{aligned}
\bar{\sigma}^{\mu} \partial_{\mu} \psi_L & = & -\frac{h_{\times}}{2} \sigma^1 \partial_v \psi_L - \frac{h_{\times}}{2}\sigma^2 \partial_u \psi_L -i \frac{1}{8} \partial_t h_{\times}\sigma^3 \psi_L - \frac{3}{8} i \partial_z h_{\times}\psi_L \\
\sigma^{\mu}\partial_{\mu} \psi_R &=&  \frac{h_{\times}}{2} \sigma^1\partial_v \psi_R + \frac{h_{\times}}{2} \sigma^2 \partial_u  \psi_R - i \frac{1}{8} \partial_t h_{\times} \sigma^3 \psi_R + \frac{3}{8} i \partial_z h \psi_R\end{aligned}$$

As previously, we can assume perturbative solutions and solve the decoupled inhomogeneous equations of the D’ Alembertians. However, in this particular case, there are resonance modes as when the frequency and the direction of the Fermions and the gravity waves coincide, there is still a non-zero inhomogeneous term. If we assume $\psi_{L(R)}= \psi^0_{L(R)}+ A_{\times}\tilde\psi_{L(R)}$ as previously, then the perturbation equations precisely yield equations which are of the form \_[L1]{} &= & ((z-t)) \^[0]{}\_[L1]{}\
\_[L2]{} &=& - i \^2  ((z-t))\^0\_[L2]{}\
\_[R1]{} &=&+ i \^2   ((z-t)) \^0\_[R1]{}\
\_[R2]{} &=& -   ((z-t)) \^0\_[R2]{}

We use the same solution for the unperturbed fermion as in (\[eqn:solution\]) with $x,y\rightarrow u,v$, this is due to the different tangent frame, and solve the perturbations. These will have corrections of the order of $\omega/\omega_0$ of the $+$ polarized gravity waves. If we abinitio calculate the the equations for the example where the directions and frequency of the two interacting waves coincide the inhomogeneous terms are proportional to ${\mathop{}\!\mathbin\Box}\tilde\psi_{L1(L2)} = -i\omega^2/2(\exp(2i\omega(z-t)) \pm 1)$ and give rise to linear and quadratically growing perturbations which can be interpreted as resonant modes. The explicit perturbation solutions are \_[L1(L2)]{} \~( i + t e\^[2i (z-t)]{} i 2 \^2 t\^2)

As obvious the fermion density built from the bilinear will have runaway behaviour with time as expected in resonance. The physical implication of this can be realized in an experimental set up. Though, the perturbation approximation might break down, but definitely this behaviour is non-trivial. If we observe the bilinear it will be proportional to $A_+ \omega t$, so for $\omega t\sim 10^{19}-10^{20}$ the corrections will be fractions and thus can account for the neutrino background density. However a careful analysis of early universe cosmology is required for this.

Conclusion
==========

In this section we briefly summarize the results, and try to find the implications for real physical systems. As we showed in the section on scalar waves, perturbations on the scalar wave could range from $10^{-3}-10^{-13}$ depending on the frequency for the gravitational waves with which some primordial inflaton interacts. What this would imply for early universe cosmology is work in progress. Next we examine the implications of our results for the cosmic neutrino background density. To study the effect of the ‘scattering’ of the neutrino by the gravity wave, we have to compute a quantization of the scattered neutrino. Naively though, we can calculate the change in neutrino density as $\delta \rho$ (\[eqn:den\]) with the estimates for neutrino density in the CMB. We estimate that the flux of the neutrino will differ by a perturbation up to oscillatory factors. The corrections are proportional to $10^{-21}$ and thus their contribution to the effective degrees of freedom are very tiny, not in the range of the fractional corrections obtained in other papers [@neut1; @neut2]. From the cross polarization though, if the gravitational waves have high frequency in the early universe, fractional corrections might appear. At the time of neutrino decoupling the gravity wave frequencies are quite small [@nuovo], however, there can be some earlier inflationary epoch where the gravity wave frequencies were bigger. The relic of such an interaction would propagate with a perturbation and contribute to the non-thermal fluctuations of the background neutrino effective degrees of freedom at a later time [@adg].

With the closed-form solutions found in this paper we could compute (numerically) second order solutions and so on. It would be interesting to compare our analysis with the one made in [@pprod] and look for effects such as particle production and therefore energy lost. This could explain why gravity waves have a very small amplitude when observed on earth.

A more plausible use of the new calculation is for detection of gravity waves on earth in neutrino experiments. As obvious, neutrino flux gets perturbed by the gravity wave into oscillatory (note not in flavour) energy density, this is of the order of sensitivity of detectors built on earth. In particular the resonant mode of the fermion perturbation for cross polarization of the gravity wave has to be investigated further.

[**Acknowledgements:**]{} S.E.M.G. would like to thank MITACS Globalink 2017 for a summer internship fellowship. Majority of the work was completed in May-August 2017.

Appendix 1: Integral for the massless fields
============================================

$$\begin{aligned}
I &= \int_{\bar{B}(0,t)} \mathrm{e}^{i(\textbf{k}_0 \cdot (\textbf{x} + \textbf{r}) - \omega_0(t-r)+\delta_0)} \frac{\mathrm{e}^{i(\omega(z' - (t-r)) + \delta)}+\mathrm{e}^{-i(\omega(z' - (t-r)) + \delta)}}{2r} r^2\sin\theta d\phi d\theta dr \end{aligned}$$

We can calculate the integral without changing the axis and just using that $z' = z + r \cos \theta$, as in Equation (\[eqn:define\]) then the integral would be

$$\begin{aligned}
I &= \frac{\mathrm{e}^{i(\textbf{k}_0 \cdot\textbf{x}- \omega_0 t+\delta_0 )}}{2} \int_{\bar{B}(0,t)} \mathrm{e}^{i( \textbf{k}_0 \cdot \textbf{r} +\omega_0 r)} \left( \mathrm{e}^{i(\omega(r \cos\theta +r) + \omega(z-t)+\delta)} +\mathrm{e}^{-i(\omega(r \cos\theta +r) + \omega(z-t)+\delta)}\right) r \sin \theta d\phi d\theta dr \\
&= \frac{\mathrm{e}^{i(\textbf{k}_0 \cdot\textbf{x}- \omega_0 t+\delta_0 )}}{2} \int_{0}^{t} \int_{0}^{\pi}  \left( \mathrm{e}^{i(\omega(z-t)+\delta)}\mathrm{e}^{i(\omega+\omega_0)r} \mathrm{e}^{i\omega r \cos\theta} +  \mathrm{e}^{-i(\omega(z-t)+\delta)}\mathrm{e}^{i(\omega_0 -\omega)r} \mathrm{e}^{-i\omega r \cos\theta}\right) \left( \int_{0}^{2\pi} \mathrm{e}^{i\textbf{k}_0 \cdot \textbf{r}} d\phi \right) r \sin\theta d\theta dr\end{aligned}$$

Since the only part dependant of $\phi$ is $\textbf{k}_0 \cdot \textbf{r} = k_{0x} r \sin\theta \cos \phi + k_{0y} r \sin\theta \sin \phi + k_{0z}r\cos\theta$, the integral can be easily calculated again as a Bessel function of the first kind. Breaking again the integral into the $+$ and $-$ terms of the exponential form of the $\cos(\omega(z-t)+\delta)$, the first integral would be

$$\begin{aligned}
I_1 &= 2\pi\mathrm{e}^{i(\omega(z - t) + \delta)} \int_{0}^{t} \int_{0}^{\pi} \mathrm{e}^{i(\omega+\omega_0)r} \mathrm{e}^{i(\omega  +k_{0z} )r\cos\theta} J_0 \left( \sqrt{k_{0x}^2 + k_{0y}^2} r \sin\theta\right) r\sin\theta dr \\
&= 2\pi\mathrm{e}^{i(\omega(z - t) + \delta)} \int_{0}^{t}  \mathrm{e}^{i(\omega+\omega_0)r} \left( \int_{-1}^{
+1} \mathrm{e}^{i(\omega  +k_{0z} )r\cos\theta} J_0 \left( \sqrt{k_{0x}^2 + k_{0y}^2} r \sin\theta\right) d(\cos\theta)\right) r dr \\
&= 2\pi \mathrm{e}^{i(\omega(z - t) + \delta)} \int_{0}^{t}  \mathrm{e}^{i(\omega+\omega_0)r} \left( \int_{-1}^{+1} 2\cos((\omega  +k_{0z} )rx) J_0 \left( \sqrt{k_{0x}^2 + k_{0y}^2} r \sqrt{1-x^2}\right) dx\right) r dr \\
&\mbox{and since $k_{0x}^2+k_{0y}^2 + (k_{0z}+\omega)^2 = k_{0x}^2+k_{0y}^2 + k_{0z}^2+\omega^2 + 2k_{0z}\omega = \omega_0^2 + \omega^2 + 2k_{0z}\omega  $  }  \\
&\mbox{$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ = \omega_{0}^2 +\omega^2  + 2\textbf{k}\cdot \textbf{k}_0 =\lVert \textbf{k}+ \textbf{k}_0 \rVert  $, then we have} \\
&= 2\pi\mathrm{e}^{i(\omega(z - t) + \delta)} \int_{0}^{t} \mathrm{e}^{i(\omega+\omega_0)r} \left( 2\frac{\sin[\lVert \textbf{k}+ \textbf{k}_0 \rVert r]}{\lVert \textbf{k}+ \textbf{k}_0 \rVert r}\right)rdr \\
&=\frac{4\pi \mathrm{e}^{i(\omega(z - t) + \delta)}}{\lVert \textbf{k} + \textbf{k}_0\rVert} \int_{0}^t  \mathrm{e}^{i r(\omega + \omega_0)} \sin(r\lVert \textbf{k} + \textbf{k}_0\rVert) dr \\
&= \frac{4\pi \mathrm{e}^{i(\omega(z - t) + \delta)}}{\lVert \textbf{k} + \textbf{k}_0\rVert} \frac{ \left(-\lVert \textbf{k} + \textbf{k}_0\rVert+e^{i t \left(\omega +\omega _0\right)} \left(\lVert \textbf {k} + \textbf{k}_0\rVert \cos (\lVert \textbf{k} + \textbf{k}_0\rVert t)-i \left(\omega +\omega _0\right) \sin
   (\lVert \textbf{k} + \textbf{k}_0\rVert t)\right)\right)}{\left(-\lVert \textbf{k} + \textbf{k}_0\rVert+\omega +\omega _0\right) \left(\lVert \textbf{k} + \textbf{k}_0\rVert+\omega +\omega _0\right)} \\
&= \frac{4\pi \mathrm{e}^{i(\omega(z - t) + \delta)}}{\lVert \textbf{k} + \textbf{k}_0\rVert} \frac{ \left(-\lVert \textbf{k} + \textbf{k}_0\rVert+e^{i t \left(\omega +\omega _0\right)} \left(\lVert \textbf{k} + \textbf{k}_0\rVert \cos (\lVert \textbf{k} + \textbf{k}_0\rVert t)-i \left(\omega +\omega _0\right) \sin
   (\lVert \textbf{k} + \textbf{k}_0\rVert t)\right)\right)}{2\omega \omega_0 (1-\cos\beta)}\end{aligned}$$

Appendix 2: Integral for the massive scalar
===========================================

Isolating the $t'$ dependent terms in the integral (\[eqn:mass\]) for the massive perturbation we get \_[-]{}\^[(t-r)]{} e\^[-i (+\_0) t’]{} dt’ + (-) Writing $|\textbf{x}-\textbf{x'}| =r $ and $t-t'= \tau$, and observing that $t-t'> r$, the integral is e\^[-i (+\_0)t]{} \_[r]{}\^ (-i(+ \_0))   d+ (-) In the above, we keep only the $\cos(\tilde\omega t)$ as the inclusion of the imaginary part of the exponential, the sin function gives a spurious pole at $2 m^2$ in the final answer. = - (- i t)   (2( - ) )(2 (+ ) ) \[eqn:time\] we introduce $\Lambda$ to encode the above =2r (-i t) (r, ) where $\tilde{\omega}=\omega+\omega_0$ The remaining exponential integrals follow the same steps as in equations (81-87). If we label the expression in equation(\[eqn:time\]) as $\Lambda$ then in the end the perturbation is proportional to 2   (k\_x\^2-k\_y\^2) e\^[-i t]{} e\^[i (**k**+ **k**\_0)**x**]{} \_0\^t (r,) dr given as && 4(k\_[0x]{}\^2-k\_[0y]{}\^2) Similarly the correction in the transverse sector can be calculated.

[99]{} B. P. Abbott et al Phys. Rev. Lett [**116**]{} 061102 (2016). B. S. Sathyaprakash, B. F. Schutz, Living Reviews in Relativity (2009). M. Maggiore Physics Reports [**331**]{} 283 (2000).. P. Jones, et al Physics Letters [**B 781**]{} 621(2017); arXiv:1706.09402. V. Villalba, Physics Letters A [**136**]{} 197 (1989). L. Husdal, Galaxies, [**4**]{} 4, 78 (2016). G. Mangano et al, Nucl, Phys, [**B 729**]{} 221 (2005). Work in progress E. E. Flanagan, S. A. Hughes, New. J. Phys. [**7**]{} 204 (2005). A. V. Bitsadze [*Partial Differential Equations*]{} World Scientific (1994). E. Zauderer, [*Partial Differential Equations of Applied Mathematics,*]{} Wiley Interscience (2006). I. S. Gradshteyn, I. M. Rhyzik, [*Table of Integrals, Series and Products*]{} Academic Press, Elsevier, (2007). M. C. Guzetti et al, Riv. D. Nuovo. Cim. [**39**]{} (9) 399 (2016). P. Jones, et al Phys. Rev. [**D 95**]{}, 065010 (2017); arXiv:1610.02973.

[^1]: Here we have taken the peak frequency of the CMB as an example; one component of the EM field satisfies similar inhomogeneous wave equation; (work to appear)
