Higher order analysis of the gravitational wave velocity memory effect between two free-falling gyroscopes in the plane wave spacetime

The theory of exact plane gravitational waves differs fundamentally from that of linear gravitational waves. For exact plane waves, the derivation relies on identifying geometries that share the same symmetries as plane waves in flat spacetime, independent of Einstein’s equations [12]. The resulting metric, commonly expressed in Brinkmann coordinates $(U,V,X^{2},X^{3})$, is as follows:

$$g=2dUdV+\delta_{ij}dX^{i}dX^{j}+DdU^{2}\,,$$

(1)

where $D$ is a scalar function of coordinates ($U,X^{1},X^{2}$) with the form of

$$D=K_{ij}X^{i}X^{j}=\frac{1}{2}A_{+}(U)\bigg{[}(X^{2})^{2}-(X^{3})^{2}\bigg{]}+A_{\times}(U)X^{2}X^{3}\,.$$

(2)

Here $K_{ij}$ is a traceless, symmetric $2\times 2$ matrix that describes exact plane waves in vacuum general relativity. When the metric describe a solution to Einstein equations, it reduces to $tr(K)=0$.

The line element in the local Lorentz frame $(\hat{t},\hat{x},\hat{y},\hat{z})$ is given by

$$ds^{2}=\left(-1+\frac{\Phi}{2}\right)d\hat{t}^{2}+d\bar{x}^{2}+d\hat{y}^{2}+\left(1+\frac{\Phi}{2}\right)d\hat{z}^{2}+\Phi d\hat{t}d\hat{z}\,,$$

(3)

where we have introduced the gravitational wave potential $\Phi(\hat{t}+\hat{z},\hat{x},\hat{y})\equiv\frac{1}{2}\ddot{h}_{+}(\hat{x}^{2}-\hat{y}^{2})+\ddot{h}_{\times}\hat{x}\hat{y}$. The transformation between the local Lorentz coordinates and the Brinkmann coordinates is given by[16]

	$\displaystyle U$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{2}}\left(\hat{t}+\hat{z}\right)\,,$
	$\displaystyle V$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{2}}\left(-\hat{t}+\hat{z}\right)\,,$
	$\displaystyle X^{2}$	$\displaystyle=$	$\displaystyle\hat{x}\,,\quad X^{3}=\hat{y}\,,$

and

$\displaystyle A_{+,{\times}}(U)$

$\displaystyle=$

$\displaystyle\ddot{h}_{+,{\times}}\left(\hat{t}+\hat{z}\right)\,.$

(4)

The locally Lorentz frame can approximate the gravitational wave effects described by Brinkmann coordinates. Our analysis proceeds in local Lorentz coordinates, in line with Refs.[16, 1]. The transformation between these two coordinate systems, as defined in Eq.(4), is crucial for the calculations.

In this section, we review the more concise general solution to the geodesic equation provided by [1], building upon the results of [17]. We consider an geodesic with tangent vector $\bm{u}$,the geodesic equation for coordinate $U$ is

$\displaystyle\frac{\mathrm{d^{2}}U}{\mathrm{d}\tau^{2}}=0\,,$

(5)

Given the initial position $U(\tau_{0})=U_{0}$,the solution to Eq.(5) is given by:

$\displaystyle U=\gamma(\tau-\tau_{0})+U_{0}\,.$

(6)

The observer proper velocity reads $\bm{u}=\frac{\mathrm{d}x^{\mu}}{\mathrm{d}\tau}\partial_{\mu}=\gamma(\partial_{U}+v^{1}\partial_{V}+v^{a}\partial_{X^{a}})$ in Brinkmann cooedinates,where the ’time dilation factor’ $\gamma$ defined as

$\displaystyle\gamma=\frac{\mathrm{d}U}{\mathrm{d}\tau}=\frac{1}{\sqrt{-2v^{1}-D-v^{a}v_{a}}}\,,$

(7)

normalizes $\bm{u}$ with respect to the plane wave metric (1), i.e. $g_{\mu\nu}u^{\mu}u^{\nu}=-1$. Since there is a linear relationship between $U$ and $\tau$,we will use the coordinate $U$ instead of the geodesic affine parameter $\tau$ in the rest paper. Then, the geodesic equation of the other three coordinates is as follows [1]:

$\displaystyle\bm{\ddot{X}}=\bm{K}(U)\bm{X}\,,$

(8)

$\displaystyle\ddot{V}=-\frac{1}{2}\bm{\dot{D}}-2\bm{\dot{X}}^{T}\bm{KX}\,.$

(9)

Hereafter the dot denotes the derivative with respect to $U$. With the initial position $x^{\mu}(U_{0})=(U_{0},V_{0},\bm{X_{0}})$ and the initial four-velocity $u^{\mu}(U_{0})=\gamma(1,\dot{V_{0}},\bm{\dot{X}_{0}})$, the general solutions to Eqs. (8) and (9) are given by:

$\displaystyle\bm{X}=\bm{P}(U)\bm{X}_{0}+(U-U_{0})\bm{H(U)}\bm{\dot{X}}_{0}\,,$

(10)

$\displaystyle V=V_{0}-\frac{1}{2}\bigg{[}\bm{X}^{T}\bm{\dot{X}}-\bm{X}_{0}^{T}\bm{\dot{X}}_{0}+\frac{1}{\gamma^{2}}(U-U_{0})\bigg{]}\,,$

(11)

and

$\displaystyle\bm{\dot{X}}=\bm{\dot{P}}\bm{X}_{0}+\bigg{[}(U-U_{0})\bm{\dot{H}}(U)+\bm{H}(U)\bigg{]}\bm{\dot{X}}_{0}\,.$

(12)

Both $\bm{P}$ and $\bm{H}$ are $2\times 2$ matrices that satisfy the following equations, respectively:

$\displaystyle\bm{\ddot{P}}=\bm{K}\bm{P}\,,$

(13)

$\displaystyle(U-U_{0})\bm{\ddot{H}}+2\bm{\dot{H}}=(U-U_{0})\bm{K}\bm{H}\,.$

(14)

At this point, the initial conditions for Eqs.(13) and (14) are given by $\bm{P}(U_{0})=\bm{I}$, $\bm{H}(U_{0})=\bm{I}$, and $\bm{\dot{P}}(U_{0})=0$, $\bm{\dot{H}}(U_{0})=0$, respectively. If we define $(U-U_{0})\bm{H}=\bm{M}$, then Eq. (14) becomes $\bm{\ddot{M}}=\bm{K}\bm{M}$, with initial conditions $\bm{M}(U_{0})=0$ and $\bm{\dot{M}}(U_{0})=\bm{I}$. It is evident that obtaining the nonlinear expansion terms for $\bm{P}$ and $\bm{M}$ depends on solving Eqs.(13) and (14), for which the initial conditions are known. These equations can be solved iteratively using the Dyson series to derive the general forms of $\bm{P}(U)$ and $\bm{M}(U)$, as discussed in Section III. In [18], Sturm-Liouville theory is identified as key to the memory effect. However, for plane waves, it may be more direct to assert that the Dyson series plays a central role in the memory effect.

In [1], a method for observing the velocity memory effect statically in a plane wave spacetime is presented by comparing the angular deviation between two freely falling gyroscopes after passing through a gravitational wave region. One gyroscope, $G_{0}$, is placed at the origin, while the other, $G_{1}$, is positioned along the $X^{2}$ axis. Initially, the gyroscopes are aligned, and after passing through the gravitational wave region, the angular deviation between the two is measured. It is shown that the precession angle is independent of the specific path, allowing the gyroscopes’ worldlines to overlap for a direct comparison of their angular deviation. A freely falling gyroscope has a spin vector $S$ that obeys Fermi-Walker transport:

$$(\bm{u}\cdot\nabla\bm{S})^{\mu}=(u^{\mu}a^{\nu}-u^{\nu}a^{\mu})S_{\nu},$$

where $\bm{e}_{\hat{\mu}}$ is a local tetrad at the gyroscope’s location, with $\bm{e}_{\hat{0}}=\bm{u}$ and $\bm{e}_{\hat{\mu}}\cdot\bm{e}_{\hat{\nu}}=\eta_{\hat{\mu}\hat{\nu}}$. This allows the parallel transport to be expressed as the precession equation.

$\displaystyle\frac{\mathrm{d}S^{\hat{i}}}{\mathrm{d}\tau}=\Omega^{\hat{i}}_{\hat{j}}S^{\hat{j}}=-u^{a}\omega_{a}^{\hat{i}\hat{j}}S^{\hat{j}}\,,$

(15)

when $\omega_{a}^{\hat{i}{j}}$ is the spin connection one-form associated with the tetrad ${\bm{e}_{\hat{\mu}}}$,and $\Omega$ canbe considered as a 2-form of angular velocity. The resulting angular deviation between the two gyroscopes is thus given by [1] :

$\displaystyle\Theta^{a}=-\arcsin\left(\frac{v^{a}}{\sqrt{2}}\right)\,.$

(16)

The angular deviation between two gyroscopes depends on $v$, and thus, it is essentially a manifestation of the velocity memory effect. The angular deviation is calculated under initial conditions of a separation displacement $L$ along the $X^{2}$-direction and a separation velocity $v_{0}$. To facilitate the observation of this angular deviation, numerical simulations are performed within the gravitational collapse model, where both $\bm{P}$ and $\bm{M}$ are expanded to linear order. The TT-gauge metric is a trivial linearization of the Baldwin-Jeffery-Rosen (BJR) coordinates [12]. The linear expansion of $P$ utilizes the relationship between Brinkmann coordinates and BJR coordinates [10]. From $a_{ij}=\delta_{ij}+h_{ij}(u)+\dots$, the expression $P_{ij}=\delta_{ij}+\frac{1}{2}h_{ij}(u)+\dots$ is obtained. The BJR coordinates allow for a comparison between linear and nonlinear gravitational waves in the TT gauge. In Section III, it will be shown that, without resorting to the B coordinates, the expressions for $\bm{P}$ and $\bm{M}$ can be derived by analogy with the Dyson series.

The second-order differential equation $\bm{\ddot{P}}(U)=\bm{K}(U)\bm{P}(U)$, along with the initial conditions $\bm{P}(U_{0})=\bm{I},\bm{\dot{P}}(U_{0})=0$, allows

$$\bm{P}(U)=\bm{I}+\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\bm{K}(U_{2})\bm{P}(U_{2})$$

to be expressed. By iteratively applying $\bm{P(U)}$, the general solution is obtained,namely[12]

$$\bm{P}(U)=\bm{I}+\sum_{j=1}^{n}\left(\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\bm{K}(U_{2})\prod_{i=2}^{j}\int_{U_{0}}^{U_{2i-2}}dU_{2i-1}\int_{U_{0}}^{U_{2i-1}}dU_{2i}\bm{K}(U_{2i})\right)\,,$$

(17)

Let $(U-U_{0})\bm{H}(U)=\bm{M}(U)$ hold, then $\bm{\dot{M}}(U)=(U-U_{0})\bm{\dot{H}}(U)+\bm{H}(U)$,$\bm{\ddot{M}}(U)=\bm{K}(U)\bm{M}(U)$ follows. With the initial conditions given by $\bm{M}(U_{0})=0$,$\bm{\dot{M}}(U_{0})=\bm{I}$, the result is then obtained,

$$\bm{M}(U)=\int_{U_{0}}^{U}dU_{1}\bm{I}+\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\bm{K}(U_{2})\bm{M}(U_{2})\,.$$

By iteratively applying $\bm{M(U)}$, the final result is obtained,

$$\bm{M}(U)=\int_{U_{0}}^{U}dU_{1}\bm{I}+\sum^{n}_{j=1}\left(\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\bm{K}(U_{2})\prod^{j}_{i=2}\int_{U_{0}}^{U_{2i-2}}dU_{2i-1}\int_{U_{0}}^{U_{2i-1}}dU_{2i}\bm{K}(U_{2i})\int_{U_{0}}^{U_{2j}}dU_{2j+1}\bm{I}\right)\,.$$

(18)

From Eqs.(17) and (18) can be expanded to arbitrary orders. Note that $\tau$ is the time-ordering operator.

When only the initial conditions are considered, with a separation displacement $L$ along the $X^{2}$-direction, $\Delta V=\dot{P}L$ can be obtained from the initial conditions $X_{0}=L$, $\dot{X}_{0}=0$, and Eq.(12). Using Eq.(16), the angular deviation is determined. For clarity, its matrix form is explicitly given as follows:

$$\begin{pmatrix}\Delta\Theta^{2}\\ \Delta\Theta^{3}\end{pmatrix}=-\arcsin\left(\frac{1}{\sqrt{2}}\begin{pmatrix}\dot{P}_{22}&\dot{P}_{23}\\ \dot{P}_{32}&\dot{P}_{33}\end{pmatrix}\begin{pmatrix}L\\ 0\end{pmatrix}\right)=-\arcsin\left(\frac{1}{\sqrt{2}}\begin{pmatrix}\dot{P}_{22}L\\ \dot{P}_{32}L\end{pmatrix}\right)\,.$$

(19)

To determine the angle difference, one must first calculate ${P}_{22}$ and ${P}_{32}$, after which the matrix expansion of ${P}$ can be written as

			$\displaystyle\begin{pmatrix}P_{22}(U)&P_{23}(U)\\ P_{32}(U)&P_{33}(U)\end{pmatrix}=\begin{pmatrix}1&0\\ 0&1\end{pmatrix}+\frac{1}{2}\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\begin{pmatrix}A_{+}(U_{2})&A_{\times}(U_{2})\\ A_{\times}(U_{2})&-A_{+}(U_{2})\end{pmatrix}$		(20)
		$\displaystyle+$	$\displaystyle\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\begin{pmatrix}A_{+}(U_{2})&A_{\times}(U_{2})\\ A_{\times}(U_{2})&-A_{+}(U_{2})\end{pmatrix}\int_{U_{0}}^{U_{2}}dU_{3}\int_{U_{0}}^{U_{3}}dU_{4}\begin{pmatrix}A_{+}(U_{4})&A_{\times}(U_{4})\\ A_{\times}(U_{4})&-A_{+}(U_{4})\end{pmatrix}+\cdots\,.$

To facilitate the observation of angular variations, we use the toy model of gravitational collapse from [19], namely:

$\displaystyle A_{+}(U)=\ddot{h}_{+}=\frac{1}{2}\frac{\mathrm{d}^{3}(e^{-U^{2}})}{\mathrm{d}U^{3}}\,.$

(21)

For simplicity, we begin by considering linear polarization, where $A_{\times}=0$. Later, the case of $A_{\times}\neq 0$ will also be examined. By considering a third-order expansion and substituting (21) into (20),we obtain:

	$\displaystyle P_{22}$	$\displaystyle=1+\frac{1}{2}h_{+}(U)+\frac{1}{8}h_{+}(U)^{2}+\frac{1}{48}h_{+}(U)^{3}$		(22)
		$\displaystyle-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\dot{h}_{+}(U_{2})^{2}-\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\dot{h}_{+}(U_{2})^{2}h_{+}(U_{2})$
		$\displaystyle-\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\int_{U_{0}}^{U_{1}}dU_{2}\ddot{h}_{+}(U_{2})\int_{U_{0}}^{U_{2}}dU_{3}\int_{U_{0}}^{U_{3}}dU_{4}\dot{h}_{+}(U_{4})^{2}+\dots\,.$

So that

	$\displaystyle\Delta\Theta^{2}$	$\displaystyle=-\arcsin\left(\frac{1}{\sqrt{2}}\dot{P}_{22}L\right)$		(23)
		$\displaystyle=-\frac{L}{\sqrt{2}}\bigg{(}\frac{1}{2}\dot{h}_{+}(U)+\frac{1}{4}\dot{h}_{+}(U)h_{+}(U)+\frac{1}{16}\dot{h}_{+}(U)h_{+}(U)^{2}-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}$
		$\displaystyle-\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}h_{+}(U_{1})-\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\int_{U_{0}}^{U_{2}}dU_{3}\dot{h}_{+}(U_{3})^{2}+\dots\bigg{)}\,,$

for small angles. From Fig.1, it is clear that when considering only the initial conditions with a separation displacement $L$ along the $X^{2}$ direction, no final angular deviation arises in the linear approximation of the simplified gravitational collapse model. However, an angular deviation emerges starting from the second-order approximation, confirming that higher-order terms introduce angular deviation [1], beginning at second order. By the third-order approximation, the results closely match numerical simulations. In our calculations, this assumption is implicitly applied for simplicity, specifically regarding $h_{+}(U_{0})=0$ and $\dot{h}_{+}(U_{0})=0$.

Considering only the initial conditions, with separation velocity along the $X^{2}$ direction at a rate of $V_{0}$, the initial state is defined as $\dot{X}_{0}=V_{0}$. To provide a clear representation of the angular deviation, we begin by expressing it in matrix form:

$$\begin{pmatrix}\Delta\Theta^{2}\\ \Delta\Theta^{3}\end{pmatrix}=-\bigg{[}\arcsin(\frac{1}{\sqrt{2}}\begin{pmatrix}V^{2}\\ V^{3}\end{pmatrix})-\arcsin(\frac{1}{\sqrt{2}}\begin{pmatrix}V_{0}\\ 0\end{pmatrix})\bigg{]}\,,$$

(24)

where $V^{2}$ and $V^{3}$ must first be calculated to determine the angular deviation:

$$\begin{pmatrix}V^{2}\\ V^{3}\end{pmatrix}=\bigg{(}(U-U_{0})\begin{pmatrix}\dot{H}_{22}&\dot{H}_{23}\\ \dot{H}_{32}&\dot{H}_{33}\end{pmatrix}+\begin{pmatrix}H_{22}&H_{23}\\ H_{32}&H_{33}\end{pmatrix}\bigg{)}\begin{pmatrix}V_{0}\\ 0\end{pmatrix}=\bigg{(}(U-U_{0})\begin{pmatrix}\dot{H}_{22}\\ \dot{H}_{32}\end{pmatrix}+\begin{pmatrix}H_{22}\\ H_{32}\end{pmatrix}\bigg{)}V_{0}=\begin{pmatrix}\dot{M}_{22}\\ \dot{M}_{32}\par\end{pmatrix}V_{0}\,.$$

(25)

At the beginning, $M_{22}$ and $M_{32}$ are first calculated. For simplicity, we initially focus on the $+$ polarization, where $M_{32}=0$ and $V^{3}=0$. Expanding $M_{22}$ to third order yields $V^{2}$ as follows:

	$\displaystyle V^{2}=$	$\displaystyle\left((U-U_{0})\dot{H}_{22}+H_{22}\right)V_{0}=\dot{M}_{22}V_{0}$		(26)
	$\displaystyle=$	$\displaystyle\Bigg{(}1-\frac{1}{2}h_{+}(U)+\frac{1}{2}\dot{h}_{+}(U)(U-U_{0})+\frac{1}{4}\dot{h}_{+}(U)h_{+}(U)(U-U_{0})-\frac{1}{8}h_{+}(U)^{2}$
		$\displaystyle-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}(U_{1}-U_{0})-\frac{1}{2}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}h_{+}(U_{2})$
		$\displaystyle+\frac{1}{8}\dot{h}_{+}(U)h_{+}(U)^{2}U-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}U_{1}-\frac{1}{24}h_{+}(U)^{3}$
		$\displaystyle-\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\dot{h}_{+}(U_{2})h_{+}(U_{2})U_{2}-\frac{3}{16}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}h_{+}(U_{2})^{2}$
		$\displaystyle-\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\int_{U_{0}}^{U_{2}}dU_{3}\dot{h}_{+}(U_{3})^{2}U_{3}-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\int_{U_{0}}^{U_{2}}dU_{3}$
		$\displaystyle\ddot{h}_{+}(U_{3})\int_{U_{0}}^{U_{3}}dU_{4}h_{+}(U_{4})-\frac{1}{16}\dot{h}_{+}(U)h_{+}(U)^{2}U_{0}+\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}h_{+}(U_{1})U_{0}$
		$\displaystyle+\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\int_{U_{0}}^{U_{2}}dU_{3}\dot{h}_{+}(U_{3})^{2}U_{0}\Bigg{)}V_{0}\,.$

Thus, we have $\Delta\Theta^{2}$ :

$$\begin{aligned} \Delta\Theta^{2}&=-\bigg{[}\arcsin\left(\frac{V^{2}}{\sqrt{2}}\right)-\arcsin\left(\frac{V_{0}}{\sqrt{2}}\right)\bigg{]}\\ &=-\frac{V_{0}}{\sqrt{2}}\Bigg{(}-\frac{1}{2}h_{+}(U)+\frac{1}{2}\dot{h}_{+}(U)(U-U_{0})+\frac{1}{4}\dot{h}_{+}(U)h_{+}(U)(U-U_{0})-\frac{1}{8}h_{+}(U)^{2}\\ &\quad-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}(U_{1}-U_{0})-\frac{1}{2}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}h_{+}(U_{2})\\ &\quad+\frac{1}{8}\dot{h}_{+}(U)h_{+}(U)^{2}U-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}U_{1}-\frac{1}{24}h_{+}(U)^{3}\\ &\quad-\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\dot{h}_{+}(U_{2})h_{+}(U_{2})U_{2}-\frac{3}{16}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}h_{+}(U_{2})^{2}\\ &\quad-\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\int_{U_{0}}^{U_{2}}dU_{3}\dot{h}_{+}(U_{3})^{2}U_{3}-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\int_{U_{0}}^{U_{2}}dU_{3}\\ &\quad\ddot{h}_{+}(U_{3})\int_{U_{0}}^{U_{3}}dU_{4}h_{+}(U_{4})-\frac{1}{16}\dot{h}_{+}(U)h_{+}(U)^{2}U_{0}+\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}h_{+}(U_{1})U_{0}\\ &\quad+\frac{1}{8}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}\int_{U_{0}}^{U_{2}}dU_{3}\dot{h}_{+}(U_{3})^{2}U_{0}\Bigg{)}\,.\end{aligned}$$

(27)

As shown in Fig.2, for initial conditions with a separation velocity $V_{0}$ only in the $X^{2}$ direction, no final angular deviation is observed under the linear approximation. However, when second-order and third-order corrections are included, a final angular deviation emerges.

When the initial conditions are considered, with both displacement and velocity separated, the analysis is limited to second order. Using Eqs. (12),(16),(23) and (27),we have

	$\displaystyle\Delta\Theta^{2}$	$\displaystyle=-\frac{1}{\sqrt{2}}(\dot{P}_{22}L+\dot{M}_{22}V_{0})$		(28)
		$\displaystyle=-\frac{L}{\sqrt{2}}\left(\frac{1}{2}\dot{h}_{+}(U)+\frac{1}{4}\dot{h}_{+}(U)h_{+}(U)-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}\right)$
		$\displaystyle-\frac{V_{0}}{\sqrt{2}}\bigg{(}-\frac{1}{2}h_{+}(U)+\frac{1}{2}\dot{h}_{+}(U)(U-U_{0})+\frac{1}{4}\dot{h}_{+}(U)h_{+}(U)(U-U_{0})-\frac{1}{8}h_{+}(U)^{2}$
		$\displaystyle\quad-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}(U_{1}-U_{0})-\frac{1}{2}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}h_{+}(U_{2})\bigg{)}\,.$

The oscillation modes are inactive, and restore the standard unit, we can get

	$\displaystyle\Delta\Theta^{2}$	$\displaystyle=\frac{L}{\sqrt{2}c}\bigg{(}\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}\bigg{)}-\frac{V_{0}}{\sqrt{2}c}\bigg{(}-\frac{1}{2}h_{+}(U)-\frac{1}{8}h_{+}(U)^{2}$		(29)
		$\displaystyle\quad-\frac{1}{4}\int_{U_{0}}^{U}dU_{1}\dot{h}_{+}(U_{1})^{2}(U_{1}-U_{0})-\frac{1}{2}\int_{U_{0}}^{U}dU_{1}\ddot{h}_{+}(U_{1})\int_{U_{0}}^{U_{1}}dU_{2}h_{+}(U_{2})\bigg{)}\,.$

For the plane wave approximation to hold, the maximum initial separation distance $L$ between two gyroscopes with a separation velocity $v_{0}\ll c$ is $10^{17}$ meters [16, 1]. Considering that the strain amplitude from the merger of two supermassive black holes is about $10^{-15}$ [20], the resulting contribution angle due to the initial separation $L$ is approximately $10^{-23}$ radians. The contribution from the initial separation displacement only becomes significant at higher orders [1], specifically at second order. This is because the first-order expansion of $P=1+\frac{1}{2}h$ with respect to the initial separation $L$ corresponds to the linear gravitational waves, which do not generate a velocity memory effect. The angular formula $\theta=-\frac{L}{2\sqrt{2}}\dot{h}_{+}(U)$, derived from the initial separation distance $L$, corresponds to the first-order expansion $P=1+\frac{1}{2}h$ of $\bm{P}$, where the oscillating term does not contribute to the memory effect.