This IPython notebook accompanies the paper "Oscillating asymmetries in waveforms from precessing black-hole binaries". There are two main purposes to this notebook:
The notebook reflects the structure of the paper itself, with sections here corresponding to sections in the paper. Essentially, this is a powerful, extended, and interactive set of footnotes. Of course, the notebook can also be read without execution, meaning that there are no software dependencies.
This code in this notebook depends on five software modules, each freely available from github:
The last one includes the other four as submodules; to install them all, just install GWFrames
. This package also includes a lightly modified version of the spinsfast package, which implements exact transformations for spin-weighted spherical harmonics.
The test data set is assumed to come from the SXS catalog, at least some of which must be locally available to run this notebook -- though any data in the same format can be used. In particular, my favorite system is catalog number SXS:BBH:0003
, which will be used for all the examples in this notebook (as in the paper). In particular, the rhOverM_Asymptotic_GeometricUnits.h5
, Horizons.h5
, and metadata.txt
files from that system are used. Other systems could be substituted with minimal changes to the code.
We start off with some basic imports and definitions.
import sys
import os.path
import re
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import h5py
import sympy
## These are optional configurations
#mpl.rcParams['figure.figsize'] = 8,6 # Set default figure size (in inches)
#%config InlineBackend.figure_format='svg' # If plotting inline, use SVG or PNG
sympy.init_printing()
## These import some of the custom code; they must be installed properly first
import GWFrames
import GWFrames.plot
import Quaternions
## These are a few convenient symbols we use below
v, nu, delta, Sigma_n, Sigma_lambda = sympy.symbols('v, nu, delta, Sigma_n, Sigma_lambda', real=True)
Read in the numerical data:
## Set up the paths correctly
PathToSXSCatalogLinks = os.path.expanduser('~')+'/Research/Data/SimulationAnnex/CatalogLinks/'
PathToSXSCatalog = os.path.expanduser('~')+'/Research/Data/SimulationAnnex/Catalog/'
MyFavoriteSystem = PathToSXSCatalog + 'q1.0/SBBH/d19.0_q1.0_s0.5_0_0_s0_0_0/Lev5/'
MyFavoriteWaveform = 'rhOverM_Asymptotic_GeometricUnits.h5/Extrapolated_N2.dir'
MyFavoriteData = MyFavoriteSystem + MyFavoriteWaveform
## Read in the waveform and adjust it
W_NR = GWFrames.ReadFromNRAR(MyFavoriteData);
t0 = -W_NR.MaxNormTime()
t1 = 200
W_NR.SetT(W_NR.T()+t0)
print('The initial time has been moved back to {0}'.format(t0))
## Get the coordinate-position data from the h5 file
with h5py.File(MyFavoriteSystem+'Horizons.h5', 'r') as f:
tA = f['AhA.dir/CoordCenterInertial.dat'][:,0]+t0
xA = f['AhA.dir/CoordCenterInertial.dat'][:,1:]
mA = f['AhA.dir/ArealMass.dat'][:,1]
chiA = f['AhA.dir/chiInertial.dat'][:,1:]
tB = f['AhB.dir/CoordCenterInertial.dat'][:,0]+t0
xB = f['AhB.dir/CoordCenterInertial.dat'][:,1:]
mB = f['AhB.dir/ArealMass.dat'][:,1]
chiB = f['AhB.dir/chiInertial.dat'][:,1:]
## Do some calculations with the coordinates
d = xA-xB
nHat = array([[q[1],q[2],q[3]] for q in Quaternions.normalized(Quaternions.QuaternionArray(d))])
dnHatdt = array([[q[1],q[2],q[3]] for q in Quaternions.QuaternionDerivative(Quaternions.QuaternionArray(nHat), tA)])
lambdaHat = array([[q[1],q[2],q[3]] for q in Quaternions.normalized(Quaternions.QuaternionArray(dnHatdt))])
Omega_orb = cross(nHat, dnHatdt)
v = array([(o[0]**2+o[1]**2+o[2]**2)**(1./6.) for o in Omega_orb])
COM = ((xA.transpose()*mA+xB.transpose()*mB)/(mA+mB)).transpose()
## Set the average initial velocity from the inital-data ADM quantities
initial_ADM_energy = 0.0
initial_ADM_linear_momentum = []
number_regex = r'[-+eE\.\d]*'
with open(MyFavoriteSystem+'metadata.txt','r') as f:
for line in f:
if line.startswith('initial-ADM-energy'):
initial_ADM_energy = float(re.search('.* = *({0})'.format(number_regex),line).group(1))
if line.startswith('initial-ADM-linear-momentum'):
initial_ADM_linear_momentum = array([float(x)
for x in re.search('.* = *({0}), *({0}), *({0})'.format(number_regex),line).groups()])
if (len(initial_ADM_linear_momentum)!=0 and initial_ADM_energy!=0.0):
break
v_avg = initial_ADM_linear_momentum/initial_ADM_energy
print("The speed of the drift we are removing is {0}".format(sqrt(v_avg.dot(v_avg))))
## Translate the waveform to remove average velocity drift
W_NR = W_NR.Translate([v_avg*(t-W_NR.T(0)) for t in W_NR.T()])
## Make a copy of the waveform, and transform it to the co-rotating frame
UseEllModes = [2,3,4] # Which ell modes to use to calculate the angular velocity and hence corotating frame
t_alignment = t0+2000. # Time at which to align the modes to the frame
W_NR_corot = GWFrames.Waveform(W_NR)
W_NR_corot.TransformToCorotatingFrame(UseEllModes)
W_NR_corot.AlignDecompositionFrameToModes(t_alignment, Quaternions.Quaternion(nHat[argmin(abs(tA-t_alignment))]), UseEllModes);
## Just briefly read in the non-precessing waveforms
W_NR_nonpreces = GWFrames.ReadFromNRAR(PathToSXSCatalogLinks + 'SXS:BBH:0002/Lev6/' + MyFavoriteWaveform);
W_NR_nonpreces.SetT(W_NR_nonpreces.T()-W_NR_nonpreces.MaxNormTime());
W_NR_nonpreces_corot = GWFrames.Waveform(W_NR_nonpreces);
W_NR_nonpreces_corot.TransformToCorotatingFrame(UseEllModes);
W_NR_nonpreces_corot.AlignDecompositionFrameToModes(t_alignment, Quaternions.Quaternion(nHat[argmin(abs(tA-t_alignment))]), UseEllModes);
W_NR_eccentric = GWFrames.ReadFromNRAR(PathToSXSCatalogLinks + 'SXS:BBH:0091/Lev5/' + MyFavoriteWaveform);
W_NR_eccentric.SetT(W_NR_eccentric.T()-W_NR_eccentric.MaxNormTime());
W_NR_eccentric_corot = GWFrames.Waveform(W_NR_eccentric);
W_NR_eccentric_corot.TransformToCorotatingFrame(UseEllModes);
W_NR_eccentric_corot.AlignDecompositionFrameToModes(t_alignment, Quaternions.Quaternion(nHat[argmin(abs(tA-t_alignment))]), UseEllModes);
The initial time has been moved back to -11789.7489449 The speed of the drift we are removing is 6.07395108534e-05
The PN-expansion parameter $v$ is essentially a velocity of the black holes, and is defined as \begin{equation} v = \left( \frac{G\, M\, \Omega} {c^3} \right)^{1/3}. \end{equation} $\Omega$ is the orbital frequency, which (with many caveats) we can approximate as $\Omega \approx \omega/2$, where $\omega$ is the angular frequency of the emitted waves. Now, since we will be referring to the frequency $f$ observed by a gravitational-wave detector, we use $\Omega = \omega/2 = \pi f$. Thus, $v$ is roughly \begin{equation} v \approx \left( \frac{G\, M\, \pi\, f} {c^3} \right)^{1/3}. \end{equation} The lowest-mass black hole currently known has mass $3.8M_\odot$; this suggests we have no reason to expect to find any BBH with total mass greater than $7.6M_\odot$. The seismic wall in advanced GW detectors is 10Hz. Plugging these numbers in, we find $v \approx 0.106$. This is the very lowest $v$ we should expect to find for a BBH in advanced detectors.
The first task is just to look at the dominant modes themselves. We have already transformed the NR data to the corotating frame above, and aligned that frame to the modes, so in some sense, this data is as simple as it can be made with rotations.
## Plot the amplitudes of the (2,2) and (2,-2) modes
W_NR_corot.plot('Abs', [[2,2],[2,-2]])
plt.gca().set_color_cycle(None)
W_NR_nonpreces_corot.plot('Abs', [[2,2],[2,-2]], ls='dashed')
plt.gca().set_color_cycle(None)
W_NR_eccentric_corot.plot('Abs', [[2,2],[2,-2]], ls='dotted')
ylim(7e-2,0.15)
xlim(t0,t1)
title(r'Comparing the $(\ell,m) = (2,\pm 2)$ mode pair amplitudes')
legend(loc='upper left');
## Save the data to file, for direct incorporation into the main LaTeX file via `pgfplot`
numpy.savetxt('data/ModeAmplitudes_NR_2_2.dat',
array([W_NR_corot.T(),
W_NR_corot.Abs(W_NR_corot.FindModeIndex(2,2)),
W_NR_corot.Abs(W_NR_corot.FindModeIndex(2,-2))]).transpose()[::60],
fmt='%.5e')
numpy.savetxt('data/ModeAmplitudes_Nonspinning_2_2.dat',
array([W_NR_nonpreces_corot.T(),
W_NR_nonpreces_corot.Abs(W_NR_nonpreces_corot.FindModeIndex(2,2)),
W_NR_nonpreces_corot.Abs(W_NR_nonpreces_corot.FindModeIndex(2,-2))]).transpose()[::100],
fmt='%.5e')
numpy.savetxt('data/ModeAmplitudes_Eccentric_2_2.dat',
array([W_NR_eccentric_corot.T(),
W_NR_eccentric_corot.Abs(W_NR_eccentric_corot.FindModeIndex(2,2)),
W_NR_eccentric_corot.Abs(W_NR_eccentric_corot.FindModeIndex(2,-2))]).transpose()[::100],
fmt='%.5e')
## Plot the phases
W_NR_corot.plot('Arg', [[2,2],[2,-2]])
plt.gca().set_color_cycle(None)
W_NR_nonpreces_corot.plot('Arg', [[2,2],[2,-2]], ls='dashed')
plt.gca().set_color_cycle(None)
W_NR_eccentric_corot.plot('Arg', [[2,2],[2,-2]], ls='dotted')
ylim(-0.06,0.06)
xlim(t0,t1)
title(r'Comparing the $(\ell,m) = (2,\pm 2)$ mode pair phases')
legend(loc='lower left');
## Save the data to file, for direct incorporation into the main LaTeX file via `pgfplot`
numpy.savetxt('data/ModePhases_NR_2_2.dat',
array([W_NR_corot.T(),
W_NR_corot.Arg(W_NR_corot.FindModeIndex(2,2)),
W_NR_corot.Arg(W_NR_corot.FindModeIndex(2,-2))]).transpose()[::60],
fmt='%.5e')
numpy.savetxt('data/ModePhases_Nonspinning_2_2.dat',
array([W_NR_nonpreces_corot.T(),
W_NR_nonpreces_corot.Arg(W_NR_nonpreces_corot.FindModeIndex(2,2)),
W_NR_nonpreces_corot.Arg(W_NR_nonpreces_corot.FindModeIndex(2,-2))]).transpose()[::100],
fmt='%.5e')
numpy.savetxt('data/ModePhases_Eccentric_2_2.dat',
array([W_NR_eccentric_corot.T(),
W_NR_eccentric_corot.Arg(W_NR_eccentric_corot.FindModeIndex(2,2)),
W_NR_eccentric_corot.Arg(W_NR_eccentric_corot.FindModeIndex(2,-2))]).transpose()[::100],
fmt='%.5e')
As discussed in the paper, it is possible to calculate the assymetric contributions to the modes within PN theory. In particular, assuming that the black holes lie along the $x$ axis, and their velocities lie in the $x$-$y$ plane, we will see that the $h^{2,\pm2}$ modes including the overall leading-order term and the leading-order asymmetric term are given by \begin{equation} h^{2,\pm 2} = 8\nu v^{2} \sqrt{\frac{\pi}{5}}\, \frac{M}{R}\, \left\{1 \mp \frac{1}{2} v^2 \left(\Sigma_{\lambda} \pm i\Sigma_n \right) \right\}. \end{equation} We can use this expression to calculate the leading-order term in the expression for \begin{equation*} \frac{\lvert h^{2,2}\rvert - \lvert h^{2,-2} \rvert} {\lvert h^{2,2} \rvert + \lvert h^{2,-2} \rvert}: \end{equation*}
v,M = sympy.symbols('v,M', real=True)
h2p2 = 1 - v**2*(Sigma_lambda + sympy.I*Sigma_n)/(2*M**2)
h2m2 = 1 + v**2*(Sigma_lambda - sympy.I*Sigma_n)/(2*M**2)
sympy.series((sympy.sqrt(h2p2*h2p2.conjugate()) - sympy.sqrt(h2m2*h2m2.conjugate()))
/(sympy.sqrt(h2p2*h2p2.conjugate()) + sympy.sqrt(h2m2*h2m2.conjugate())), x=v, x0=0, n=3)
As mentioned above, $v$ is defined as \begin{equation} v = \left( \frac{G\, M\, \Omega} {c^3} \right)^{1/3}. \end{equation} The ISCO frequency in dimensionless units is $M\Omega = 6^{-3/2}$, corresponding to $v^2/2 = 1/12 \approx 0.08$, which gives us the scale of these oscillations at ISCO.
To get a number corresponding to the seismic wall (10Hz), we need to be more careful. $\Omega$ is the orbital frequency, which (with many caveats) we can approximate as $\Omega \approx \omega/2$, where $\omega$ is the angular frequency of the emitted waves. Now, since we will be referring to the frequency $f$ observed by a gravitational-wave detector, we use $\Omega = \omega/2 = \pi f$. We will also assume that $\vec{\Sigma}$ has its largest possible magnitude: $1$. And to understand the scale of the difference, assume $\vec{\Sigma}$ is parallel to $\hat{\lambda}$. Then, the amplitude of the relative amplitude differences, as a percentage, is given by \begin{equation} 50\,\left( \frac{G\, M\, \pi\, f} {c^3} \right)^{2/3}. \end{equation} The following function calculates this quantity:
def RelativeAmplitudeDifference(TotalMassInSolarMasses, FrequencyInHertz):
"Scale of relative amplitude differences as a percentage"
from scipy.constants import c, G
from numpy import pi
SolarMass = 1.9891e30 # kilograms
return 50*(G*TotalMassInSolarMasses*SolarMass*pi*FrequencyInHertz/c**3)**(2./3.)
And we can play around with the results a little to find the mass for which the scale at the seismic wall is 1%, and see what effect that has at the higher end of LIGO's sensitive band:
RelativeAmplitudeDifference(18.275, 10.), RelativeAmplitudeDifference(18.275, 100.), RelativeAmplitudeDifference(18.275, 1000.)
For this system, of course, the higher frequency is above ISCO, so that number is not reliable -- though it does give an idea of the scale we might expect.
We can look at the same set of frequencies for an extremely (possibly unrealistically) low-mass system:
RelativeAmplitudeDifference(10., 10.), RelativeAmplitudeDifference(10., 100.), RelativeAmplitudeDifference(10., 1000.)
The calculation for the relative amplitude differences at ISCO is independent of mass, and was done above.
Using the expression for the $\bar{A}$ operator derived in Appendix B, we can prove explicitly the assertion that the (not yet normalized) antisymmetry measure \begin{equation} \label{eq:Asymmetry} \alpha^2 = \sum_{\ell,m} \lvert{h^{\ell,m} - (-1)^{\ell+m}\bar{h}^{\ell,-m}}\rvert^{2} \end{equation} is invariant under rotations. Begin with the fact that \begin{equation} \label{eq:ModeRotation} \breve{h}^{\ell, m} = \sum_{m'} h^{\ell, m'}\, \mathfrak{D}^{(\ell)}_{m',m} \left( R^{-1} \right), \end{equation} and the property of Wigner's $\mathfrak{D}^{(\ell)}$ matrices that \begin{equation} \mathfrak{D}^{(\ell)}_{-m',-m} = (-1)^{m'+m}\, \bar{\mathfrak{D}}^{(\ell)}_{m',m}. \end{equation} Then, we can re-express the transformation law for $\bar{h}^{\ell,-m}$ obtained from Eq. $\ref{eq:ModeRotation}$ by taking the conjugate and reversing the sign of $m$: \begin{align} \breve{\bar{h}}^{\ell, -m} &= \sum_{m'} \bar{h}^{\ell, m'}\, \bar{\mathfrak{D}}^{(\ell)}_{m',-m} \left( R^{-1} \right) \\ &= \sum_{m'} \bar{h}^{\ell, -m'}\, \bar{\mathfrak{D}}^{(\ell)}_{-m',-m} \left( R^{-1} \right) \\ &= \sum_{m'} (-1)^{m'+m}\, \bar{h}^{\ell, -m'}\, \mathfrak{D}^{(\ell)}_{m',m} \left( R^{-1} \right). \end{align} Finally, using the fact that $\mathfrak{D}^{(\ell)}$ is unitary (hence its conjugate transpose is its inverse), we can calculate \begin{align} \breve{\alpha}^2 &= \sum_{\ell,m} \lvert{\breve{h}^{\ell,m} - (-1)^{\ell+m}\breve{\bar{h}}^{\ell,-m}}\rvert^{2} \\ &= \sum_{\ell,m} \left( \breve{h}^{\ell,m} - (-1)^{\ell+m}\breve{\bar{h}}^{\ell,-m} \right)\, \left( \bar{\breve{h}}^{\ell,m} - (-1)^{\ell+m}\breve{h}^{\ell,-m} \right) \\ &= \sum_{\ell,m,m',m''} \left( h^{\ell,m'} - (-1)^{\ell+m'}\bar{h}^{\ell,-m'} \right)\, \mathfrak{D}^{(\ell)}_{m',m}\, \left( \bar{h}^{\ell,m''} - (-1)^{\ell+m''}h^{\ell,-m''} \right)\, \bar{\mathfrak{D}}^{(\ell)}_{m'',m} \\ &= \sum_{\ell,m,m',m''} \left( h^{\ell,m'} - (-1)^{\ell+m'}\bar{h}^{\ell,-m'} \right)\, \mathfrak{D}^{(\ell)}_{m',m}\, \mathfrak{D}^{(\ell)\,-1}_{m,m''}\, \left( \bar{h}^{\ell,m''} - (-1)^{\ell+m''}h^{\ell,-m''} \right) \\ &= \sum_{\ell,m',m''} \left( h^{\ell,m'} - (-1)^{\ell+m'}\bar{h}^{\ell,-m'} \right)\, \delta_{m',m''}\, \left( \bar{h}^{\ell,m''} - (-1)^{\ell+m''}h^{\ell,-m''} \right) \\ &= \sum_{\ell,m} \left( h^{\ell,m} - (-1)^{\ell+m}\bar{h}^{\ell,-m} \right)\, \left( \bar{h}^{\ell,m} - (-1)^{\ell+m}h^{\ell,-m} \right) \\ &= \sum_{\ell,m} \lvert{h^{\ell,m} - (-1)^{\ell+m}\bar{h}^{\ell,-m}}\rvert^{2} \\ &= \alpha^2. \end{align} More briefly, this is \begin{equation} \breve{\alpha} = \alpha; \end{equation} the asymmetry measure $\alpha$ is rotationally invariant.
The normalization factor giving us $a$ is also invariant; the proof is just a simpler version of the above.
*Caution!* The following cell can be quite slow, because of the MinimalParityViolation
function, which runs a minimization routine at each time step. Be sure to downsample the data, or be willing to go for coffee while this runs.
%%time
W_NR_corot_short = W_NR_corot.Interpolate(W_NR_corot.T()[::10])
MinimalParityViolation_NR = W_NR_corot_short.MinimalParityViolation()
Antisymmetry_NR = W_NR_corot_short.NormalizedAntisymmetry()
CPU times: user 4min 24s, sys: 2.18 s, total: 4min 26s Wall time: 5min
# Plot the two data sets
semilogy(W_NR_corot_short.T(), MinimalParityViolation_NR, label=r'$p$')
semilogy(W_NR_corot_short.T(), Antisymmetry_NR, label=r'$a$')
ylim((2.2e-3,1))
xlim((t0,t1))
xlabel(r'$(t-t_{\mathrm{merger}})/M$')
ylabel(r'Rotationally invariant antisymmetries')
legend(loc='upper left');
# Save the data to file, for direct incorporation into the main LaTeX file via `pgfplot`
numpy.savetxt('data/Asymmetry_NR.dat', array([W_NR_corot_short.T(), Antisymmetry_NR, MinimalParityViolation_NR]).transpose(),
fmt='%.5e')
# Relative antisymmetries near merger, as percentages
i_merger = W_NR_corot_short.MaxNormIndex()
Antisymmetry_NR[i_merger:i_merger+10]*100
array([ 5.34945718, 5.79995986, 6.29248169, 6.81630673, 7.35821709, 7.90396728, 8.44139518, 8.95957282, 9.44836288, 9.90084139])
It is worth noting that we can expect these figures to be roughly four times larger for more strongly spinning systems.
This plot shows the directions of three vectors measured from the waveform
and of two vectors measured from the coordinate data for the black-hole horizons
# Just plot between these times, so it's not too messy
t_1 = W_NR.T(0)+375.
t_2 = -3000.
i_1 = argmin(abs(W_NR.T()-t_1))
i_2 = argmin(abs(W_NR.T()-t_2))
# Get the waveform vectors
AV_NRHat = array([o/np.sqrt(o.dot(o)) for o in W_NR.AngularVelocityVector()[i_1:i_2,:]])
Vh_NRHat = array([-o/np.sqrt(o.dot(o)) for o in W_NR.LLDominantEigenvector()[i_1:i_2,:]])
Ldt_NRHat = array([-o/np.sqrt(o.dot(o)) for o in W_NR.LdtVector()[i_1:i_2,:]])
# Now get the vectors for the orbital motions
i_3 = argmin(abs(tA-t_1))
i_4 = argmin(abs(tA-t_2))
Omega_orbHat = array([o/np.sqrt(o.dot(o)) for o in Omega_orb])[i_3:i_4,:]
R_orb = Quaternions.UnflipRotors(
Quaternions.FrameFromXY(
Quaternions.QuaternionArray(nHat), Quaternions.QuaternionArray(lambdaHat)
)
)
Omega_tot = Quaternions.vec(Quaternions.FrameAngularVelocity(R_orb, tA))
Omega_tot = array(Omega_tot)
Omega_totHat = array([o/np.sqrt(o.dot(o)) for o in Omega_tot])[i_3:i_4,:]
# Plot the vector directions
plot(Omega_orbHat[:,0], Omega_orbHat[:,1], label=r"$\hat{\Omega}_{\mathrm{orb}}=\hat{\ell}$")
plot(Omega_totHat[:,0], Omega_totHat[:,1], label=r"$\hat{\Omega}_{\mathrm{tot}}$")
plot(Vh_NRHat[:,0], Vh_NRHat[:,1], label=r"$\hat{V}_h$")
plot(Ldt_NRHat[:,0], Ldt_NRHat[:,1], label=r"$\langle L \partial_t \rangle$")
plot(AV_NRHat[:,0], AV_NRHat[:,1], label=r"$\hat{\omega}$")
# Plot the starting points as dots
plot([Omega_orbHat[0,0],], [Omega_orbHat[0,1],], 'o')
plot([Omega_totHat[0,0],], [Omega_totHat[0,1],], 'o')
plot([Vh_NRHat[0,0],], [Vh_NRHat[0,1],], 'o')
plot([Ldt_NRHat[0,0],], [Ldt_NRHat[0,1],], 'o')
plot([AV_NRHat[0,0],], [AV_NRHat[0,1],], 'o')
xlabel(r'$x$')
ylabel(r'$y$')
xlim((-0.1, 0.3))
ylim((-0.125,0.125))
legend();
# Save the data to file, for direct incorporation into the main LaTeX file via `pgfplot`
numpy.savetxt('data/WaveformVectors_NR.dat',
array([AV_NRHat[:,0], AV_NRHat[:,1], Vh_NRHat[:,0], Vh_NRHat[:,1],
Ldt_NRHat[:,0], Ldt_NRHat[:,1],]).transpose()[::20],
fmt='%.5e')
numpy.savetxt('data/WaveformVectors_Coordinates_NR.dat',
array([Omega_orbHat[:,0], Omega_orbHat[:,1], Omega_totHat[:,0], Omega_totHat[:,1]]).transpose()[::10],
fmt='%.5e')
The implementation of PN orbital dynamics is found in the PostNewtonian package (already installed inside of the GWFrames package). In particular, the PNTerms
subdirectory contains various IPython notebooks which collect and process the expressions. The C++
subdirectory also processes those results, and outputs code that is compiled into the GWFrames
package. This allows us to easily evolve the dynamics and construct waveforms for general PN systems, as demonstrated in Sec. III. C.
In this section, we just reproduce the plots done previously for NR, using PN data.
First, we construct the PN waveform, using as initial data the horizon quantities from NR, and align the PN to the NR waveform (adjusting its time offset and attitude).
# Get the initial conditions for the PN data
t_alignment = t0+800.
i_i = argmin(abs(tA-t_alignment))
t_i = tA[i_i]
approximant = 'TaylorT1'
delta = float( (mA[i_i]-mB[i_i])/(mA[i_i]+mB[i_i]) )
chiA_i = [float(c) for c in chiA[i_i]]
chiB_i = [float(c) for c in chiB[i_i]]
Omega_orb_i = float( np.sqrt(sum([o**2 for o in Omega_orb[i_i]])) )
R_frame_i = Quaternions.FrameFromXY([Quaternions.Quaternion(nHat[i_i]),], [Quaternions.Quaternion(lambdaHat[i_i]),])[0]
# Evolve the PN system
W_PN_coorb = GWFrames.PNWaveform(approximant, delta, chiA_i, chiB_i, Omega_orb_i, 0.9*Omega_orb_i, R_frame_i)
W_PN_coorb.SetT(W_PN_coorb.T()+t0)
# Make copies in the inertial and corotating frames
W_PN = GWFrames.PNWaveform(W_PN_coorb)
W_PN.TransformToInertialFrame();
W_PN_corot = GWFrames.PNWaveform(W_PN)
W_PN_corot.TransformToCorotatingFrame();
# Align PN to NR
GWFrames.AlignWaveforms(W_NR_corot, W_PN_corot, t_alignment+250., t_alignment+2000., 100);
W_PN = GWFrames.PNWaveform(W_PN_corot)
W_PN.TransformToInertialFrame();
The following cells simply construct the plots, as above.
# Plot the amplitudes of the (2,2) and (2,-2) modes
W_PN_corot.plot('Abs', [[2,2],[2,-2]])
plt.gca().set_color_cycle(None)
W_NR_corot.plot('Abs', [[2,2],[2,-2]], ls='dotted')
ylim((7e-2,0.15))
xlim((t0,t1))
title(r'Comparing the $(\ell,m) = (2,\pm 2)$ mode pair amplitudes')
legend(loc='upper left');
# Save the data to file, for direct incorporation into the main LaTeX file via `pgfplot`
numpy.savetxt('data/ModeAmplitudes_PN_2_2.dat', array([W_PN_corot.T(), W_PN_corot.Abs(W_PN_corot.FindModeIndex(2,2)),
W_PN_corot.Abs(W_PN_corot.FindModeIndex(2,-2))]).transpose()[::180],
fmt='%.5e')
# Plot the phases
figure()
W_PN_corot.plot('Arg', [[2,2],[2,-2]])
plt.gca().set_color_cycle(None)
W_NR_corot.plot('Arg', [[2,2],[2,-2]], ls='dotted')
ylim((-0.06,0.06))
xlim((t0,t1))
title(r'Comparing the $(\ell,m) = (2,\pm 2)$ mode pair phases')
legend(loc='lower left');
# Save the data to file, for direct incorporation into the main LaTeX file via `pgfplot`
numpy.savetxt('data/ModePhases_PN_2_2.dat', array([W_PN_corot.T(), W_PN_corot.Arg(W_PN_corot.FindModeIndex(2,2)),
W_PN_corot.Arg(W_PN_corot.FindModeIndex(2,-2))]).transpose()[::180],
fmt='%.5e')
%%time
PN_short_times = W_NR_corot_short.T()[W_NR_corot_short.T()<W_PN_corot.T(W_PN_corot.NTimes()-1)]
W_PN_corot_short = W_PN_corot.Interpolate(PN_short_times)
MinimalParityViolation_PN = W_PN_corot_short.MinimalParityViolation()
Antisymmetry_PN = W_PN_corot_short.NormalizedAntisymmetry()
CPU times: user 3min 40s, sys: 1.38 s, total: 3min 42s Wall time: 4min 5s
# Plot the two data sets
missed_minimizations = [176,]
semilogy(np.delete(W_PN_corot_short.T(),missed_minimizations), np.delete(MinimalParityViolation_PN, missed_minimizations), label=r'$p$')
semilogy(W_PN_corot_short.T(), Antisymmetry_PN, label=r'$a$')
ylim((2.2e-3,1))
xlim((t0,t1))
xlabel(r'$(t-t_{\mathrm{merger}})/M$')
ylabel(r'Rotationally invariant antisymmetries')
legend(loc='upper left');
# Save the data to file, for direct incorporation into the main LaTeX file via `pgfplot`
numpy.savetxt('data/Asymmetry_PN.dat',
np.delete(array([W_PN_corot_short.T(), Antisymmetry_PN, MinimalParityViolation_PN]),
missed_minimizations, axis=1).transpose(),
fmt='%.5e')
# Just plot between these times, so it's not too messy
W_PN_vectors = GWFrames.PNWaveform(W_PN)
W_PN_vectors.InterpolateInPlace(W_NR.T()[i_1:i_2])
i_5 = argmin(abs(W_PN_corot.T()-t_1))
i_6 = argmin(abs(W_PN_corot.T()-t_2))
# Get the waveform vectors
AV_PNHat = array([o/np.sqrt(o.dot(o)) for o in W_PN_vectors.AngularVelocityVector()])
Vh_PNHat = array([o/np.sqrt(o.dot(o)) for o in W_PN_vectors.LLDominantEigenvector()])
Ldt_PNHat = array([-o/np.sqrt(o.dot(o)) for o in W_PN_vectors.LdtVector()])
# Now get the vectors for the orbital motions
Omega_orbHat_PN = array([o/np.sqrt(o.dot(o)) for o in W_PN_coorb.Omega_orb()[i_5:i_6]])
Omega_totHat_PN = array([o/np.sqrt(o.dot(o)) for o in W_PN_coorb.Omega_tot()[i_5:i_6]])
# Plot the vector directions
plot(Omega_orbHat_PN[:,0], Omega_orbHat_PN[:,1], label=r"$\hat{\Omega}_{\mathrm{orb}}=\hat{\ell}$")
plot(Omega_totHat_PN[:,0], Omega_totHat_PN[:,1], label=r"$\hat{\Omega}_{\mathrm{tot}}$")
plot(Vh_PNHat[:,0], Vh_PNHat[:,1], label=r"$\hat{V}_h$")
plot(Ldt_PNHat[:,0], Ldt_PNHat[:,1], label=r"$\langle L \partial_t \rangle$")
plot(AV_PNHat[:,0], AV_PNHat[:,1], label=r"$\hat{\omega}$")
# Plot the starting points as dots
plot([Omega_orbHat_PN[0,0],], [Omega_orbHat_PN[0,1],], 'o')
plot([Omega_totHat_PN[0,0],], [Omega_totHat_PN[0,1],], 'o')
plot([Vh_PNHat[0,0],], [Vh_PNHat[0,1],], 'o')
plot([Ldt_PNHat[0,0],], [Ldt_PNHat[0,1],], 'o')
plot([AV_PNHat[0,0],], [AV_PNHat[0,1],], 'o')
xlabel(r'$x$')
ylabel(r'$y$')
xlim((-0.1, 0.3))
ylim((-0.125,0.125))
legend();
# Save the data to file, for direct incorporation into the main LaTeX file via `pgfplot`
numpy.savetxt('data/WaveformVectors_PN.dat',
array([AV_PNHat[:,0], AV_PNHat[:,1], Vh_PNHat[:,0], Vh_PNHat[:,1],
Ldt_PNHat[:,0], Ldt_PNHat[:,1],]).transpose()[::50],
fmt='%.5e')
numpy.savetxt('data/WaveformVectors_Coordinates_PN.dat',
array([Omega_orbHat_PN[:,0], Omega_orbHat_PN[:,1], Omega_totHat_PN[:,0], Omega_totHat_PN[:,1]]).transpose()[::50],
fmt='%.5e')
In his Appendix C, Eq. (C.7), Alcubierre gives the relation \begin{equation} {}_{s}Y_{\ell,m}(\hat{n})\, {}_{s'}Y_{\ell',m'}(\hat{n}) = \sum_{\ell'',s'',m''} {}_{s''}\bar{Y}_{\ell'',m''}(\hat{n})\, \sqrt{\frac{(2\ell+1)(2\ell'+1)(2\ell''+1)}{4\pi}} \begin{pmatrix} \ell & \ell' & \ell''\\ m & m' & m'' \end{pmatrix} \begin{pmatrix} \ell & \ell' & \ell''\\ -s & -s' & -s'' \end{pmatrix}. \end{equation} Thus, using a property of the SWSHs in the second factor, and simplifying with various properties of the $3j$ symbols, we can write \begin{align} {}_{-2}Y_{\ell,m}(\hat{n})\, {}_{-2}\bar{Y}_{\ell',m'}(\hat{n}) &= (-1)^{m'} {}_{-2}Y_{\ell,m}(\hat{n})\, {}_{2}Y_{\ell',-m'}(\hat{n}) \\ &= (-1)^{m'} \sum_{\ell'',s'',m''} {}_{s''}\bar{Y}_{\ell'',m''}(\hat{n})\, \sqrt{\frac{(2\ell+1)(2\ell'+1)(2\ell''+1)}{4\pi}} \begin{pmatrix} \ell & \ell' & \ell''\\ m & -m' & m'' \end{pmatrix} \begin{pmatrix} \ell & \ell' & \ell''\\ 2 & -2 & -s'' \end{pmatrix}. \\ &= (-1)^{m'} \sum_{\ell''} \bar{Y}_{\ell'',m'-m}(\hat{n})\, \sqrt{\frac{(2\ell+1)(2\ell'+1)(2\ell''+1)}{4\pi}} \begin{pmatrix} \ell & \ell' & \ell''\\ m & -m' & m'-m \end{pmatrix} \begin{pmatrix} \ell & \ell' & \ell''\\ 2 & -2 & 0 \end{pmatrix}. \end{align} Applied to the case of the linear momentum flux integral, where \begin{align} \hat{n}_j &= \left( \sin\vartheta\, \cos\varphi, \sin\vartheta\, \sin\varphi, \cos \vartheta \right)_j \\ &= \left( \sqrt{\frac{2 \pi}{3}} (Y_{1, -1} - Y_{1, 1}), i\sqrt{\frac{2 \pi}{3}} (Y_{1, -1} + Y_{1, 1}), \sqrt{\frac{4\pi}{3}} Y_{1,0} \right)_j, \end{align} we obtain \begin{align} \frac{dp_j}{dt} &= \frac{1}{16\pi}\, \sum_{\ell,\ell',m,m',m'''} \hat{n}_j^{1,m'''} \dot{h}^{\ell,m}\, \bar{\dot{h}}^{\ell',m'} \int Y_{1,m'''}\, {}_{-2}Y_{\ell,m}\, {}_{-2}\bar{Y}_{\ell',m'} d\Omega \\ &= \frac{1}{16\pi}\, \sum_{\ell,\ell',m,m',m'''} \hat{n}_j^{1,m'''} \dot{h}^{\ell,m}\, \bar{\dot{h}}^{\ell',m'} (-1)^{m'} \sum_{\ell''} \int Y_{1,m'''}\, \bar{Y}_{\ell'',m'-m}d\Omega\, \sqrt{\frac{(2\ell+1)(2\ell'+1)(2\ell''+1)}{4\pi}} \begin{pmatrix} \ell & \ell' & \ell''\\ m & -m' & m'-m \end{pmatrix} \begin{pmatrix} \ell & \ell' & \ell''\\ 2 & -2 & 0 \end{pmatrix} \\ &= \frac{1}{16\pi}\, \sum_{\ell,\ell',m,m'} \hat{n}_j^{1,m'-m} \dot{h}^{\ell,m}\, \bar{\dot{h}}^{\ell',m'} (-1)^{m'}\, \sqrt{\frac{3(2\ell+1)(2\ell'+1)}{4\pi}} \begin{pmatrix} \ell & \ell' & 1\\ m & -m' & m'-m \end{pmatrix} \begin{pmatrix} \ell & \ell' & 1\\ 2 & -2 & 0 \end{pmatrix}. \end{align} As a special, and particularly interesting case, we have \begin{align} \frac{dp_z}{dt} &= \frac{1}{16\pi}\, \sum_{\ell,\ell',m} \dot{h}^{\ell,m}\, \bar{\dot{h}}^{\ell',m} (-1)^{m}\, \sqrt{(2\ell+1)(2\ell'+1)} \begin{pmatrix} \ell & \ell' & 1\\ m & -m & 0 \end{pmatrix} \begin{pmatrix} \ell & \ell' & 1\\ 2 & -2 & 0 \end{pmatrix}. \end{align} It will be interesting to look at all the terms involving $\dot{h}^{2,\pm2}$, which may or may not include the lowest-order contributions: \begin{align} &\quad \sum_{\ell'} \dot{h}^{2,2}, \bar{\dot{h}}^{\ell',2} \sqrt{5(2\ell'+1)} \begin{pmatrix} 2 & \ell' & 1\\ 2 & -2 & 0 \end{pmatrix} \begin{pmatrix} 2 & \ell' & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad+\sum_{\ell'} \dot{h}^{2,-2}, \bar{\dot{h}}^{\ell',-2}, \sqrt{5(2\ell'+1)} \begin{pmatrix} 2 & \ell' & 1\\ -2 & 2 & 0 \end{pmatrix} \begin{pmatrix} 2 & \ell' & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad +\sum_{\ell} \dot{h}^{\ell,2}, \bar{\dot{h}}^{2,2}, \sqrt{(2\ell+1)5} \begin{pmatrix} \ell & 2 & 1\\ 2 & -2 & 0 \end{pmatrix} \begin{pmatrix} \ell & 2 & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad +\sum_{\ell} \dot{h}^{\ell,-2}, \bar{\dot{h}}^{2,-2}, \sqrt{(2\ell+1)5} \begin{pmatrix} \ell & 2 & 1\\ -2 & 2 & 0 \end{pmatrix} \begin{pmatrix} \ell & 2 & 1\\ 2 & -2 & 0 \end{pmatrix}\ % % &= \lvert \dot{h}^{2,2}\rvert^2, 5 \begin{pmatrix} 2 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix} \begin{pmatrix} 2 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad +\dot{h}^{2,-2}, \bar{\dot{h}}^{2,-2}, 5 \begin{pmatrix} 2 & 2 & 1\\ -2 & 2 & 0 \end{pmatrix} \begin{pmatrix} 2 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad +\dot{h}^{2,2}, \bar{\dot{h}}^{2,2}, 5 \begin{pmatrix} 2 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix} \begin{pmatrix} 2 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad +\dot{h}^{2,-2}, \bar{\dot{h}}^{2,-2}, 5 \begin{pmatrix} 2 & 2 & 1\\ -2 & 2 & 0 \end{pmatrix} \begin{pmatrix} 2 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix}\ % ell'=3 &\quad +\dot{h}^{2,2}, \bar{\dot{h}}^{3,2} \sqrt{35} \begin{pmatrix} 2 & 3 & 1\\ 2 & -2 & 0 \end{pmatrix} \begin{pmatrix} 2 & 3 & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad +\dot{h}^{2,-2}, \bar{\dot{h}}^{3,-2}, \sqrt{35} \begin{pmatrix} 2 & 3 & 1\\ -2 & 2 & 0 \end{pmatrix} \begin{pmatrix} 2 & 3 & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad +\dot{h}^{3,2}, \bar{\dot{h}}^{2,2}, \sqrt{35} \begin{pmatrix} 3 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix} \begin{pmatrix} 3 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix}\ &\quad +\dot{h}^{3,-2}, \bar{\dot{h}}^{2,-2}, \sqrt{35} \begin{pmatrix} 3 & 2 & 1\\ -2 & 2 & 0 \end{pmatrix} \begin{pmatrix} 3 & 2 & 1\\ 2 & -2 & 0 \end{pmatrix} \ % % &= \frac{4}{3} \left(\lvert\dot{h}^{2,2}\rvert^2 - \lvert \dot{h}^{2,-2}\rvert^{2} \right)
\end{align} \begin{equation} \frac{dp_z}{dt} = \frac{1}{16\pi}, \frac{4}{3} \left(\lvert\dot{h}^{2,2}\rvert^2 - \lvert \dot{h}^{2,-2}\rvert^{2} \right)
\end{equation} These look a lot like the terms in Kidder's Eq. (3.29), so we might call these the two lowest-order terms. Also, I would be hard-pressed to call that first term "interference".
We can calculate the number of evaluations needed to naively evaluate a waveform (given as a mode decomposition in a rotating frame) at a single (inertial) point.
ell = sympy.Symbol('ell', integer=True)
ell_Max = sympy.Symbol(r'\ell_{\mathrm{max}}', integer=True)
# At each time step, we would need these many elements of the Wigner D matrices:
WignerDEvaluations = sympy.summation((2*ell+1)**2, (ell,2,ell_Max)).simplify()
display(WignerDEvaluations)
for lmax in range(2,8+1):
print('Including up to ell={0} would require {1} Wigner D evaluations.'.format(lmax,WignerDEvaluations.subs(ell_Max,lmax)))
print('\n')
# We would also need these many evaluations of the SWSHs, but just once:
SWSHEvaluations = sympy.summation((2*ell+1), (ell,2,ell_Max)).simplify()
display(SWSHEvaluations)
for lmax in range(2,8+1):
print('Including up to ell={0} would require {1} SWSH evaluations.'.format(lmax,SWSHEvaluations.subs(ell_Max,lmax)))
Including up to ell=2 would require 25 Wigner D evaluations. Including up to ell=3 would require 74 Wigner D evaluations. Including up to ell=4 would require 155 Wigner D evaluations. Including up to ell=5 would require 276 Wigner D evaluations. Including up to ell=6 would require 445 Wigner D evaluations. Including up to ell=7 would require 670 Wigner D evaluations. Including up to ell=8 would require 959 Wigner D evaluations.
Including up to ell=2 would require 5 SWSH evaluations. Including up to ell=3 would require 12 SWSH evaluations. Including up to ell=4 would require 21 SWSH evaluations. Including up to ell=5 would require 32 SWSH evaluations. Including up to ell=6 would require 45 SWSH evaluations. Including up to ell=7 would require 60 SWSH evaluations. Including up to ell=8 would require 77 SWSH evaluations.
Once again, we emphasize that using Euler angles is typically a very poor method. All waveform manipulations described in this paper are better done using quaternions, and the included GWFrames
module provides all the necessary code to do so. However, in case it would be useful to use this particular method with older code, we provide the necessary formulas.
We can write down the general expression for $R_{\breve{\vartheta}, \breve{\varphi}}\, e^{\breve{\gamma}\, \hat{z}/2}$. Then, assuming that we know what the result should be (in this case, because we know the product $R_{\text{f}}\, R_{\vartheta, \varphi}$), we can solve for the angles.
# Import the basics for symbolic calculations with quaternions
from sympy.galgebra.ga import *
vartheta, varphi, gamma = sympy.symbols('varthetabreve, varphibreve, gammabreve')
X = (x,y,z) = sympy.symbols('x y z')
(ex,ey,ez,grad) = MV.setup('e_x e_y e_z',metric='[1,1,1]',coords=X)
xHat = ey^ez
yHat = -ex^ez
zHat = ex^ey
def detan(e):
"This is just a utility function to simplify atan(tan(x)) to x."
from sympy import Wild, atan, tan
a = Wild('a')
return e.replace( atan(tan(a)), a )
def decos(e):
"This is just a utility function to simplify acos(cos(x)) to x."
from sympy import Wild, acos, cos
a = Wild('a')
return e.replace( acos(cos(a)), a )
Now, we define $R_{\text{tot}} := R_{\text{f}}\, R_{\vartheta, \varphi} = R_{\breve{\vartheta}, \breve{\varphi}}\, e^{\breve{\gamma}\, \hat{z}/2}$. This is the expression in terms of the latter:
R_tot = zHat.exp(varphi/2)*yHat.exp(vartheta/2)*zHat.exp(gamma/2)
We define its quaternion components:
r_0 = R_tot.scalar().simplify()
r_1 = R_tot.coef(xHat).simplify()
r_2 = (R_tot.coef(ex^ez)).simplify()
r_3 = R_tot.coef(zHat).simplify()
display(r_0, r_1, r_2, r_3)
By inspection, we can find simple formulas for the angles: \begin{align} \breve{\vartheta} &= 2 \arccos \sqrt{r_{0}^{2} + r_{3}^{2}}, \\ \breve{\varphi} &= \arctan \frac{r_{3}} {r_{0}} - \arctan \frac{-r_{1}} {r_{2}}, \\ \breve{\gamma} &= \arctan \frac{r_{3}} {r_{0}} + \arctan \frac{-r_{1}} {r_{2}}. \end{align} The following cell shows that these formulas do correctly reduce to the desired angles.
r_a = detan(sympy.atan((-r_1/r_2).simplify()))
r_b = detan(sympy.atan((r_3/r_0).simplify()))
decos(2*sympy.acos(sympy.sqrt(r_0**2 + r_3**2).simplify())), sympy.simplify(r_b+r_a), sympy.simplify(r_b-r_a)
And we can, of course, verify that this gives us the correct results for random numbers using the Quaternions
module. We choose 1,000,000 sets of random angles in the appropriate ranges, convert them to rotations, deduce the original angles, and check that the error is small:
angles = np.array([[2*np.pi*a[0]-np.pi, np.pi*a[1], 2*np.pi*a[2]-np.pi] for a in rand(1000000,3)])
rotors = np.array([Quaternions.Quaternion(*a) for a in angles])
deduce = np.array([[np.arctan2(r[3],r[0])+np.arctan2(-r[1],r[2]), # varphi
2*np.arccos(np.sqrt(r[0]**2 + r[3]**2)), # vartheta
np.arctan2(r[3],r[0])-np.arctan2(-r[1],r[2])] # gamma
for r in rotors])
errors = np.array([np.linalg.norm(a-d) for a,d in zip(angles,deduce)])
max(errors)
Ugh. Just ugh. Euler angles are slower and more prone to numerical errors; use quaternions if at all possible.
We first need an explicit expression for the SWSHs. Just to be completely sure that our conventions are consistent, we re-derive them using the expressions used in the code of the SphericalFunctions
module.
We begin with the following definition of Wigner's $\mathfrak{D}$ matrices in terms of the quaternion rotor: \begin{multline} \mathfrak{D}^{(\ell)}_{m',m}(R) = \sqrt{ \frac{(\ell+m)!\, (\ell-m)!} {(\ell+m')!\, (\ell-m')!} } \, |R_a|^{2\ell-2m} \, R_a^{m+m'} \, R_b^{m-m'} \\ \times \sum_\rho (-1)^\rho \binom{\ell+m'}{\rho}\, \binom{\ell-m'} {\ell-\rho-m}\, \frac{|R_b|^{2\rho}}{|R_a|^{2\rho}}. \end{multline} Here, $R_a = R_0 + i\, R_3$ and $R_b = R_2 + i\, R_1$ are just convenient ways of decomposing the quaternion rotor.
The SWSHs defined in terms of the $\mathfrak{D}$ matrices as \begin{equation} \label{eq:SWSHsAsQuaternions} {}_{s}Y_{\ell,m}(\vartheta, \varphi) = (-1)^s \, \sqrt{ \frac{2\ell+1} {4 \pi} } \, \mathfrak{D}^{(\ell)}_{m, -s} \left( e^{\varphi\, \hat{z}/2}\, e^{\vartheta\, \hat{y}/2}\ \right). \end{equation}
Here, the rotor is \begin{align*} R &= \left(\cos \frac{\varphi}{2} + \hat{z} \sin \frac{\varphi}{2} \right) \left(\cos \frac{\vartheta}{2} + \hat{y} \sin \frac{\vartheta}{2} \right) \\ &= \left(\cos \frac{\varphi}{2}\cos \frac{\vartheta}{2} + \hat{y} \cos \frac{\varphi}{2}\sin \frac{\vartheta}{2} \right) + \left(\sin \frac{\varphi}{2}\cos \frac{\vartheta}{2} \hat{z} - \hat{x} \sin \frac{\varphi}{2}\sin \frac{\vartheta}{2} \right), \end{align*} so its components are \begin{align*} R_a &= \left(\cos \frac{\varphi}{2} + i\, \sin \frac{\varphi}{2} \right)\, \cos \frac{\vartheta}{2}, \\ R_b &= \left(\cos \frac{\varphi}{2} - i\, \sin \frac{\varphi}{2} \right)\, \sin \frac{\vartheta}{2}. \end{align*}
We then expand Eq. $\ref{eq:SWSHsAsQuaternions}$ to find an explicit formula for the SWSHs: \begin{align*} {}_{s}Y_{\ell,m}(\vartheta, \varphi) &= (-1)^s \, \sqrt{ \frac{2\ell+1} {4 \pi} } \, \mathfrak{D}^{(\ell)}_{m, -s} \left( e^{\varphi\, \hat{z}/2}\, e^{\vartheta\, \hat{y}/2}\ \right) \\ &= (-1)^s \, \sqrt{ \frac{2\ell+1} {4 \pi} } \, \sqrt{ \frac{(\ell-s)!\, (\ell+s)!} {(\ell+m)!\, (\ell-m)!} } \, \cos^{2\ell+2s} \frac{\vartheta}{2} \\&\qquad \times e^{i(m-s)\varphi/2} \, \cos^{m-s}\frac{\vartheta}{2} \, e^{i(m+s)\varphi/2} \, \sin^{-s-m}\frac{\vartheta}{2} \\&\qquad \times \sum_\rho (-1)^\rho \binom{\ell+m}{\rho}\, \binom{\ell-m} {\ell-\rho+s}\, \tan^{2\rho} \frac{\vartheta}{2} \\ &= (-1)^s \, \sqrt{ \frac{2\ell+1} {4 \pi} } \, \sqrt{ \frac{(\ell-s)!\, (\ell+s)!} {(\ell+m)!\, (\ell-m)!} } \\&\qquad \times \cos^{2\ell+s+m}\frac{\vartheta}{2} \, \sin^{-s-m}\frac{\vartheta}{2} \, e^{i m \varphi} \\&\qquad \times \sum_\rho (-1)^\rho \binom{\ell+m}{\rho}\, \binom{\ell-m} {\ell-\rho+s}\, \tan^{2\rho} \frac{\vartheta}{2} \\ &= (-1)^s \, \sqrt{ \frac{2\ell+1} {4 \pi} } \, \sqrt{ \frac{(\ell-s)!\, (\ell+s)!} {(\ell+m)!\, (\ell-m)!} }\, e^{i m \varphi} \\&\qquad \times \sum_\rho (-1)^\rho \binom{\ell+m}{\rho}\, \binom{\ell-m} {\ell-\rho+s}\, \cos^{2\ell-2\rho+s+m}\frac{\vartheta}{2} \, \sin^{2\rho-s-m}\frac{\vartheta}{2} \end{align*}
We have the following coordinate forms for $\vec{r}$ and its various involutions: \begin{align*} \vec{r} &= (\sin\vartheta\, \cos\varphi, \sin\vartheta\, \sin\varphi, \cos\vartheta), \\ P_x \vec{r} &= (-\sin\vartheta\, \cos\varphi, \sin\vartheta\, \sin\varphi, \cos\vartheta), \\ P_y \vec{r} &= (\sin\vartheta\, \cos\varphi, -\sin\vartheta\, \sin\varphi, \cos\vartheta), \\ P_z \vec{r} &= (\sin\vartheta\, \cos\varphi, \sin\vartheta\, \sin\varphi, -\cos\vartheta), \\ A \vec{r} &= (-\sin\vartheta\, \cos\varphi, -\sin\vartheta\, \sin\varphi, -\cos\vartheta). \end{align*} These imply the following actions of those involutions on the coordinates: \begin{align*} P_x&: (\vartheta, \varphi) \mapsto (\vartheta, \pi-\varphi), \\ P_y&: (\vartheta, \varphi) \mapsto (\vartheta, -\varphi), \\ P_z&: (\vartheta, \varphi) \mapsto (\pi-\vartheta, \varphi), \\ A&: (\vartheta, \varphi) \mapsto (\pi-\vartheta, \pi+\varphi). \end{align*}
The effects of $P_x$ and $P_y$ are fairly simple, because they only affect the $e^{im\varphi}$ term. In particular, negating $\varphi$ is essentially the same as taking the complex conjugate; adding $\pi$ just multiplies the whole thing by $(-1)^m$.
The $\vartheta \mapsto \pi-\vartheta$ transformations, however, are somewhat more complicated because they effectively swap $\cos$ and $\sin$. But we can keep the same functional form if we simultaneously take $m \mapsto -m$ (and take the complex conjugate so that the exponential is unaffected) and $\rho \mapsto \ell - \rho + s$. The latter transformation primarily induces a change to the summation variable, except that the factor $(1)^\rho$ also results in an additional sign change.
Putting the above together, we obtain \begin{align*} {}_sY_{\ell,m}(P_x \vec{r}) &= (-1)^m {}_s\bar{Y}_{\ell,m}(\vec{r}), \\ {}_sY_{\ell,m}(P_y \vec{r}) &= {}_s\bar{Y}_{\ell,m}(\vec{r}), \\ {}_sY_{\ell,m}(P_z \vec{r}) &= (-1)^{\ell+s} {}_s\bar{Y}_{\ell,-m}(\vec{r}), \\ {}_sY_{\ell,m}(A \vec{r}) &= (-1)^{\ell+s+m} {}_s\bar{Y}_{\ell,-m}(\vec{r}). \end{align*}
Given values of the field $h(\hat{r})$ in any direction $\hat{r}$, and the corresponding modes $h^{\ell,m}$, we wish to derive an expression for $h(-\hat{r})$. From the above considerations, we see that \begin{equation} {}_{s}Y_{\ell,m} (\pi-\vartheta, \pi+\varphi) = (-1)^{\ell+s+m} {}_{s}\bar{Y}_{\ell,-m} (\vartheta, \varphi). \end{equation} Then, we can calculate \begin{align} h(-\hat{r}) &= h(\pi-\vartheta, \pi+\varphi) \\ &= \sum_{\ell,m} h^{\ell,m} {}_{-2}Y_{\ell,m} (\pi-\vartheta, \pi+\varphi) \\ &= \sum_{\ell,m} (-1)^{\ell+m} h^{\ell,m} {}_{-2}\bar{Y}_{\ell,-m} (\vartheta, \varphi) \\ &= \sum_{\ell,m} (-1)^{\ell+m} h^{\ell,-m} {}_{-2}\bar{Y}_{\ell,m} (\vartheta, \varphi). \end{align} Thus, we can also write \begin{equation} \bar{h}(-\hat{r}) = \sum_{\ell,m} (-1)^{\ell+m} \bar{h}^{\ell,-m} {}_{-2}Y_{\ell,m} (\vartheta, \varphi). \end{equation} If we express conjugation and evaluation at the antipodes as the operator $\bar{A}$, we can therefore express the new modes in terms of the old as \begin{equation} \bar{A}\{h\}^{\ell,m} = (-1)^{\ell+m}\, \bar{h}^{\ell,-m}. \end{equation} Note that conjugation is needed so that both fields are expressed as an expansion in SWSHs, rather than conjugated SWSHs.