\section{Introduction}
It was just 100 years ago in Prague when Einstein wrote the paper \cite{Ei} in which he, for the first time, expressed his understanding of Mach's Principle. Within his pre-general relativity  theory in which there was only one metric function he considered a mass point inside a shell accelerated ``upwards" and found that the mass point is dragged along by the shell.

Many formulations and studies of Mach's Principle have appeared during the last 100 years, most of them were analyzed in the T\"ubingen conference in 1993 which led to the remarkable volume \cite{BP} containing lectures as well as valuable discussions. We studied Machian effects in various contexts, both in asymptotically flat spacetimes
 and within cosmological perturbation  theory -- see, e.g., \cite{BKL07}, and number of references therein; later, cf. Schmid \cite{CSchmid}.

More recently, we investigated the subtle question whether dragging of inertial frames  should also be attributed also to gravitational waves. After the discovery of binary pulsars losing energy and angular momentum as a consequence of emitting
 gravitational radiation it would be surprising if gravitational waves did not have an influence on local inertial frames. However, there are still doubts about the status of gravitational stress-energy as compared with stress-energy tensor $T_{\mu\nu}$ of matter in relation to Machian ideas (see, e.g., \cite{BP}, p. 83).

 In \cite{BKL08cqg1} and \cite{BKL08cqg2}, we studied gravitational waves with translational symmetry that rotate about the axis of cylindrical polar coordinates. We solved the Einstein equations to first order in the wave amplitude and formed a pulse. By solving the relevant Einstein equation  to second order we calculated and illustrated graphically the rotation of a local  inertial frame  at the axis which is caused by the rotating wave pulse.

 In the present work we investigate the effects of rotating gravitational waves  in a more general, asymptotically  flat setting, without assuming any symmetry. We again start out from linearized theory and construct an ingoing rotating pulse of radiation which later transforms into an outgoing pulse. While in the translational symmetric case our waves were characterized by just one harmonic index $m$ governing the number of wave crests in $\varphi$, now the situation becomes considerably richer involving both spherical harmonic indices $l$ and $m$.

 Our problem now resembles, in some respects, the study by Lindblom and Brill \cite{LiBr} of the gravitational collapse of a slowly rotating, massive spherical shell freely falling under its own gravity. We reconsidered and extended their analysis in \cite{KLB98}. The construction of our waves starts out from (the linearized version of) a Bonnor-Weber-Wheeler time-symmetric, ingoing and outgoing pulse of a typical width $a$ and a typical amplitude $C$. One can make a pulse that is rotating, but still localized in the radial direction so that it resembles a  falling and rotating shell. Although we do not find the complete metric up to the second order in the amplitude, the resulting spacetime  will be asymptotically  flat at spatial infinity since the first-order perturbations  decay as $h_{\mu\nu}^{(1)}\thicksim r^{-(l+1)}$ there, $l$ being the multipolar index. Hence, one can define ``fixed stars" at infinity at rest in an asymptotically  Minkowski system given in asymptotically  spherical coordinates $t, r,
 \theta, \varphi $. The waves have mass-energy, a nonvanishing total angular momentum, and a vanishing total linear momentum. Since we consider weak waves it is plausible to assume that a global coordinate system can be extended over the whole spacetime. At infinity it becomes the rest frame of the fixed stars.

Near the origin, the first-order metric of our waves behaves as $r^l$, so the region around the origin will be very nearly flat for $l$ sufficiently large. When, however, a local inertial frame  is introduced at the origin, we find that its axes rotate with respect to  the lines $\varphi=\rm const$ of the global frame, i.e., with respect to fixed stars at infinity.

 Near the origin the congruence of the worldlines $\varphi=\rm const$ twists and the observers attached to these lines experience Euler acceleration proportional to $d_t \omega_0$, where $\omega_0$ is the angular velocity of the inertial frame  near the origin.
 The angular velocity  $\omega_0$ enters {\it the second-order odd-parity dipole $l=1$ perturbation } of the metric, $g_{t\varphi }^{(2)}= - \omega_0\; r^2\sin^2\theta$. (The centrifugal acceleration is higher order in the angular velocity.) The situation thus indeed resembles the interior of a collapsing slowly rotating shell -- see \cite{KLB98} where the vorticity of the lines $\varphi=\rm const$ is given in covariant form. In \cite{KLB98} we also calculated how the fixed  stars at infinity rotate with respect to  the inertial frame  at the origin by considering photons emitted radially inwards from the stars.

 In the present work we first consider gravitational waves  with odd parity in the Regge-Wheeler gauge in the linearized Einstein theory. The metric components $h_{\mu\nu}^{(1)}$ describing the waves can all be derived from one potential-type function satisfying the flat-space scalar wave equation. We choose this function in the form of a generalized Bonnor-Weber-Wheeler ingoing and outgoing pulse characterized by specific harmonic multipolar orders $l$ and $m, |m|\le l$. For $m\ne0$, the field represents the superposition of {\it rotating} waves which are ingoing from infinity towards the origin, stay regular all the time, and become outgoing waves again to decay at infinity.

 Turning to the second-order metric $h_{\mu\nu}^{(2)}$, we employ the second-order perturbative Einstein equations  which are well known to be again linear differential equations  that contain as ``source terms" expressions quadratic in the first-order perturbations  $h_{\mu\nu}^{(1)}$ and their first- and second-order derivatives (see, e.g., \cite{NiGP}, \cite{GNPP}). Since we look for the second-order rotational perturbations  we concentrate on odd-parity again. To determine the influence of gravitational waves  on the rotation of local inertial frames  at the axis of symmetry we do not need to solve the second-order Einstein equations  in general, because the high $l$ terms do not have significant amplitude at the origin. The terms that give the inertial frame near the origin are $l=1, m=0$, so we can average the ``source terms" over azimuth $\varphi $, and the problem becomes axisymmetric. This simplifies the calculations. In addition, an important result proved and discussed in \cite{NiGP}, \cite{GNPP}
states that the Regge-Wheeler gauge can be introduced in both first and second order, and it is unique. As emphasized in \cite{GNPP}, perturbations  in the Regge-Wheeler gauge can be considered to be gauge-invariant expressions.

 The angular velocity  $\omega_0$ of rotation of the local inertial frames at the origin, whose neighborhood remains as flat as we wish for all times by choosing very high $l$, is found by solving a simple linear elliptic equation  for the second-order dipole odd-parity perturbation  with a complicated source term given by the averaged component, $\langle R_{t \varphi}^{(2)}[h^{(1)}, h^{(1)}]\rangle $, of the Ricci tensor, quadratic in $h^{(1)}$ and its derivatives.

 A lengthy, complicated form of $\langle R_{t \varphi}^{(2)}\rangle $ allows us to find $\omega_0$ explicitly in a useful form only at particular values of time. However, by using computer algebra system we produce nice ``bell-shaped" curves showing how $\omega_0$ increases as the waves approach the origin and decreases as the waves recede back to infinity. From $\omega_0$  plotted for different azimuthal harmonic index $m$, with $l$ fixed, it is easily seen how an increase of the azimuthal angular~momentum  of the waves implies a stronger dragging of inertial frames  (gyroscopes) in nearly flat regions surrounded by the waves.

 We also calculate how the celestial sphere
will look when observed from the inertial frame  at the origin. This is influenced both by the passage of the stellar light through gravitational waves  and by the dragging of the inertial frame  due to their rotation. The resulting figures describing the apparent motion of fixed stars on the celestial sphere nicely illustrate the role of the multipolar harmonic indices $l$ and $m$ on the pattern of the motion.

The paper is organized as follows.
Since odd-parity perturbations  follow from a function satisfying the wave equation, in the next Sections II and III we study rotating scalar-field waves in Minkowski space. We evaluate their energy density and angular~momentum  density averaged over the azimuth $\varphi $ and in every direction, after averaging also over the latitude $\theta$.
We turn to rotating linearized gravitational waves  in Sections IV and V. Recalling some basic properties of the Regge-Wheeler formalism and tensor harmonics, and fixing our convention, we calculate the averaged energy density of the rotating linearized gravitational waves  and compare it with that of scalar waves. We find that for ``rapidly rotating" waves with $m\thicksim l$ the scalar field part dominates.

The core of the paper lies in Sections VI and VII. Here we consider the second-order perturbations. After expanding in tensor spherical harmonics, and taking the axially symmetric part, we find the second-order Einstein equations  for the radial parts of the second-order odd-parity metric components in the Regge-Wheeler gauge. For general $l$ we get linear wave-type equations  with source terms given by combinations of averaged components $\langle R_{\mu\nu}^{(2)}[h^{((1)}, h^{(1)}]\rangle $. For the dipole perturbations  we obtain an elliptic equation which we solve by variation of constants. The angular velocity   $\omega_0$ of an inertial frame  near the origin is determined and analyzed both analytically and numerically in Section VII. Here also the figures demonstrating the dragging of inertial frames
by the waves and the apparent motion of fixed stars as seen from the inertial frame at the origin are presented.
In the concluding remarks we discuss the global, instantaneous character of the dragging of inertial frames as it follows from our results and show how it can be seen in diverse situations, including cosmological perturbations. We also indicate some open questions.
Some lengthy expressions are written down in Appendices. There, also various integrals over products of Legendre functions and non-trivial radial integrals are evaluated.

\section{Rotating scalar waves}
 In linearized Einstein theory the basic Regge-Wheeler equation \cite{RW} for odd parity waves becomes the usual flat-space 1-dimensional wave equation  for the radial part. It is thus natural to first present rotating waves in Minkowski  space.

 Consider the wave equation
 \begin{equation}
 - \frac{\partial^2\psi}{\partial t^2}+ \Delta \psi=0
 \label{s2eq1}
 \end{equation}
in spherical coordinates. Assume
\begin{align}
\psi_{lm}&
\frac{1}{2} \left[Q_l(t,r) Y_{lm}(\theta,\varphi )+\overline Q_l(t,r) \overline Y_{lm}(\theta,\varphi ) \right],
\label{s2eq2}
\\
Y_{lm}&=
N^m_l P^m_l(\cos\theta)e^{im\varphi },~~~l\ge0,~~|m|\le l;
\label{s2eq3}
\end{align}
bar denotes complex conjugation, $Y_{lm}$ are scalar spherical harmonics satisfying the standard normalization and orthogonality conditions over a unit sphere (see, e.g., \cite{Ja}),
\begin{align}
N^m_l = \sqrt{\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!}}~.
\end{align}
The angular part of Eq.$\,$(\ref{s2eq1}) then separates and the radial part satisfies
\begin{equation}
-\ddot Q_l+Q''_l+\frac{2}{r}Q'_l - \frac{l(l+1)}{r^2}Q_l=0.
\label{s2eq4}
\end{equation}
Hereafter the dot denotes $\partial_t$, the prime $\partial_r$.

\begin{figure}[hb]
\begin{center}
\includegraphics[angle = 0, width=6.4cm,bb=159 133 1052 1193,clip=true]{tminus.eps}
\\
\includegraphics[angle = 0, width=6.4cm,bb=159 133 1052 1193,clip=true]{tplus.eps}
\end{center}
\caption{\label{Figure1} Snapshots of function $\psi_{lm}$ for $l=27, m=5$, and times $\tilde t=-6,-4,-2, 0$ (upper image) or $\tilde t=0, 2, 4, 6$ (lower image).
For high values of $l$ the wave is concentrated near a shell with radius $r^2=a^2+t^2$,
so we can combine waves at several moments into one figure.
The wave rotates anticlockwise around the $z$ axis. A careful observation reveals that the wave comes inwards in the form of a leading spiral, i.e.,
with the outside ahead of the inside (see ``white arms'' at $\tilde t=-2$ and $-4$); at $\tilde t=0$ the spiral structure
has changed to a cartwheel but rotation keeps the wave away from the origin; at $\tilde t>0$ the spiral becomes trailing
(see "white arms" at $\tilde t=2$ and $4$). At large $\tilde t$ the spiral becomes tighter (cf. snapshots at $\tilde t=6$).
}
\end{figure}

Assume that $Q_l(t,r)=\hat Q_l(\omega,r)e^{i\omega t}$ and substitute
\begin{equation}
\hat Q_l(\omega,r)=\frac{1}{\sqrt{r}}u_l(\omega,r).
\label{25}
\end{equation}
Equation (\ref{s2eq4}) is then transformed into
\begin{equation}
u''_l+\frac{1}{r}u'_l - \frac{(l+{\textstyle{\frac{1}{2} }\!})^2}{r^2}u_l+\omega^2u_l=0,
\label{26}
\end{equation}
which is the Bessel equation  with index $\nu=l+{\textstyle{\frac{1}{2} }\!}$. The independent solutions of (\ref{26}) are given in terms of spherical Bessel and Hankel functions $j_l(\omega r), n_l(\omega r)$ and $h_l(\omega r)$. We want solutions regular at the origin, so we construct our radial waves as superpositions of
\begin{equation}
j_l(\omega r)=\sqrt{\frac{\pi}{2\omega r}} J_{l+{\textstyle{\frac{1}{2} }\!}}(\omega r),
\label{27}
\end{equation}
$J_l$ is a Bessel function of the first kind. Inspired by the example of the Bonnor \cite{Bo}, Weber-Wheeler \cite{WW} cylindrical pulse and by our recent work \cite{BKL08cqg1}, \cite{BKL08cqg2}, in which time symmetric incoming and outgoing rotating waves are smooth and finite everywhere at all time, we consider, in spherical polar coordinates, the superposition
\begin{align}
\label{s2eq8}
Q_l(t,r)&=B_l  \sqrt{\frac{\pi a^3}{2r}}
\int_0^\infty(a\omega)^{l+{\textstyle{\frac{1}{2} }\!}}\;e^{-a\omega} \; e^{-i\omega t} J_{l+{\textstyle{\frac{1}{2} }\!}}(\omega r) \;d \omega,
\end{align}
where the real amplitude $B_l$ and the characteristic width of the pulse $a$ are constant.
It is convenient to introduce
\begin{equation}
\alpha=\alpha(t)=a+it.
\label{29}
\end{equation}

Employing then formula (8.6.5) in \cite{BE} (see also formula 6.623 in \cite{GR}), we find the integral in (\ref{s2eq8}) to yield
\begin{equation}
Q_l(t,r)= B_l2^ll! \; \frac{(r/a)^l}{[(\alpha^2+r^2)/a^2]^{l+1}}.
\label{s2eq10}
\end{equation}

To construct {\it rotating} waves coming in and then going out in a time-symmetric manner we just multiply by $Y_{lm}(\theta,\varphi )$ containing the factor $e^{im\varphi}$. Clearly, the field is the superposition of rotating waves: it is apparent that for each $\omega$ the wave contains the factor
\begin{equation}
e^{im\varphi }e^{-i\omega t}=e^{im(\varphi -\frac{\omega}{m}t)}
\label{s2eq12}
\end{equation}
so that the wave pattern rotates with the rate $\Omega_p=\omega/m$. The rotation and the pulse character of the wave as well as its regularity is easily seen when $\psi_{lm}$ is written explicitly in real terms:
\begin{equation}
\psi_{lm} = \frac{B_l\, 2^l\,l!\, N^m_l P^m_l(\cos\theta)~\tilde r^l}{[((1+\tilde r^2-\tilde t^2)^2+4\tilde t^2]^{(l+1)/2}}\cos\left[  m\varphi  - \lambda(t,r)  \right],
\label{s2eq13}
\end{equation}
where
\begin{equation}
\tilde r=\frac{r}{a}
~,~~~
\tilde t=\frac{t}{a}
\label{s2eq14}
\end{equation}
and
\begin{align}
\label{s2eq15}
\lambda(t,r) &= (l+1) \arctan\frac{2\tilde t}{1+{\tilde r}^2-{\tilde t}^2}
\\\nonumber
	&=  {\rm arg}~\left[{\tilde r}^2+(1+i{\tilde t})^2\right]^{(l+1)}.
\end{align}

Some snapshots of the rotating waves profiles are illustrated in Figure \ref{Figure1}.

\section{Energy and angular~momentum  in rotating scalar waves}
We start with some simple properties of scalar fields. The energy density of the scalar field can be derived from the scalar field Lagrangian density
\begin{equation}
 {\cal L}= - {\textstyle{\frac{1}{2} }\!}\sqrt{-\overline g}\,\overline g^{\mu\nu}\partial_\mu\psi\partial_\nu\psi,
\label{31}
\end{equation}
where we use the background  Minkowski  metric in spherical coordinates
\begin{equation}
\overline g_{\mu\nu}=\{1, -\overline\gamma_{kl}\}=\{1, -1, -r^2,-r^2\sin^2\theta\}.
\label{32}
\end{equation}
(Hereafter, the Greek indices $\lambda, \mu, \nu = 0, 1, 2, 3$, the Latin indices $i,j ,k = 1, 2, 3$; $l,m$ are reserved for harmonic indices associated with $Y_{lm}$.)
The energy density, the time-time component of the stress-energy-momentum  tensor
\begin{equation}
\sqrt{-\overline g}~T^{\mu\nu}=\frac{\delta{\cal L}}{\delta\overline g_{\mu\nu}},
\label{33}
\end{equation}
reads
\begin{equation}
\varepsilon_S={{\frac{1}{2} }\!}\!\left(\dot\psi^2+\overline\gamma^{jk}\overline\nabla_j\psi\overline\nabla_k\psi\right)\!,
\end{equation}
$\overline\nabla_k\psi $ are covariant derivatives associated with $\overline\gamma_{jk}$; for scalar $\psi$, $\overline\nabla_k\psi =\partial_k\psi$, so in spherical coordinates
\begin{equation}
\varepsilon_S={{\frac{1}{2} }\!}\!\left[  \dot\psi^2+\psi'^2+\frac{1}{r^2}(\partial_\theta\psi)^2 +  \frac{1}{r^2\sin^2\theta}(\partial_\varphi \psi)^2  \right]\!.
\label{s3eq5}
\end{equation}
The energy density satisfies the conservation law which can be written as
\begin{align}
\label{35}
&\frac{d}{dt}\int_V{\textstyle{\frac{1}{2} }\!}\!\left(\dot\psi^2+\overline\gamma^{jk}\overline\nabla_j\psi\overline\nabla_k\psi\right) dV
\\\nonumber
&=\int_V\overline\nabla^i(\dot\psi\overline\nabla_i\psi)dV=\int_S\dot\psi\overline\nabla_i\psi n^idS.
\end{align}
The momentum density $p_k= - \dot\psi\overline\nabla_k\psi$  of the scalar field obeys a conservation law as well:
\begin{align}
\label{36}
&\frac{d}{dt}\int_Vp_kdV = \frac{d}{dt}\int_V -\dot\psi\overline\nabla_k\psi dV
\\\nonumber
&= - \int_V\overline\nabla_i\Big{\{}\overline\nabla^i\psi\overline\nabla_k\psi
- {\textstyle{\frac{1}{2} }}\delta^i_k \left[(\overline\nabla_j\psi)(\overline\nabla^j\psi) -\dot\psi^2\right]\Big{\}}dV\nonumber
\\\nonumber
&= -\int_ST^i_kn_idS.
\end{align}

Similarly, the angular~momentum  in the field obeys the conservation law
\begin{equation}
\frac{d}{dt}\int_V\varepsilon_{ijk}x^jp^kdV= - \int_S\varepsilon_{ijk}x^jT^k_sn^sdS.
\label{37}
\end{equation}
The angular~momentum  flux is given by the components $\varepsilon_{ijk}x^jT^{ks}$ through the surface $S$. The angular~momentum  density along the $z$-axis is
\begin{equation}
j_S= - \dot\psi\partial_\varphi \psi.
\label{s3eq9}
\end{equation}

Let us calculate the energy  density and the angular~momentum  density for rotating scalar waves. Consider a single $(l,m)$ component, like in (\ref{s2eq2}), written now in the form
\begin{equation}
\psi_{lm}={{\frac{1}{2} }}N_l^m\left[  Q_{lm}(t,r)e^{im\varphi }+ \overline Q_{lm}(t,r)e^{-im\varphi } \right] P^m_l(\mu)~,
\label{s3eq10}
\end{equation}
with $\mu=\cos\theta$, $P_l^m(\mu)$ are real Legendre functions.
The energy density (\ref{s3eq5}) of the scalar field becomes
 \begin{align}
\label{s3eq11}
&\varepsilon_S=\frac{1}{4} (N_l^m)^2 \Bigg\{ \left(  \dot Q_{lm}\dot{\overline Q}_{lm} +  Q'_{lm} \overline Q'_{lm}     \right)(P^m_l)^2
\\\nonumber
&+\frac{Q_{lm} \overline Q_{lm}}{r^2}\left[   \frac{m^2}{1-\mu^2} (P^m_l)^2+ (1-\mu^2)({P^m_l}')^2 \right]\nonumber
\\\nonumber
&+\frac{1}{2}\Big[   \big( \dot Q_{lm}^2  +{Q'}^2_{lm}  \big)  e^{2im\varphi }
+\big(  \dot{\overline Q}_{lm}^2 +\overline{Q'}^2_{lm}  \big)  e^{-2im\varphi }  \Big] (P^m_l)^2\nonumber
\\\nonumber
&+\frac{1}{2r^2}\left(  Q^2_{lm}e^{2im\varphi } +\overline Q^2_{lm}e^{-2im\varphi }  \right) \times
\\\nonumber
&~~~~~~~~~~~~~~~~~~~\left[  - \frac{m^2}{1-\mu^2}  (P^m_l)^2+\textstyle{(1-\mu^2)} ({P^m_l}')^2 \right] \Bigg\},
 \end{align}
 where we write ${P^m_l}'$ for $(d/d\mu) P^m_l$. If we average over $\varphi $ terms, we obtain
 \begin{align}
&\langle \varepsilon_S\rangle = \frac{1}{2\pi}\int_0^{2\pi}\!\!\!\!\!\varepsilon_Sd\varphi
\\\nonumber
&~=\frac{1}{4} (N_l^m)^2 \Bigg\{ (\dot Q_{lm}\dot{\overline Q}_{lm} + Q'_{lm} \overline Q'_{lm})(P^m_l)^2
\\\nonumber
&+\frac{Q_{lm} \overline Q_{lm}}{r^2}
\left[  \frac{m^2}{1-\mu^2} (P^m_l)^2+\textstyle{(1-\mu^2)}({P^m_l}')^2 \right] \Bigg\}.
 \label{s3eq12}
 \end{align}

 If, in addition, we average over the latitude $\theta$ by taking ${\textstyle{\frac{1}{2} }\!}\int_{-1}^1\langle\epsilon_S\rangle d\mu$,  and use the integrals of $P'$'s evaluated in Appendix A, we arrive at the energy density,  averaged in every direction, as a function of $r$ in the form
 \begin{align}
\langle\!\langle\epsilon_S\rangle\!\rangle  = \frac{1}{16\pi} \left[
\dot Q_{lm} \dot{\overline Q}_{lm} +Q_{lm}' {\overline Q}_{lm}' +l (l+1) \frac{Q_{lm} \overline Q_{lm}}{r^2}
\right].
 \end{align}

The $z$-component of the angular momentum density is still simpler:
 \begin{align}
 j_S=\frac{im}{4}  (N_l^m)^2 \Big(\dot Q_{lm}\overline Q_{lm} - \dot{Q}_{lm} Q_{lm}e^{2im\varphi } \!- {\rm c.c.} \Big)(P^m_l)^2\!,
 \label{313}
 \end{align}
``c.c.'' denotes ``complex conjugation''.

So, after averaging over $\varphi $ and over $\theta$ (see Appendix A)  we get
\begin{equation}
\langle\!\langle j_S\rangle\!\rangle =
\frac{im}{16\pi}
\left(  \dot Q_{lm}\overline Q_{lm} - Q_{lm}\dot{\overline Q}_{lm}   \right)\!.
\label{314}
\end{equation}

\section{ Rotating linearized gravitational waves}
We consider first the linearized theory of gravity in a flat background  in spherical coordinates $(t,r,\theta,\varphi )$. It is natural to decompose the metric   perturbations  into tensor harmonics. These are given in terms of standard spherical harmonics (\ref{s2eq3}) and their derivatives, see, e.g., \cite{RW}, \cite{Ze}, \cite{NI}. For any symmetric second-rank covariant tensor $S_{\mu\nu}$ we define the expansion
\begin{align}
\label{s4eq1}
  &S_{\mu\nu}={\sum_{l,m}}\Big[  {\cal A}_{0lm}(t,r)a_{0lm\;\mu\nu} + {\cal A}_{1lm}(t,r)a_{1lm\;\mu\nu}
\\\nonumber   &+  {\cal A}_{lm}(t,r)a_{lm\;\mu\nu} +  {\cal B}_{0lm}(t,r)b_{0lm\;\mu\nu} +  {\cal B}_{lm}(t,r)b_{lm\;\mu\nu}
\\\nonumber   & + {\cal Q }_{0lm} (t,r)c_{0lm\,\mu\nu} +  {\cal Q }_{lm} (t,r)c_{lm\,\mu\nu}+{\cal D}_{lm}(t,r)d_{lm\,\mu\nu}
\\\nonumber   &+ {\cal G}_{lm}(t,r)g_{lm\;\mu\nu}+{\cal F}_{lm}(t,r)f_{lm\;\mu\nu}           \Big],
\end{align}
where $a_{0lm\mu\nu}(\theta, \varphi ),\, ...,\,f_{lm\mu\nu} (\theta, \varphi )$ are tensor harmonics; we follow the notation and normalizations used by \cite{NI}. The harmonics can be classified into even- and odd-parities -- those with even-parity transform by the parity $(-1)^l$ under the inversion, i.e., the transformation $(\theta, \varphi ) ~\rightarrow~(\pi - \theta, \pi+\varphi )$, in the same way as $Y_{lm}(\theta, \varphi )$ -- while odd-parity harmonics are defined by the parity $(-1)^{l+1}$. The harmonics are mutually orthogonal and normalized over the sphere of radius $r$ in the following sense:
\begin{equation}
\int_0^{2\pi}\!\!\!\!\int_0^\pi\!\! \overline X_{l'm'\mu\nu}\delta^{\mu\rho}\delta^{\nu\sigma}X_{lm\rho\sigma}\sin\theta \,d\theta d\varphi=\delta_{l'l}\delta_{m'm},
\label{42}
\end{equation}
where
\begin{equation}
\delta^{\nu\sigma}={\rm diag}\{1,1,1/r^2,1/(r\sin\theta)^2\},
\end{equation}
 $X$ stands for any of the harmonics $a, ... , f$. The orthogonality properties enable us to find the ``coefficients" ${\cal A}_{0lm}$, etc. in the decomposition (\ref{s4eq1}) of the general covariant tensor $S_{\mu\nu}$. For example, in the following section we shall need the coefficient ${\cal Q }_{0lm}(t,r)$ multiplying $c_{0lm}$.
Using the orthogonality property of the tensor harmonics and multiplying $S_{\mu\nu}$ by $\overline c_{0lm}$, we find
 \begin{equation}
 {\cal Q }_{0lm}(t,r)=\int_0^{2\pi}\!\!\!\!\int_0^\pi\delta^{\mu\rho}\delta^{\nu\sigma}S_{\mu\nu}\overline c_{0lm\,\rho\sigma}\sin\theta d\theta d\varphi .
 \label{s4eq4}
 \end{equation}
It is well known that the resulting field equations  for the radial factors may be discussed separately for odd- and even-parity perturbations of the Schwarzschild  black hole, see, for instance, \cite{RW}, \cite{GNPP}, \cite{Bi}, \cite{NI}, which of course remains true for flat space. For our purpose it is sufficient to consider odd-parity waves. In the expansions (\ref{s4eq1}) the odd-parity harmonics are $c_{0lm}, c_{lm}$ and $d_{lm}$, with explicit forms
of the non-zero components
being
\begin{align}
\label{s4eq5}
c_{0lm\;t\theta} &= \frac{r}{\sqrt{2l(l+1)}}\frac{1}{\sin\theta}\partial_\varphi  Y_{lm},
\\
\label{s4eq6}
c_{0lm\;t\varphi} &= -\frac{r}{\sqrt{2l(l+1)}}  \sin\theta\partial_\theta Y_{lm},
\\
c_{lm\;r\theta} &= \frac{ir}{\sqrt{2l(l+1)}}  \frac{1}{\sin\theta}\partial_\varphi  Y_{lm},
\\
c_{lm\;r\phi} &= - \frac{ir}{\sqrt{2l(l+1)}}   \sin\theta\partial_\theta Y_{lm},
\\
d_{lm\;\theta\theta} &= \frac{ir^2}{\sqrt{2l(l+1)(l-1)(l+2)}}  \frac{1}{\sin\theta}X_{lm},
\end{align}
\begin{align}
d_{lm\;\varphi\varphi} &= \frac{-ir^2}{\sqrt{2l(l+1)(l-1)(l+2)}}  \sin\theta X_{lm},
 \\
\label{s4eq11}
 d_{lm\;\theta\varphi} &= \frac{-ir^2}{\sqrt{2l(l+1)(l-1)(l+2)}}  \sin\theta W_{lm},
\end{align}
where
\begin{align}
\label{s4eq12}
X_{lm}&=2\partial_\varphi \left(\partial_\theta - \cot\theta\right) Y_{lm}~,
\\\nonumber
W_{lm}&=\left(\partial^2_\theta  - \cot\theta\partial_\theta - \frac{1}{\sin^2\theta}\partial^2_\varphi    \right) Y_{lm}.
\end{align}
For the explicit forms of the other tensor harmonics, see \cite{RW}, \cite{Ze}, \cite{Bi} and \cite{NI}. As mentioned above, our convention for harmonics is  that used in \cite{NI} which differs slightly from \cite{RW} or \cite{Ze}. However, following \cite{NI} we assume the expansion of the odd-parity metric perturbations  $h^{(i)}_{\mu\nu}$ ($(i)=(1),(2)$ denotes the first- and second-order perturbations) to have the form
\begin{align}
\label{s4eq13}
\nonumber
h^{(i)}_{\mu\nu}=\sum_{lm}\Bigg[&-\frac{\sqrt{2l(l+1)}}{r}  h^{(i)}_{0lm}(t,r)c_{0lm\,\mu\nu}
\\
&+ i \frac{\sqrt{2l(l+1)}}{r}  h^{(i)}_{1lm}(t,r)c_{lm\,\mu\nu}
\\\nonumber
&+i \frac{\sqrt{2l(l+1)(l-1)(l+2)}}{2r^2} h^{(i)}_{2lm}(t,r)d_{lm\,\mu\nu} \Bigg],
\end{align}
so that the factors multiplying the harmonics (\ref{s4eq5}-\ref{s4eq11}) cancel the multiplicative factors at $h^{(i)}_{0lm }, h^{(i)}_{1lm }, h^{(i)}_{2lm }$ in (\ref{s4eq13}). Hence the radial functions, the $h$'s, are identical in \cite{RW}, \cite{Ze}, \cite{Bi} and \cite{NI}.

In the literature time-dependent perturbations  are usually Fourier analyzed and perturbations  with even parity discussed. In
\cite{Bi} a complete account of the Hamiltonian, Regge-Wheeler-Zerilli and Newman-Penrose formalism for perturbations (with both parities) of a Reissner-Nordstr\"om black hole is given and the dipole $l=1$ perturbations  are also considered. Putting $M=e=0$ there and choosing the Regge-Wheeler gauge condition for odd-parity perturbations, $h^{(i)}_{2lm\mu\nu}=0$, we obtain the radial functions $ h^{(i)}_{0lm\mu\nu}(t,r)$ and $ h^{(i)}_{1lm\mu\nu}(t,r)$ in terms of a function $R$ satisfying the wave equation\footnote{One can get these equations  from the original paper \cite{RW} as well as after going back from their Fourier decompositions in $e^{-ikt}$ to the time-dependence and after correcting several misprints in their equations (22).} as follows (discarding labels):
\begin{align}
&h_0 = - \frac{1}{(l-1)(l+2)}(rR)'~,
\label{s4eq14}
\\
&h_1 =- \frac{r}{(l-1)(l+2)}\dot R~,
\\
&- \ddot R+R''-\frac{l(l+1)}{r^2}R=0.
\end{align}
We see that last equation can be put in the form of Eq. (\ref{s2eq4}) for the scalar field amplitude $Q_l(t,r)$ through the substitution
\begin{equation}
R=rq.
\label{48}
\end{equation}
So in terms of $q(t,r)$ the odd-parity metric perturbations are given by
\begin{align}
h_0&= - \frac{1}{(l-1)(l+2)}(r^2q)',
\label{s4eq18}
\\
h_1&= - \frac{1}{(l-1)(l+2)}r^2\dot q.
\label{s4eq19}
\end{align}
These formulas are valid for $l\ge 2$; in the linear theory in vacuum the dipole odd-parity perturbations ($l=1$) can be transformed away by a simple gauge transformation  \cite{Ze,Bi}.
Later, in Section VI, we shall substitute for $q(t,r)$ the solution given in (\ref{s2eq8}) and consider (\ref{s4eq18}-\ref{s4eq19})
to be the first-order perturbations  $h^{(1)}$ generating the second-order perturbations.

\section{ Energy in gravitational waves }
 In \cite{KLB06} we defined gravitational energy  on $x^0$=const slices in stationary asymptotically  flat spacetimes  by taking the difference between the total mass-energy and the total ``mechanical" energy arising from the matter stress tensor. Many formulas given in \cite{KLB06} are true in general, non-stationary spacetimes $\!\!$.

Consider the slice $x^0$=const and the metric in the $1+3$ decomposition associated with the Hamiltonian formulation (see, e.g., \cite{MTW})
\begin{align}
ds^2&=g_{00}dt^2+2g_{0k}dtdx^k+g_{ik}dx^idx^k
\label{s5eq1}
\\\nonumber
&=(w^2-W^2)dt^2+2W_kdtdx^k-\gamma_{ik}dx^idx^k,
\end{align}
 $i, k=1, 2, 3$, $w=\sqrt{g_{00}+W^2}$ and $W^2=\gamma_{ik}W^iW^k$ (cf. (2.36) in \cite{KLB06}). The energy of matter can be calculated using the Einstein constraint equation  for
 \begin{equation}
 2wG^{00}={\cal R} + K^2 - K^{ik}K_{ik},
 \label{52}
 \end{equation}
 where ${\cal R}$ is the scalar curvature of the 3-metric $\gamma_{ik}$ and the external curvature of $x^0$=const slices is
 \begin{equation}
K_{ik}=w^{-1}\left(\nabla_{(i}W_{k)}   -{\textstyle{\frac{1}{2} }\!} \partial_t\gamma_{ik}   \right)~~,~~K=\gamma^{ik}K_{ik},
\label{s5eq3}
 \end{equation}
 $\nabla_k$ is a $\gamma$-covariant derivative. The total (gravitational  and matter) mass-energy as measured at infinity is $M$. The gravitational  field energy on the $x^0$=const slice can be defined as
 \begin{equation}
 E_G=M -\frac{1}{2\kappa}\int({\cal R}+K^2-K^{ik}K_{ik})\sqrt{\gamma}d^3x~,
 \label{s5eq4}
 \end{equation}
$\kappa$ is Einstein's gravitational constant, $\gamma=\det\gamma_{ik}$, all volume integrals here  and below are taken over $x^0$=const.

It is often advantageous to use more general spatial coordinates  than Cartesian coordinates  at infinity, so assume that the metric of the background  Minkowski  spacetime  which at infinity is ``the limiting spacetime" to the physical spacetime, has the general form\footnote{Details on how suitable background metrics and mappings between backgrounds  and physical spacetime  may be introduced can be found in \cite{BK05}, \cite{KBL97}.}:
\begin{equation}
 d\overline s^2= dt^2 - \overline\gamma_{ik}dx^idx^k.
 \label{55}
\end{equation}
Denoting $\overline\gamma$-covariant derivatives by $\overline\nabla_k$ and introducing tensorial quantities
\begin{align}
\Delta^m_{~ik}&=\Gamma^m_{~ik} - {\overline\Gamma}^m_{~ik}
\\\nonumber
&={\textstyle{\frac{1}{2} }\!}\gamma^{mn}\left( \overline\nabla_i\gamma_{nk} +\overline\nabla_k\gamma_{ni} - \overline\nabla_n\gamma_{ik}\right)\!,
\label{56}
\end{align}
one can write the scalar 3-curvature density entering (\ref{s5eq4}) as follows:
\begin{equation}
\frac{1}{2\kappa}\sqrt{\gamma}{\cal R}=\frac{1}{2\kappa}\sqrt{\gamma}{\cal L} - \frac{1}{2\kappa}\partial_i(\sqrt{\gamma}k^i).
\label{s5eq7}
\end{equation}
Here ${\cal L}$ is the Einstein Lagrangian
\begin{equation}
{\cal L}=-\gamma^{ik}\left(  \Delta^m_{~ik}\Delta^n_{~mn} - \Delta^m_{~in}\Delta^n_{~mk} \right),
\label{s5eq8}
\end{equation}
and $k^i$ are given by (cf. (3.16) in  \cite{KLB06}):
\begin{equation}
k^i=\gamma^{-1}\overline\nabla_k(\gamma\gamma^{ik}).
\label{s5eq9}
\end{equation}
In an asymptotically  flat spacetime  the integral of the second term on the right-hand side of (\ref{s5eq7}) gives again the mass,
\begin{equation}
- \frac{1}{2\kappa}\int\partial_i(\sqrt{\gamma}k^i)d^3x= - \frac{1}{2\kappa}\oint \sqrt{\gamma}k^i dS_i=M.
\label{s5eq10}
\end{equation}
This is the total mass of spacetime  as measured at spatial infinity. In our case of pure gravitational waves  the surface integral gives the total mass-energy of
the waves. In a general case of radiation going through spatial infinity for all times, time-dependent radiative terms $\thicksim 1/r$ would also contribute to the surface integral (\ref{s5eq10}), however for pulse-type radiation considered in this work only the {\it static} mass-energy results. Returning to $E_G$ in (\ref{s5eq4}), substituting from (\ref{s5eq7}) and (\ref{s5eq10}) we find the gravitational energy as a volume integral over the slice $x^0$=const:
\begin{align}
\label{s5eq11}
E_G&=\frac{1}{2\kappa}\int\left( K^{ik}K_{ik} - K^2-{\cal L}\right)\sqrt{\gamma}\;d^3x
\\\nonumber
&= \int\varepsilon_G\sqrt{\gamma}\;d^3x,
\end{align}
with ${\cal L}$ and $K_{ik}$ given in (\ref{s5eq8}) and (\ref{s5eq3}). The energy density $\varepsilon_G$ is a scalar density on $x^0$=const due to the use of the background metric\footnote{It can be shown that (\ref{s5eq11}) is gauge invariant provided that the rules of the Brill-Hartle averaging employed in the classical work of Isaacson \cite{Is68-2} are used (for example, under integrals divergences become small);  in the TT-gauge $\langle \varepsilon_G\rangle $ becomes $\thicksim \dot h_{ik}\dot h^{ik}$ of Isaacson.}.

Assume now the physical metric $\gamma_{ik}$ to be close to the background  metric $\overline\gamma_{ik}$
and $g_{00}=1+h_{00}$ so that
\begin{equation}
\gamma_{ik}=\overline\gamma_{ik}+h_{ik}~~~,~~~W_k=h_{0k}~~~,~~~w=1+\frac{1}{2}h_{00}.
\end{equation}
The terms proportional to ${\cal O}(h^2)$ will be neglected. The extrinsic curvature (\ref{s5eq3}) then becomes
\begin{equation}
K_{ik}=\overline\nabla_{(i}h_{k)0} - {{\frac{1}{2} }}\dot h_{ik},
\label{s5eq13}
\end{equation}
the Lagrangian (\ref{s5eq9}) is
\begin{align}
\label{s5eq14}
- {\cal L}= &{\textstyle{\frac{1}{2} }\!}\overline\nabla_mh\overline\nabla_kh^{km} - \textstyle{\frac{1}{4}}\overline\nabla_mh\overline\nabla^mh
\\\nonumber
& -{\textstyle{\frac{1}{2} }\!} \overline\nabla^i h^{km}\overline\nabla_mh_{ik}+ \textstyle{\frac{1}{4}} \overline\nabla^i h^{km}\overline\nabla_ih_{km}.
\end{align}

Since we are considering waves by using tensor spherical harmonics it is natural to choose the background  metric $\overline\gamma_{ik}$ in spherical coordinates:
\begin{equation}
\overline\gamma_{11}=1~~~,~~~\overline\gamma_{22}= r^2~~~,~~~\overline\gamma_{33}=r^2\sin^2\theta.
\end{equation}
We start from the metric for odd-parity perturbations  in the general form (\ref{s4eq13}) with harmonics given by (\ref{s4eq5}-\ref{s4eq12}). However, since in the odd-parity case only four (non-diagonal) components are non-vanishing, the following simple notation will turn out to be useful:
\begin{align}
\nonumber
h_{02}&=h_0f_2~,~h_{03}=h_0f_3~,~h_{12}=h_1f_2~,~h_{13}=h_1f_3~,
\\
f_2 & = -\frac{1}{\sin\theta}\partial_\varphi  Y_{lm}~,~f_3=\sin\theta\partial_\theta Y_{lm}.
\label{s5eq16}
\end{align}
Starting from these expressions, calculating $K_{ik}$ and ${\cal L}$ from (\ref{s5eq13}), (\ref{s5eq14}) using $\overline\nabla_k$ in flat-space spherical coordinates, we arrive, after straightforward - though not short - calculations at the following expression for the gravitational energy  density of odd-parity waves
\begin{equation}
\varepsilon_G=\frac{1}{2\kappa}\left( K^{ik}K_{ik}-K^2 - {\cal L}\right)\!,
\end{equation}
where
\begin{align}
\nonumber
K^{ik}K_{ik}&-K^2=
\frac{h_0^2}{2r^4\sin^2\theta}\Big[ \left(    \partial_\varphi  f_2+\partial_\theta f_3-2\cot\theta f_3   \right)^2
\\
&+\frac{1}{2r^4}  \left( r\dot h_1 - rh_0' +2h_0\right)^2 \left( f_2^2+\frac{1}{\sin^2\theta}f_3^2\right)\nonumber
\\
& - 4\partial_\theta f_2\left(   \sin\theta\cos\theta f_2 +\partial_\varphi  f_3  \right)\Big]\!,
\label{s5eq18}
\end{align}
and
\begin{align}
\nonumber
-{\cal L}=\frac{1}{2r^4}\Big[& 4h_1(rh'_1 - h_1)   \left( f_2^2+\frac{1}{\sin^2\theta}f_3^2\right)
\\
&  +\frac{h_1^2}{\sin^2\theta}\left(  \partial_\varphi  f_2 - \partial_\theta f_3     \right)^2 \Big]\!.
\label{s5eq19}
\end{align}

 We now express the radial parts of perturbations  from (\ref{s4eq18}-\ref{s4eq19}) in terms of the gauge invariant quantity $q$ satisfying the wave equation  (\ref{s2eq4})  ($Q$ there is denoted by $q$ in the gravitational case), and functions $f_2, f_3$ and their derivatives in terms of spherical harmonics according to (\ref{s5eq16}). The individual terms in (\ref{s5eq18}) to (\ref{s5eq19}) are processed as follows. Each term is quadratic in ``$hf$", say $h_0 f_2 h_1 f_3$. To get real results, we should take\footnote{\label{foononcom}Notice that the terms cannot be commuted since, for example, $Re(h_0f_2)Re(h_1f_3)\ne Re(h_0f_3)Re(h_1f_2)$ so that when evaluating (\ref{s5eq18}), (\ref{s5eq19}) the {\it corresponding terms} must be grouped together.}
\begin{align}
\nonumber
&Re(h_0f_2) Re(h_1f_3)={\textstyle{\frac{1}{2} }\!}(h_0f_2+\overline h_0\overline f_2) {\textstyle{\frac{1}{2} }\!}(h_1f_3+\overline h_1 f_3)
\\
&=\textstyle{\frac{1}{4}}(h_0f_2\overline h_1\overline  f_3+\overline h_0\overline f_2h_1f_3)+{\rm terms ~in}~e^{\pm2im\varphi }.
 \label{519}
\end{align}
The last terms drop out after averaging over $\varphi $ -- see the scalar-field case, equation   (\ref{s3eq11}) -- so we retain only the first terms with ($\textstyle{\frac{1}{4}}$) expressed in terms of $q_{lm}(t,r)$ and $P^m_l(\mu)$. In this way we obtain, after somewhat lengthy but straightforward calculations and arrangements, the $\varphi $-averaged energy density of odd-parity rotating gravitational waves:
\begin{align}
\nonumber
& 2\kappa\langle\varepsilon_G\rangle ={\frac{1}{4}}[l(l+1)]^2{\Bigg\{}  (\dot q\dot{\overline q}+q'\overline q' )P^2
\\\nonumber
+&\frac{q\overline q}{r^2} \left[    \frac{m^2}{1-\mu^2}P^2   +{(1-\mu^2)}{P'}^2     \right]                {\Bigg\}}
+\frac{[l(l+1)]^2}{2r^3}  ( r^2q\overline q)'    P^2
\\\nonumber
&+\frac{1}{2r}(r^2\dot q\dot{\overline q})' \left[   { \frac{m^2}{1-\mu^2}}P^2   +\textstyle{(1-\mu^2)}{P'}^2     \right]
\\ \nonumber
&+\left[ q'{\overline q}' +\frac{2}{r}(q{\overline q})' \right]\Bigg[  { \frac{m^2}{1-\mu^2}\left(    \frac{m^2+\mu^2}{1 - \mu^2} -l(l+1)           \right)  }P^2
\\
\nonumber
&+ \left({ \frac{4m^2}{1-\mu^2} }- l(l+1)        \right)\mu PP' +\textstyle{(\mu^2+m^2)}{P'}^2    \Bigg]
\\ \nonumber
&+\frac{q\overline q}{r^2}\Bigg{\{}\!{\frac{m^2}{1-\mu^2}}\!\left[  \!1\! - \!5l(l+1)\!+ \! { \frac{4(m^2+\mu^2)}{1-\mu^2}  }  \!\right] \!P^2
\\ \nonumber
&+\left(  \!{\frac{4m^2}{1-\mu^2}} \!-\textstyle{l(l+1)} \! \right)\!\mu PP'
\\
&+\left[{1\!+\!m^2-(1\!-\!\mu^2)l(l+1) }\right] {P'}^2 \Bigg{\}}. ~~~~~~~~~~
\end{align}
Here $P$ and $P'$ denote $P^m_l(\mu)$ and $d_\mu P^m_l(\mu)$; $q$ stands for $q_l(t,r)$.

This result has some noteworthy properties. The first terms in the first curly brackets coincide precisely with the energy (\ref{s3eq12})  in scalar waves averaged over $\varphi $. Together with the next term in $[l(l+1)]^2$ they dominate for big $l$. Nevertheless, for rotating waves when $m\thicksim l$ the scalar field part dominates. The very first term $\thicksim [l(l+1)]^2\dot q\dot{\overline q}$ will dominate for high frequency waves $(\dot q\dot{\overline q} \thicksim \omega^2)$ and high $l$ for both scalar and gravitational waves.

\section{ Second-order odd parity dipole perturbations}
In general the second-order metric perturbations  $h^{(2)}$  can be obtained by solving the equations
\begin{equation}
G_{\mu\nu}^{(1)}[{h^{(2)}}]= - G_{\mu\nu}^{(2)}[h^{(1)} ,h^{(1)}],
\label{s6eq1}
\end{equation}
where $h^{(i)}$ represent quantities given in (\ref{s4eq13}) with all indices seen explicitly.
The right-hand side   of (\ref{s6eq1}) is the source term in the form of an effective energy-momentum  tensor quadratic in the first-order perturbations  $h^{(1)}$ and their derivatives known from the solutions of $G_{\mu\nu}^{(1)}[h^{(1)}]=0$. In principle
writing down the right-hand side of Eq. (\ref{s6eq1}) and solving it is a clear task, in practise a formidable one.

However, to determine the influence of gravitational waves  on the rotation of local inertial frames  at the axis of symmetry we do not need to solve (\ref{s6eq1}) in general. First, following our method in the translational symmetric case \cite{BKL08cqg1}, \cite{BKL08cqg2}, we shall average over $\varphi $ so that the source term on the right-hand side  will become axisymmetric. Before averaging it will contain terms of the form $\partial_th^{(1)}\partial_\varphi  h^{(1)}$ each of which involves two terms independent of $\varphi $ and two terms varying like $\sin 2m\varphi  $ or $\cos 2m\varphi $. Those terms cannot cause any rotation of the inertial frames  on the axis. Hence, we shall concentrate on the terms that are independent of $\varphi $ and denote them again, as in previous sections, by the average symbol $\langle \rangle $. In addition, as in the translational symmetric case the dragging will be fully determined by solutions of just one equation  in (\ref{s6eq1}) -- that with $\mu=t, \nu=\varphi $.

Therefore, we consider the equation
\begin{equation}
G_{t\varphi }^{(1)}[h^{(2)}] = - \langle G_{t\varphi }^{(2)}[h^{(1)} ,h^{(1)}]\rangle,
\label{s6eq2}
\end{equation}
or equivalently,
\begin{equation}
R_{t\varphi }^{(1)}[h^{(2)}] = - \langle R_{t\varphi }^{(2)}[h^{(1)} ,h^{(1)}]\rangle,
\label{s6eq3}
\end{equation}
which follows from (\ref{s6eq2}) because in the background $\overline g_{\mu\nu}$ is the flat metric in spherical coordinates and $R_{\mu\nu}^{(1)}[h^{(1)}] =0$.
To solve these equations  we shall expand both sides in tensor spherical harmonics. Any symmetric second-rank covariant tensor can be expanded as in (\ref{s4eq1}). We are interested in the axisymmetric $(t,\varphi )$ component, so the only harmonic entering the left-hand side  of (\ref{s6eq3}) linear in $h^{(2)}$ is the odd-parity harmonic $c_{0lm}$ given in (\ref{s4eq6}).
In order to extract the contribution proportional to the same harmonic on the right-hand side  of (\ref{s6eq3})
 we imagine it to be expanded as the general covariant tensor $S_{\mu\nu}$ in (\ref{s4eq1}).
The coefficient ${\cal Q }_{0lm}(t,r)$ multiplying $c_{0lm\,\mu\nu}$ can be derived using formula (\ref{s4eq4}). In our case of the averaged  -- and thus axially symmetric -- component ~~~~~~~~~~$\langle R_{t\varphi }^{(2)}[h^{(1)},h^{(1)}]\rangle $ we put $m=0$ and $(\mu\nu)=(t\varphi )$ to obtain
\begin{align}
{\cal Q }_{0l0}(t,r)&= \int_0^{2\pi}\!\!\! \int_0^\pi -2\langle R_{t\varphi }^{(2)}\rangle \overline c_{0l0\;t\varphi } \frac{\sin\theta d\theta}{r^2\sin^2\theta} ~ d\varphi
\nonumber\\
&=\frac{4\pi}{\sqrt{2l(l+1)}r}\int_0^\pi\langle R_{t\varphi}^{(2)}\rangle \partial_\theta Y_{l0}\;d\theta ,
\label{65}
\end{align}
where we substituted for $\overline c_{0l0\;t\varphi }$ from (\ref{s4eq6}). Hence, the right-hand side  of equation  (\ref{s6eq3}) reads
\begin{equation}
2\pi\Bigg{\{}\!\!\int_0^\pi\!\!\!\!\langle R_{t\varphi }^{(2)}[h^{(1)}\!,h^{(1)}]\rangle \partial_\theta Y_{l0}d\theta\Bigg{\}}\left(   -\frac{\sin\theta\partial_\theta Y_{l0}}{l(l+1)}\right).
\label{s6eq6}
\end{equation}

The left-hand side  of (\ref{s6eq1}) in flat background  (in spherical coordinates) is formed from the second derivatives of the second-order perturbations:
\begin{align}
G_{\mu\nu}^{(1)}[h^{(2)}]= -{{\frac{1}{2} }\!} \Big[ &
  {h^{(2);\alpha}_{\mu\nu;\alpha}}
+ {h^{(2);\alpha}_{\mu\alpha;\nu}}
+ {h^{(2);\alpha}_{\nu\alpha;\mu}}  - h^{(2)\alpha}_{~\;\alpha\;;\mu\nu}
\nonumber\\
 & -\overline g_{\mu\nu} \!\left(   h^{(2)~;\alpha\beta}_{~\beta\alpha}  -  h_{~\beta~~~;\alpha}^{(2)\beta;\alpha}\right)   \Big].
\label{s6eq7}
\end{align}
We expand $G_{\mu\nu}^{(1)}$ in tensor spherical harmonics  and realize that we need only the axisymmetric component $G_{t\varphi }^{(1)}$; the result will thus be proportional to $c_{0l0}$. It is now important  that by using the  gauge transformation $x'\rightarrow x+\xi$ to the second order we can choose the Regge-Wheeler gauge also for the second-order perturbations  \cite{GNPP}, \cite{NI}, so for odd-parity perturbations  we can again set $h^{\,(2)}_2=0$ and keep just $h^{\,(2)}_0$  and $h^{\,(2)}_1$.

Calculations of $G_{t\varphi }^{(1)}[h^{(2)}]$ from (\ref{s6eq7}) lead us to the following expression:
\begin{equation}
 {{\frac{1}{2} }\!} \left[
 -{h_0^{(2)}}''\!\! +\frac{l(l+1)}{r^2}h^{(2)}_0\!
+{{\dot h}_1^{(2)}}{}'\!+\!\frac{2}{r}{\dot h^{(2)}}_1 \right]\sin\theta\partial_\theta Y_{l0}.
\end{equation}

Combining this with the right-hand side (\ref{s6eq6}) we arrive at the differential equation  implied by the field equation  (\ref{s6eq2})
for radial parts in the form
\begin{align}
\label{s6eq9}
&{{\frac{1}{2} }\!}\left[   {h_0^{(2)}}''\!\!\! - \frac{l(l+1)}{r^2} {h_0^{(2)}} -{{\dot h}_1^{(2)}}{}'\!-\frac{2}{r} {{\dot h}_1^{(2)}} \right]
\\
\nonumber
&~~~~~~~~~~~~~~~~ =\frac{2\pi}{l(l+1)}\int_0^\pi\!\!\!\!  \langle R_{t\varphi }^{\,(2)}[h^{(1)} ,h^{(1)}]\rangle \partial_\theta Y_{l0}  \,\, d\theta.
\end{align}
This equation  suggests that we need to consider more field equations  to determine both second-order perturbations
$h^{(2)}_0$  and $h^{(2)}_1\!\!$ from the first-order ones. Indeed, if we wish to calculate both
$h^{(2)}_0$  and $h^{(2)}_1$ for all $l$'s we have to take into account equations  (\ref{s6eq1}) also for $\textstyle{(\mu\nu)=(\theta,\varphi ),(r,\varphi )}$ and by combining them we can achieve just one wave-type equation  for $h^{(2)}_1$:
\begin{equation}
r \left(\partial_{t}^{2}-\partial_{r}^{2}+\frac{l(l+1)}{r^2}\right) \frac{h^{(2)}_{1}}{r} =S_{r\varphi }+S_{\theta\varphi }' - \frac{2}{r}S_{\theta\varphi },
\label{610}
\end{equation}
where components $S_{\mu\nu}$ stand for  the averaged components $\langle R_{\mu\nu}^{(2)}[h^{(1)} ,h^{(1)}]\rangle $. This is a hyperbolic equation  indicating that for general $l$ the effects of the first-order terms are perturbations  of the second-order terms that are non-instantaneous. However, in all our previous work on dragging of inertial frames  due to angular~momentum  of the sources (see, for example, \cite{BKL07} and references therein), the effects were instantaneous, at least at this lowest order. In fact, the same situation arises here. Inertial frames at the origin will be influenced primarily by the {\it dipole} perturbations  since all waves behave as $r^l$ there. Now it is well known that for $l=1$ one can achieve $h_1=0$ by an appropriate gauge transformation \cite{RW}, \cite{Bi}. For the dipole second-order perturbations  the equation  (\ref{s6eq9})
becomes
\begin{equation}
{h^{(2)}_{\;0}}'' -  \frac{2}{r^2}  h^{(2)}_{\;0} = g(t,r) ,
\label{s6eq11}
\end{equation}
where
\begin{equation}
g(t,r) = -
\sqrt{ \frac{3}{4\pi} }
\int_0^\pi\!\!
\int_0^{2\pi}\!\!\!\!\!\!R_{t\varphi }^{(2)}[h^{(1)} ,h^{(1)}]
\sin\theta \,d\theta d\varphi.
\label{s6eq12}
\end{equation}
Here we indicated the $\varphi $-averaging explicitly (which cancelled the factor $2\pi$) and we substituted for $\partial_\theta Y_{10}$. With $g(t,r)$
known, the solution of (\ref{s6eq11}) can easily be found by variation of constants from the simple solutions proportional to $1/r$ and $r^2$ of the homogeneous equation. Putting the arbitrary constants entering the homogeneous solution equal to zero, the complete solution for $h^{(2)}_{\;0}(t,r)$ reads
\begin{equation}
h^{(2)}_{\;0}= - \frac{1}{3r}\int_0^r \!\!g(t,r')\,r'^2 dr' - \frac{r^2}{3}\int_r^\infty\!\! g(t,r')\,\frac{dr'}{r'}.
\label{s6eq13}
\end{equation}

\section{How rotating waves drag inertial frames}
Since we are interested in the dragging of inertial frames  near the origin we wish to express
\begin{equation}
h^{(2)}_{0}(t,r)\simeq - \frac{r^2}{3}\int_0^\infty g(t,r')\frac{1}{r'}dr',
\end{equation}
with $g(t,r)$ given by (\ref{s6eq12}) and $\langle R^{(2)}_{t\varphi}\rangle$ in Appendix B.

The general form of the odd-parity metric perturbations  (\ref{s4eq13}), with the tensor spherical harmonics  (\ref{s4eq12}), for the special case of $l=1,m=0$ perturbations  with only $h^{\!^{(2)}}_{0\,10}\ne0$ yields the metric component $g_{t\varphi }$ near the origin in the form
\begin{equation}
g_{t\varphi }^{(2)}= - \sqrt{\frac{3}{4\pi}} h^{(2)}_{0}(t,r)\sin^2\theta = - \omega_0 \, r^2\sin^2\theta,
\end{equation}
where the rotation of an inertial frame  (of a gyroscope) located near the origin has the angular velocity
\begin{equation}
\omega_0=  \frac{1}{4 \pi} \int_0^\infty  \int R_{t\varphi }^{(2)}[h^{(1)} ,h^{(1)}]~ d\Omega~ \frac{dr}{r},
\label{Omega0}
\end{equation}
where $d\Omega = \sin \theta\, d\theta\, d\varphi$.

To evaluate $R_{t\varphi }^{(2)}[h^{(1)} ,h^{(1)}]$, which has more complicated structure than any expression calculated
so far, we found it more convenient to start from the real form of the solution of the original scalar wave equation (\ref{s2eq1}) and
use it to express the gravitational wave perturbations -- see Eqs. (\ref{s7eq4}), (\ref{s7eq5}) below.
As a consequence, also the metric coefficients will be real from the start, rather than of the
form (\ref{s5eq16}). With real expressions one does not have to be careful when grouping together properly
the corresponding terms in the quadratic expressions like we had to do in calculating the energy
densities (cf. footnote 4) and it is much easier to
find the results using the standard (commutative) computer algebra codes such as {\tt tensor} package in {Maple}.
In fact, we first proceeded by using the same methods as in the case of energies and, after going to the real
variables, we checked that both results agree.

So now we start from the metric perturbations written as follows
\begin{align}
      h_{t \theta}&=\frac{1}{\sin \theta} \frac{\partial \chi}{\partial r \partial \varphi}
~,    &h_{t \varphi}&=-{\sin\theta} \frac{\partial \chi}{\partial r \partial \theta}~,
\label{s7eq4}
\\
      h_{r \theta}&=\frac{1}{\sin\theta} \frac{\partial \chi}{\partial t \partial \varphi}
~,    &h_{r \varphi}&=-{\sin\theta} \frac{\partial \chi}{\partial t \partial \theta}~,
\label{s7eq5}
\end{align}
where function $\chi(t,r,\theta,\varphi) = r^2 \psi_{lm}$, ($\psi$ satisfies the scalar wave equation (\ref{s2eq1}) in flat
space); introducing $\chi$ instead of $\psi$ turns out to simplify the second-order Ricci tensor. Using the
general result for the the component $R^{(2)}_{t\varphi}$ (see (\ref{appBeq2}) in Appendix B) and substituting for $\chi$ from
(\ref{s7eq4}-\ref{s7eq5}) we arrive at
\begin{align}
\label{s7eq6}
&R^{(2)}_{t\varphi} = -\frac{1}{2}\, \frac {1}{{r}^{2}}\Bigg(
2\,\chi_{,\theta {r} {r} \theta}  \chi_{,\varphi t}-\chi_{,r \varphi r}  \chi_{,\theta t \theta}
\\\nonumber
&+2\,\chi_{,t r \varphi \theta}  \chi_{,r \theta}-\chi_{,{r} {r} \varphi \theta}  \chi_{,t \theta}-\chi_{,\theta t r \theta}  \chi_{,\varphi r}+\chi_{,{r} {r} \theta}  \chi_{,t \varphi \theta}
\\\nonumber
&
+\chi_{,\theta r \theta}  \chi_{,r \varphi t}
-\chi_{,t \theta}  \chi_{,{t} {t} \varphi \theta}-\chi_{,{t} {t} \theta}  \chi_{,t \varphi \theta}
-\chi_{,\varphi t}  \chi_{,\theta {t} {t} \theta} \Bigg)
\\\nonumber
&
+{\frac {\left (\chi_{,r \theta}  \chi_{,r \varphi t} -\chi_{,t r \theta}  \chi_{,\varphi r} -\chi_{,r \varphi r}  \chi_{,t \theta}
+\chi_{,{t} {t} \theta}  \chi_{,\varphi t}\right )\cos \theta}{2{r}^{2}\sin\theta}}
\\\nonumber
&-{\frac {\chi_{,r {\varphi} {\varphi} t}  \chi_{,\varphi r}
+\chi_{,{\varphi} {\varphi} r}  \chi_{,r \varphi t}-\chi_{,t \varphi t}  \chi_{,{\varphi} {\varphi} t} }{2{r}^{2} \sin^2\theta }}
+{\frac {\chi_{,r \theta} \chi_{,\varphi t}\cos\theta}{{r}^{3}\sin\theta}}
\\\nonumber
&
+{\frac {\chi_{,\theta r \theta}  \chi_{,\varphi t}}{{r}^{3}}}
+{\frac {\chi_{,\varphi t}}{2 {r}^{3}\sin^2\theta }}\left( 2\chi_{,{\varphi} {\varphi} r}+2r\chi_{,{\varphi} {\varphi} t t}  - r \chi_{,{\varphi} {\varphi} r r}\right)  ~.
\end{align}
We shall now write the real part of the scalar potential $\chi(t,r,\theta,\varphi)$ in the form corresponding to the form (\ref{s2eq13}) as the product
\begin{align}
\label{chi_fce}
 \chi
= \tilde B_l N^m_l\, \kappa(t,r) P_{l}^m(\cos\theta) \cos( m\varphi-\lambda(t,r) )~,
\end{align}
where $\tilde B_l=B_l\, 2^l l!$, $N^m_l$ and $\lambda$ are given in (\ref{s2eq14}) and (\ref{s2eq15}), and
we have introduced function
\begin{align}
 \label{s7eq8}
 \kappa(t,r) &=  \frac{\tilde r^{l+2} }{[ (1+\tilde r^2-\tilde t^2)^2 + 4 \tilde t^2 ]^{(l+1)/2} }~.
\end{align}
Inserting the function $\chi$ from (\ref{chi_fce}) into (\ref{s7eq6}) and employing the formulas for the integrals of the products
of the Legendre functions given in Appendix A to find the average over the spherical angles we first
find
\begin{align}
\label{s7eq9}
\nonumber
&\int R^{(2)}_{t \varphi}\, d\Omega = -
\tilde B_l^2 \frac{m}{r^2} \Big [
2\,\kappa\lambda_{{,r}}\kappa_{{,r t}}
-2\,{\kappa_{{,r}}}^{2}\lambda_{{,t}}
-2\,\kappa_{{,r}}\kappa\lambda_{{,r t}}
\\\nonumber
+&2\left ({l}^{2}+l+1\right )\kappa_{{,t}}\kappa_{{,r}}\lambda_{{,r}}
-\frac{2}{r}\,{{l\left (l+1\right )\kappa\left (\lambda_{{,r}}\kappa_{{,t}}-\kappa_{{,r}}\lambda_{{,t}}\right )}}
\\\nonumber
+&\left ({l}^{2}+l+2\right ) \kappa \left (\kappa_{{,t}}\lambda_{{,r r}}-\lambda_{{,t}}\kappa_{{,r r}}\right )
- l\left (l+1\right )\kappa^{2}{\lambda_{{,t}}}^{3}
\\
+&l\left (l+1\right )\left ( \kappa ^{2}\lambda_{{,t}}{\lambda_{{,r}}}^{2}
-2\,{\kappa_{{,t}}}^{2}\lambda_{{,t}}+\kappa\kappa_{{,t t}}\lambda_{
{t}}-\kappa\lambda_{{,t t}}\kappa_{{,t}}\right )\Big ].
\end{align}

\begin{figure}[b]
\includegraphics[angle = 270, width=8.4cm,clip=true,bb=65 104 546 744]{bell-1-10.eps}
\caption{\label{Figure2} Angular velocity of the central inertial frame $\omega_0(l,m;t)$ for $l=10$ and $m=1,2,... 10$ (from bottom to top).
The vertical axis is scaled in units of $\omega_0(10,1;0)$.}
\end{figure}
\begin{figure}[t]
\includegraphics[angle = 270, width=8.2cm, clip=true, bb=65 104 546 744]{bell2310etc.eps}
\caption{\label{Figure3} The dependence of normalized angular velocity of the central inertial frame $\omega_0(l,1;t)/\omega_0(l,1;t=0)$
on the parameter $l=2,$ $3, 10, 20, 30$ (from inside to out).
Also the function $(1+\tilde t^2)^{-3/2}$ is shown as a dashed line to indicate the limit for large $l$ (cf. Eq. (\ref{s7eq14})). }
\end{figure}

Substituting for $\kappa$ and $\lambda$ the expressions (\ref{s7eq8}) and (\ref{s2eq15}) and integrating over the radial
coordinate we obtain the angular velocity $\omega_0$ in the form
\begin{align}
\label{s7eq10}
 & \omega_0(t) =
\frac{\tilde B_l^2}{2\pi}\,\frac {m\left (l+1\right )\left (l+2\right )}{a\left (l+3\right )l}
\times
\\\nonumber
&~~~\left[\left( U_l -V_l {\tilde t\,}^2 \right)  \left (1+{\tilde t\,}^{2}\right ) I_{l+3}^2(\tilde t) + \left (U_l+V_l {\tilde t\,}^{2}\right ) I_{l+3}^1(\tilde t)\right],
\end{align}
where
\begin{align}
\label{s7eq11}
U_l &= \left (2\,{l}^{5}+7\,{l}^{4}+4\,{l}^{3}-7\,{l}^{2}+24\,l+36\right ),
\\\nonumber
V_l &= 3\left ({l}^{4}+2\,{l}^{3}+3\,{l}^{2}+8\,l+12\right )~.
\end{align}
Here the radial integrals $I_L^M$ are defined and analyzed in Appendix C.
Although the last result has indeed an analytic, explicit character, it is not very practical for
calculations because the integrals $I^{1,2}_{l+3}$ are determined by high-order derivatives for large $l$.
Hence, we describe their properties in detail in Appendix C.

Turning to the numerical evaluation of the angular velocity of dragging we obtain nice
``bell-shaped'' figures for $\omega_0$ for various values of $m$ at a given $l$. In Figure \ref{Figure2},
$\omega_0$ is plotted for $l=10$ and $m = 1,2,3, ... 10$.

Also, the integrals can be explicitly calculated at the moment of time-symmetry when
the dragging near the origin is maximal since the rotating waves are at nearest distance.
The integrals have simple forms also at $t = \pm a$, but these values are of no special
significance though they can serve for checking the numerical results. The central and maximal value reads
\begin{align}
\omega_0(t=0) =
\frac{\tilde B_l^2}{4\pi}\frac{m}{a}  {\frac {\left (l+1\right )!\,\left (l+2\right )!}{l \left (2\,l+3\right)!}} U_l~.
\end{align}

Finally, at late times the analytic behavior of $I^M_L$ implies
\begin{align}
\frac{\omega_0(t)}{\omega_0(t=0)} \sim  \left| \frac{a}{t} \right|^3  + O\left(\left| \frac{a}{t} \right|^5 \right)~,
\end{align}
and in Figure \ref{Figure3} we illustrate that $\omega(t)$ for large $l$ and given $m$ can be approximated by function
\begin{equation}
\label{s7eq14}
\omega_{l\rightarrow \infty} = \frac{ \omega_0(t=0) } {\left[1+\tilde t^2\right]^{\frac{3}{2}} } ~.
\end{equation}

\section{Observing stars through gravitational waves}
Since our waves are sufficiently weak near infinity, we can consider static stars there and ask how they would move across the sky as seen by an observer fixed in the flat region at the origin. We shall evaluate the first-order effects of the waves on the propagation of photons, the effects due to dragging of observer's frame are of second order and will be discussed subsequently. The influence of gravitational waves on the passing photons can be computed as a perturbation of the ingoing radial null geodesic
\begin{align}
 x^{(0)\mu} = [T + \lambda, -\lambda, \theta, \varphi], ~~\lambda \in (-\infty, 0]~,
\end{align}
which describes a ray from the star with celestial coordinates $\theta, \varphi$ which reaches the center at time $T$.
The geodesic equation can be written in the form
\begin{align}
 \frac{d}{d\lambda} x^\mu =& p^{(0)\mu} +  p^{(1)\mu} + ...~,
\\
 \frac{d}{d\lambda} p^{(1)\mu} =&
 - \Gamma^{(1)\mu}_{\rho\sigma}\left(x^\nu= x^{(0)\nu}(\lambda) \right) p^{(0)\rho}p^{(0)\sigma}
\\\nonumber
 &- 2\Gamma^{(0)\mu}_{\rho\sigma}\left(x^\nu= x^{(0)\nu}(\lambda) \right) p^{(1)\rho}p^{(0)\sigma}~.
\end{align}
To solve these differential equations we use covariant components of the perturbed photon momenta for which the effects of curvilinear coordinates cancel out; we omit superscript ${}^{(0)}$ for simplicity. Then we can write
\begin{align}
&\frac{d}{d\lambda} p^{(1)}_{\theta}
=
-r^2 \frac{d}{d\lambda} p^{(1) \theta}
- 2r p^{(1) \theta} p^r
\\\nonumber
&~~~~=
-({p^t})^{2}{h_{t\theta,t}}
-({p^r})^{2}{h_{r\theta,r}}
- {p^t}\,{p^r}\left[ {h_{r\theta,t}}+{h_{t\theta,r}} \right]
\\
&\frac{d}{d\lambda} p^{(1)}_{\varphi}
=
-r^2 \sin^2 \theta \frac{d}{d\lambda} p^{(1) \varphi}
-2 r p^{(1) \varphi} \sin^2 \theta p^r
\\\nonumber
&~~~~=
-({p^t})^{2}{h_{t\varphi,t}}
-({p^r})^{2}{h_{r\varphi,r}}
-{p^t}{p^r}\left[ {h_{r\varphi,t}}+ {h_{t\varphi,r}} \right].
\end{align}
\begin{figure}[hb]
\begin{center}
\includegraphics[angle = 270, width=4.2cm,clip=true,bb=64 172 541 650]{spir03.eps}
\includegraphics[angle = 270, width=4.2cm,clip=true,bb=64 172 541 650]{spir13.eps}
\\
~$l=3$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\;
$l=13$
\\
\includegraphics[angle = 270, width=4.2cm,clip=true,bb=64 172 541 650]{spir23.eps}
\includegraphics[angle = 270, width=4.2cm,clip=true,bb=64 172 541 650]{spir33.eps}
\\
~$l=23$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
$l=33$
\caption{\label{Figure4}
Since light from distant stars is influenced by the gravitational waves the observed positions of the stars change.
An observer at the origin can record the apparent position of the stars on the celestial sphere on a photographic plate.
When appropriately scaled and rotated, the trajectories of all stars are the same.
A star starts at the origin of the plate ($x=y=0$ in the planes above) and moves along closed trajectories the structure
of which becomes more complicated with increasing $l$. At the moment of time symmetry, $\tilde t=0$, the image is located
at maximal value of $x$ (with $y=0$).
This curve is an image of a straight line in the complex plane $z=1+is$, $s\in(-\infty,\infty)$ mapped by the function  $f(z) = z^{-l-2}$
(cf. Eq. (\ref{s8eq15})). Together with the trajectory, positions at time $\tilde t = 0,\pm 0.05,\pm 0.1,...$ are shown as circles.
}
\end{center}
\end{figure}

We see that along the radial rays (for which $\sin \theta^{(0)}=\rm const$) the covariant components of photons momenta $p^{(1)}_\theta$ and $p^{(1)}_\varphi$
can be written as integrals along the {\it un}perturbed ray and we do not need to solve the geodesic equation as set of differential equations.
Similarly, the differential equations for ${\theta}^{(1)}$ and ${\varphi}^{(1)}$ can be written as follows
\begin{align}
\frac{d}{d\lambda} {\theta}^{(1)}
&= -\frac{1}{r^2}   p^{(1)}_{\theta}
\\\nonumber
&= \frac{1}{p^{r}} \frac{d}{d\lambda} \left(\frac{1}{r}  p^{(1)}_{\theta} \right) - \frac{1}{p^{r}} \frac{1}{r} \frac{d}{d\lambda} p^{(1)}_{\theta}~,
\\
\frac{d}{d\lambda} {\varphi}^{(1)}
&= -\frac{1}{r^2\sin^2 \theta} p^{(1)}_{\varphi}
\\\nonumber
&= \frac{1}{p^{r}} \frac{d}{d\lambda} \left(\frac{1}{r\sin^2 \theta}  p^{(1)}_{\varphi} \right) - \frac{1}{p^{r}} \frac{1}{r\sin^2 \theta} \frac{d}{d\lambda} p^{(1)}_{\varphi}
~.
\end{align}
For photons which reach the origin the boundary terms arising after integration by parts vanish, so we can write
\begin{align}
\label{s8eq8}
\delta\theta(T) &= \left[ {\theta}^{(1)}\right]^{r=0}_{r=\infty}
= -\int_{-\infty}^0  \frac{1}{p^{r}} \frac{1}{r} \frac{d}{d\lambda} p^{(1)}_{\theta} d\lambda
\nonumber\\
=& -\int_{\infty}^0  \frac{1}{(p^{r})^2} \frac{1}{r} \frac{d}{d\lambda} p^{(1)}_{\theta} dr
\\\nonumber
=& \int_0^\infty \left[{h_{r\theta,t}}+{h_{t\theta,r} - {h_{t\theta,t}} -{h_{r\theta,r}} } \right]_{| t=T-r} \frac{dr}{r}
~,
\end{align}
\begin{align}
\label{s8eq9}
\nonumber
\delta\varphi(T) &= \left[ {\varphi}^{(1)}\right]^{r=0}_{r=\infty}
= -\int_{-\infty}^0 \frac{1}{p^{r}} \frac{1}{r\sin^2 \theta} \frac{d}{d\lambda} p^{(1)}_{\varphi} d\lambda
\\
=& -\int_{\infty}^0 \frac{1}{(p^{r})^2} \frac{1}{r\sin^2 \theta} \frac{d}{d\lambda} p^{(1)}_{\varphi} dr
\\\nonumber
=& \int_0^\infty \!\!\!\left[  {h_{r\varphi,t}}+ {h_{t\varphi,r}}  -{h_{t\varphi,t}} -{h_{r\varphi,r}} \right]_{| t=T-r}\frac{dr}{r\sin^2 \theta}
~.
\end{align}
Since now we work in the linear approximation  to find the light deflection we can go back to the simpler, complex-valued formulas (\ref{s2eq10}) and (\ref{s4eq14}-\ref{s4eq19}).
When we introduce amplitudes
\begin{align}
\Delta\theta &=
 {B_l N^m_l} 2(l-1)! \frac{m P^m_l(\cos \theta) }{\sin \theta},
\\
\Delta\varphi &=
- {B_l N^m_l}2(l-1)! P'^m_{~l}(\cos \theta),
\end{align}
the integrals (\ref{s8eq8}) and (\ref{s8eq9}) simplify to
\begin{align}
\delta\theta(T) &= {Re}\big[\Delta\theta\; \rho(a+iT) \big],
\\
\delta\varphi(T) &= {Re}\big[\Delta\varphi\; i \rho(a+iT) \big],
\end{align}
where the complex function $\rho(a+iT)$ reads
\begin{align}
\rho&= -i \;  e^{im\varphi} \frac{l(l+1)(l+2)}{\left (a+ i T\right )^{l-1}}\int^{\infty}_0{\frac {a^{l+2}\;{(2r)}^{l-1}}{\left (a+ i T-2\, i r\right )^{l+3}}}~dr
\nonumber\\
&= -  e^{im\varphi} i^{l+1} \frac{1}{\left(1+i\frac{T}{a}\right)^{l+2}}~.
\end{align}

If we assume a telescope pointed towards the star's initial position at $T\rightarrow -\infty$,
the quantities $\delta \theta$ and $\sin \theta \delta \varphi$ are proportional to star's coordinates
on telescope's photographic plate; star's trajectory can be described by a simple formula
\begin{align}
  \label{s8eq15}
  \frac{\delta\varphi}{\Delta\varphi} +  i \frac{\delta\theta}{\Delta\theta} = \frac{i^l\;e^{im\varphi}}{\left(1+i \frac{T}{a} \right)^{l+2}}~.
\end{align}
We see that the same image is replicated for all stars (Figure \ref{Figure4}).
Depending on the star's position the image is rotated by $e^{ i m\varphi}$ and
then scaled by factors $\Delta\theta$ in the latitudinal and by $\Delta\varphi$ in the longitudinal direction. Surprisingly, the images are time-symmetric with the maximal deflection occurring at $T=0$ despite the fact that the deflection is calculated using retarded integrals (\ref{s8eq8}) and (\ref{s8eq9}).

Since the calculation of the second-order metric everywhere would be very complicated we cannot, of course, describe all second-order effects.
However we can describe the second-order effect of the total rotation of the central observer's telescope due to the dragging of his inertial frames given above, cf. (\ref{s7eq10}), (\ref{s7eq11}). Since the pulse starts and ends at infinity at $t=\pm \infty$, the total rotation of the telescope is clearly gauge invariant.
As a consequence, the curves in Figure \ref{Figure4} will no longer be closed.
The stellar initial position at $T \rightarrow -\infty$ will differ from the final position at $T \rightarrow +\infty$
by an amount proportional to the area under the graph of $\omega_0(t)$.

\section{Concluding remarks}
The results obtained here extend significantly our past studies \cite{BKL08cqg1,BKL08cqg2} of gravitational waves with the translational symmetry that rotate about the axis of cylindrical symmetry and cause a rotation of local inertial frames at the axis. In the present work we investigated the effects of the rotating gravitational waves in a general, asymptotically flat situation. The main conclusions following from these studies demonstrate explicitly that any statement of Mach's principle that attributes
all dragging of inertial frames solely to the distribution of energy and momentum of matter, characterized by $T_{\mu\nu}$, as the origin of inertia is false --
gravitomagnetic effects are also caused by gravitational waves satisfying Einstein's field equations in vacuum for which a local energy and momentum cannot in general be defined.

We found that the dragging of inertial frames at the origin occurs {\it instantaneously}, with no time-delay  depending on the position of the pulse. The dipole perturbations  determining the dragging follow from the constraint equation  which is elliptic (see equation  (\ref{s6eq11})). Gravity is a tensor field, so gravitational waves  cannot have a dipolar part. This is the case with our first-order rotating pulse. It is the non-linearity of Einstein's field equations  causing a self-interaction  of the waves which produces an effective source for the second-order perturbations. For a general multipole moment $l\ge 2$, the resulting equations  are again wave-type equations  indicating that the effects of the first-order terms cause the second-order terms in a non-instantaneous way. However, inertial frames  at the origin will only be influenced by the dipole perturbations  since waves decay as $r^l$ at the origin. And for the {\it dipole} $(l=1)$ second-order perturbations, the equation  becomes elliptical,
 with the source term given by a component of the second-order Ricci tensor $R_{t \varphi}^{(2)}$ quadratic in the first-order metric and its derivatives.

The instantaneous character of the dragging of inertial frames  by angular momentum  of sources implied by the Einstein constraint equation  can, in fact, be seen in diverse situations, including those in cosmological perturbation  theory. We encountered the effect in the study of local inertial frames  affected by the rotations beyond the cosmological horizon \cite{BLK1}, \cite{BLK2}, or in the general analysis of ``rotations and accelerations" of local inertial frames  in the linearly perturbed Friedmann-Robertson-Walker  spacetimes  \cite{BKL07}. In the cosmological context we introduced desirable gauges - we call them the ``Machian gauges" - as those in which local inertial frames  can be determined instantaneously from the distributions of energy  and momentum in the universe by means of the perturbed Einstein constraint field equations. The ``uniform-Hubble-expansion-gauge" or the ``minimal-shear hypersurface" gauge combined with ``the minimal distortion" gauge on the spatial metric are examples of the
Machian gauges; the synchronous gauge or the generalized Lorentz-de Donder (harmonic) gauge are not. (See \cite{BKL07} for detailed analysis and preferences of the Machian gauges.)

Likewise, in the context of gravitational waves  we should bear in mind that such physical effects like the dragging of inertial frames are of a global nature and require the introduction of suitable coordinates in which the waves are described\footnote{This point has been emphasized also in contributions and discussions in the excellent book \cite{BP} (see, for example, remarks by D. Brill, p. 213).}.

In the present work we also investigated for the first time, as far as we are aware, the effects of the rotating gravitational waves on the propagation of light from distant fixed stars at (asymptotically flat) infinity.

There are various open problems in these themes which deserve further study. In order to have more natural understanding of the effect of rotating gravitational waves one should consider the waves within a cosmological (not asymptotically flat) context. Although in \cite{BKL07} we emphasized the role of the instantaneous Machian gauges, it would be also interesting to study these effects, in particular those connected with gravitational waves, in other type of gauges like the
Lorenz-de Donder gauge (cf. Appendix D in \cite{BKL07}).

Of course, any simple model problem illustrating the Machian effects {\it exactly} within general relativity, like that employing the conformal stationary metrics in our recent note on the linear gravitational dragging \cite{cqg12}, deepens considerably our understanding of such effects.

\section*{Acknowledgments}
J.B. and J.K acknowledge the hospitality of the Institute of Astronomy in Cambridge. J.B. and T.L. acknowledge the partial support from the Grant GA CR 202/09/00772 of the Czech Republic and the Grant No MSM0021620860 of the Ministry of Education.


\section{Motivation}
\label{lbl:intro}
The axion is a hypothetical elementary particle, which emerged originally from a proposal by Peccei
and Quinn, intended to solve the strong CP problem \cite{src:WISPy} in theoretical physics.
The axion is neutral, only interacts weakly with matter, has a low mass ($\approx 10^{-4} eV/c^2$), spin zero, and a natural decay constant (to 2 photons) in the order of $10^{17}$ years.
The axion belongs to the family of Weakly Interacting Sub-eV Particles (WISP). Another WISP, closely related to the axion is the paraphoton or hidden photon.
The existence of these WISPs could not be confirmed yet and all experimental efforts to date have so far produced only exclusion results. Nevertheless there is strong motivation to advance the experimental ``low energy frontier'' as the axion is the most popular solution for the strong CP-problem. Many WISPs are also excellent candidates for dark matter and explain numerous astrophysical phenomena.

\section{Experimental setup}
WISPs can be probed in the laboratory by ``Light Shining through the Wall'' (LSW) experiments.
They exploit the very weak coupling to photons, allowing an indirect proof of the otherwise hidden particles without relying on any cosmological assumptions.
Previous LSW experiments have been carried out with optical laser light at DESY (ALPS), CERN (OSQAR) and Fermilab (GammeV).

The concept of an optical LSW experiment can be adapted to microwaves \cite{src:hoogUW,src:moiMykonos}. A block diagram of the setup is shown in Fig.~\ref{fig:ovrBlock}, it consists of two identical low loss microwave cavities with a diameter of 140 mm, a height of 120 mm and a spacing between them of 150 mm.
One serves as WISP emitter and is excited by an external microwave source. It develops a strong electromagnetic (EM) field, which corresponds to a large amount of microwave photons $\gamma$.
Theory predicts that some of these photons convert to paraphotons $\gamma'$ by kinetic mixing (similar to neutrino oscillations) or -- if the cavities are placed in a strong static magnetic field -- to axion-like particles by the Primakoff effect~\cite{src:WISPy}.
Both particles only interact very weakly with matter (similar to neutrinos in this respect) and thereby, in contrast to the photons, can traverse the cavity walls. Some WISPs propagate towards the detection cavity, which is connected to a very sensitive microwave receiver. The reciprocal conversion process transforms WISPs to microwave photons, which can be observed as an excitation of the seemingly empty and well shielded detection cavity.
\begin{figure}[t]
\begin{center}
	\includegraphics[width=0.49\textwidth]{ovr2.pdf}
	\caption{Block diagram of the experiment}
	\label{fig:ovrBlock}
\end{center}
\vspace{-0.75 cm}
\end{figure}
Since there is no energy loss associated with the WISP conversion process, the regenerated photons in the detecting cavity have exactly the same energy as the photons in the emitting cavity. Thus, the signal which is coupled out from the detection cavity has the same frequency as the one which is generated on the emitting side, making a narrowband receiving concept feasible.

This paper will focus on the latest exclusion results for \textbf{paraphotons} from the microwave WISP search at CERN. In a future upgrade, an additional magnet will allow the search for axions.

Considering current exclusion limits, it takes $> 10^{24}$ photons on the emitting side to generate one photon on the detection side, making this the most challenging aspect of an LSW experiment.
The expected output power (or photon flux) from the detecting cavity towards the microwave receiver due to paraphotons is given by Eq.~\ref{equ:power},
\begin{align}
	\label{equ:power}
	P_{\mathrm{det}} &= \chi^4 \left(\frac{m_{\gamma'} c^2}{f_{\mathrm{sys}} h}\right)^8 |G|^2 Q_{\mathrm{em}} Q_{\mathrm{det}} P_{\mathrm{em}}
\end{align}
where $Q_{\mathrm{em}}$ and $Q_{\mathrm{det}}$ are the loaded Q factors of emitting and detection cavity, $f_{\mathrm{sys}}$ is the frequency where the experiment is carried out (and to which the cavities are tuned),
$h$ is Planck's constant and $G$ is a dimensionless geometric form factor in the order of 1, describing the position, shape and resonating mode of the cavities \cite{src:JaCaRi}.
The rest mass of hidden photons is a priori unknown and given by $m_{\gamma'}$.
The kinetic mixing parameter $\chi$ describes the likeliness of paraphoton - photon oscillations. A previous examination of Coloumb's law indicates that $\chi < 3 \cdot 10^{-8}$ in this energy range.

If there is no significant signal detected, an exclusion result can be produced by
determining $\chi$ from the other known values. This provides a convenient way to compare the achieved sensitivity to other experiments.

The parameters of the paraphoton experiment as it has been set up and carried out at CERN in March 2012, are summarized in Table~\ref{tbl:param}. As no paraphotons were observed, the corresponding exclusion limit in  comparison to other experiments is shown in Fig.~\ref{fig:exclPlot}.
\begin{table}[t]
\centering
\caption{parameters of the paraphoton run in March 2012}
\label{tbl:param}
\begin{tabular}{c}
\toprule
$f_{\mathrm{sys}} = 2.9565$ GHz \quad $Q_{\mathrm{det}} = 23620$ \quad $Q_{\mathrm{em}} = 23416$\\[0.1 cm]
$P_{\mathrm{det}} = 8.51 \cdot 10^{-25}$ W \quad $P_{\mathrm{em}} = 37$ W \quad $|G| = 0.222$ \\
\bottomrule
\end{tabular}
\vspace{-0.3 cm}
\end{table}
\begin{figure}[b]
\vspace{-0.5 cm}
\begin{center}
	\includegraphics[width=0.49\textwidth]{exclusion.pdf}
	\caption{Exclusion limit for paraphotons as a result of the measurement-run at CERN in March 2012, compared to other experiments (details in \cite{src:JaCaRi})}
	\label{fig:exclPlot}
\end{center}
\vspace{-0.5 cm}
\end{figure}

\section{Engineering aspects}
On the left side of Fig.~\ref{fig:ovrBlock}, a commercial microwave source is shown, which generates a signal on $f_{\mathrm{sys}}$ (see Table~\ref{tbl:param}) which is amplified up to 50 W and drives the emitting cavity on its resonant frequency.

Power is coupled in and out of each cavity with a small inductive coupling loop, adjusted for critical coupling to the TE$_{011}$ mode. This mode has been chosen for its high Q-factor and reasonable coupling to paraphotons compared to other modes. The loaded Q-factor of the silver coated brass cavities has been determined by a network analyzer, their 3 dB bandwidth is $BW_{3 dB} \approx 126$~kHz. A tuning stub perturbing the H-field allows to compensate manufacturing tolerances within a bandwidth of $\approx 10$~MHz.

Shielding is required around the detecting cavity and the microwave receiver to eliminate ambient electromagnetic interference (EMI) and to mitigate coupling to the emitting cavity by simple EM leakage. This would generate false positive results as a signal originating from leakage can not be distinguished from a signal propagating by WISP conversion. Within 15 cm, the field strength must be reduced by at least a factor of $7.7 \cdot 10^{12} = 258 \mathrm{~dB}$ to get meaningful results\footnotemark. \footnotetext{For comparison, the signals from the Voyager 1 space probe, 14 light hours from earth, are currently attenuated by a factor of 240 dB (comparing the field strength in front of sending and receiving antenna)}

The shielding box has been built from a straight piece of WR-2300 waveguide with the inside dimensions 58x29x100~cm. Feeding all RF signals over optical fibres by analog transceivers prevents the propagation of EMI trough ordinary transmission lines. An optical ethernet link is used to remote control the signal analyser.

Shielding effectiveness has been measured by comparing the electric field strength in- and outside the wall with a calibrated electric field probe and a spectrum analyzer. The shielding enclosure provides $\approx 90$~dB and each of the cavities provides an additional $\approx$~110 dB of shielding. The combined EM attenuation is $\approx 310$~dB, making thermal noise the limiting factor for the minimum detectable signal.

For an exclusion result it is necessary to prove the detector is working and actually able to pick up any WISP related signals.
Detection sensitivity will be limited if the resonant frequency of one or both cavities does not equal the system frequency $f_{\mathrm{sys}}$.
This is especially delicate as the cavities are sensitive to temperature variations; their resonant frequency is inversely proportional to the thermal expansion coefficient of their wall material. For brass, a temperature change of $\Delta T = +2$~K leads to a detuning of $\Delta f_{\mathrm{res}} = -112$~kHz and a reflection of around half the input power back towards the amplifier.

Thermal drift is more critical for the emitting cavity as it has to dissipate up to 50 W of heat by forced air cooling without any external temperature stabilization. Before data taking, the cavity was operated at full power and its resonant frequency was kept constant manually. After around 1 h, the cavity reached thermal equilibrium and no further tuning was necessary.
The reflected power $P_{\mathrm{refl}}$ indicates how far the emitting cavity is off tune, and is minimized for $f_{\mathrm{res}} = f_{\mathrm{sys}}$. The relation is given in Eq.~\ref{equ:cavReflFreq},
\begin{align}
\label{equ:cavReflFreq}
\frac{P_{\mathrm{refl}}}{P_{\mathrm{inc}}} = \left|\Gamma\right|^2 = \frac{f_{\mathrm{n}}^2}{4 + f_{\mathrm{n}}^2}
\quad
f_{\mathrm{n}} = Q_{\mathrm{em}} \left(\frac{f_{\mathrm{sys}}}{f_{\mathrm{res}}} - \frac{f_{\mathrm{res}}}{f_{\mathrm{sys}}}\right)
\end{align}
where $P_{\mathrm{inc}}$ is the constant ($\pm 1 \%$) and known incident RF power at the frequency $f_{\mathrm{sys}}$. Critical coupling has been assumed, with coupling losses and impedance mismatch not considered. $P_{\mathrm{refl}}$ is measured on a directional coupler, placed between power amplifier and cavity.

Tuning the detecting cavity once before the measurement run is sufficient, as it does not dissipate any power and changes in the laboratory's ambient temperature are small enough.
Its resonant frequency was estimated by evaluating the spectral noise power density $n_{\mathrm{o}}$ at the output of the low noise amplifier (LNA), given by Eq.~\ref{eq:thNoise}
\begin{align}
\label{eq:thNoise}
n_{\mathrm{o}} = G k_B \left[
T_{\mathrm{cav}} \left(1 - \left|\Gamma\right|^2\right) +
T_{\mathrm{LNA}} +
T_{\mathrm{LNA}}' \left|\Gamma\right|^2
\right]
\end{align}

where $G = 44.7$~dB is the LNA's gain and $k_B$ is the boltzmann constant. Three noise temperature terms have to be considered, $T_{\mathrm{cav}} \left(1 - \left|\Gamma\right|^2\right)$ describes the noise from the detecting cavity itself. Its noise temperature is frequency dependent by the reflection coefficient $\Gamma$. The maximum occurs at $f_{\mathrm{res}}$ and is equal to ambient temperature. $T_{\mathrm{LNA}} = 32.4$~K describes the intrinsic noise temperature of the amplifier. $T_{\mathrm{LNA}}' \left|\Gamma\right|^2$ describes the noise temperature of the amplifier input, which emits a noise wave towards the cavity where it is reflected back.
As the transmission line between amplifier and cavity is short ($\approx$ 1 cm) and as the noise temperature of the cavity walls is significantly higher than the noise wave transmitted from the amplifier's input, a good estimate of the resonant frequency can be determined from the maximum of the noise power spectrum. This has been done before and after the actual experimental run, the results are shown in Fig.~\ref{fig:diagTune} and indicate no significant (bigger than $BW_{3 \mathrm{dB}}$) drift of the emitting cavity's resonant frequency.
\begin{figure}[b]
\vspace{-0.5cm}
\begin{center}
	\includegraphics[width=0.5\textwidth]{tune.png}
	\caption{Measured spectral noise power density before and after the 11.5 h measurement run, indicating $f_{\mathrm{res}}$}
	\label{fig:diagTune}
\end{center}
\vspace{-0.5cm}
\end{figure}

The output signal of the detecting cavity is coupled out, amplified by the LNA ($G = 44.7$~dB) and then further processed by an Agilent EXA N9010A signal analyzer. The center frequency was set to $f_{\mathrm{sys}} + 4~\mathrm{Hz}$ to avoid internal spurious signals appearing at interesting parts of the spectrum. The center frequency is shifted to baseband and the complex IQ signal is digitized with 20 Hz bandwidth. The 11.5 h long time record is stored for further processing.
As the capture memory of the analyzer limits the maximum number of continuously acquired samples to $\approx 10^6$, a trade-off between recorded bandwidth and recording length had to be found. For this reason the signal analyzer will be replaced by a dedicated processing chain without this limitation for the future.

For offline data processing the spectral power of the recorded noise like signals is estimated by Welch's method implemented in a python script. The time record is read and divided into segments overlapping by $\approx 90 \%$. Each segment is multiplied by a Hann window and the complex spectra are calculated by a fast Fourier transform. Averaging these spectra trades resolution bandwidth for less variance in the noise floor. For a 11.5 h long time trace, resolution bandwidths as narrow as $BW_{\mathrm{res}} = 24~\mu$Hz can be achieved. The average noise floor is determined by $Pn = BW_{\mathrm{res}} \cdot n_o$.

The frequency error of RF-source and signal analyzer needs to be within $BW_{\mathrm{res}}$ during the experimental run, otherwise the signal power will spread out over several bins in the spectrum, degrading the signal to noise ratio. While absolute frequency drifts are unavoidable, phase-locking RF-source and signal analyzer to a common 10 MHz frequency reference allows to achieve a good relative frequency stability. This has been explained and successfully demonstrated down to $BW_{\mathrm{res}} = 10~\mu$Hz in \cite{src:narrowband}.

Immediately before the actual measurement, a test run was conducted where the shielding box was left open to provoke EM leakage between the cavities.
The resulting power spectrum contained a single peak, clearly above the noise floor, spanning only one single bin. Its absolute position on the frequency axis was offset by $\approx 3$~mHz due to the finite resolution of the RF source. As no parameters were changed after the test run and only the shielding box was closed, the WISP related signal is expected to appear at the same position in the spectrum.
\begin{figure}[t]
\begin{center}
	\includegraphics[width=0.48\textwidth]{spectra.png}
	\caption{Spectral noise power from the detecting cavity. The left diagram shows the min. and max. peaks of the recorded span, while the right diagram is zoomed}
	\label{fig:spectra}
\end{center}
\vspace{-0.75 cm}
\end{figure}
The peaks within a window of $\pm 1.5$~mHz around the expected signal do not exceed the peaks in other parts of the spectrum, as shown in Fig.~\ref{fig:spectra}. Therefore an exclusion result is produced by setting the minimum detectable power ($P_{\mathrm{det}} = -211$~dBm) to the maximum peak within the frequency window.

\section{Conclusion}
No paraphotons were observed in the first measurement-run of the microwave WISP search at CERN, improving the existing exclusion limits. Several technical challenges, like $> 300$~dB EM shielding between the cavities, keeping them on tune during the 11.5 h measurement run and filtering the signal with a bandwidth of $BW_{\mathrm{res}} = 24~\mu$Hz, had to be overcome.
\\
{\small
We are grateful for the practical hints and assistance from M.~Gasior and M.~Thumm.
Thanks to R.~Jones, E.~Jensen and the BE department management for support.
}