\documentclass[prd,aps,nofootinbib,showpacs,notitlepage,twocolumn]{revtex4-1} %\usepackage[letterpaper,margin=1in,footskip=0.5in]{geometry} \usepackage{amsfonts} \usepackage{amsmath} \usepackage{amssymb} \usepackage{bm} \usepackage{dcolumn} \usepackage[dvips]{graphicx} \usepackage{graphics} \usepackage[latin1]{inputenc} \usepackage{latexsym} \usepackage{rotating} \usepackage[colorlinks=true]{hyperref} \usepackage{xspace} % Sensible space treatment at end of simple macros \usepackage[usenames]{color} \usepackage{mathrsfs} \usepackage{multirow} \usepackage{pifont} \usepackage{enumitem} %\usepackage{ulem} \linespread{1.0} %% Try to control orphans, widows, and extra whitespace \widowpenalty=1000 \clubpenalty=1000 \raggedbottom \definecolor {darkgreen}{rgb}{0.2,0.7,0.2} \definecolor{purple}{rgb}{0.5,0,0.5} \newcommand{\red}{\textcolor{red}} \newcommand{\blue}{\textcolor{blue}} \newcommand{\green}{\textcolor{green}} \newcommand{\cyan}{\textcolor{cyan}} \newcommand{\ny}[1]{\textcolor{blue}{\it{\textbf{ny: #1}}} } \newcommand{\kc}[1]{\textcolor{cyan}{\textbf{kc: #1}} } \newcommand{\comm}{\textcolor{red}} \newcommand{\eq}{\begin{equation}} \newcommand{\be}{\begin{equation}} \newcommand{\eeq}{\end{equation}} \newcommand{\ee}{\end{equation}} \newcommand\T{\rule{0pt}{3ex}} \newcommand\B{\rule[-2ex]{0pt}{0pt}} \newcommand{\mnras}{Mon.\ Not.\ Roy.\ Astron.\ Soc.\ } \newcommand{\apjl}{Astrophys.\ J.\ Lett.\ } \newcommand{\aap}{Astron.\ Astrophys.\ } \newcommand{\cqg}{Class.\ Quantum Grav.\ } \newcommand{\araa}{Ann.\ Rev.\ Astron.\ Astrophys.\ } \newcommand{\pasp}{Publ. Astron. Soc. Pac.\ } \newcommand{\NZ}{{\mbox{\tiny NZ}}} \newcommand{\GW}{{\mbox{\tiny GW}}} \newcommand{\ppE}{{\mbox{\tiny ppE}}} \newcommand{\IM}{{\mbox{\tiny IM}}} \newcommand{\IMR}{{\mbox{\tiny IMR}}} \newcommand{\MR}{{\mbox{\tiny MR}}} \newcommand{\RD}{{\mbox{\tiny RD}}} \newcommand{\FZ}{{\mbox{\tiny FZ}}} \newcommand{\plusonetotwo}{+(1\leftrightarrow 2)} \newcommand{\GR}{{\mbox{\tiny GR}}} \newcommand{\dip}{{\mbox{\tiny Dip}}} \newcommand{\ED}{{\mbox{\tiny ED}}} \newcommand{\EA}{{\mbox{\tiny EA}}} \newcommand{\KG}{{\mbox{\tiny KG}}} \newcommand{\MG}{{\mbox{\tiny MG}}} \newcommand{\MDR}{{\mbox{\tiny MDR}}} \newcommand{\low}{{\mbox{\tiny low}}} \newcommand{\hi}{{\mbox{\tiny high}}} \newcommand{\hicut}{{\mbox{\tiny high-cut}}} \newcommand{\locut}{{\mbox{\tiny low-cut}}} \newcommand{\ratiolo}{{\mbox{\tiny lratio}}} \newcommand{\ratiohi}{{\mbox{\tiny hratio}}} \newcommand{\isco}{{\mbox{\tiny isco}}} \newcommand{\yrs}{{\mbox{\tiny 3 years}}} \newcommand{\spc}{{\mbox{\tiny space}}} \newcommand{\grnd}{{\mbox{\tiny ground}}} \newcommand{\fdamp}{{f_{\mbox{\tiny DAMP}}}} \newcommand{\damp}{{\mbox{\tiny DAMP}}} \newcommand{\MAT}{{\mbox{\tiny mat}}} \newcommand{\nn}{\nonumber} \newcommand{\rsep}{r_{12}} \newcommand{\vsep}{v_{12}} \newcommand{\hDef}{\mathfrak{h}} \newcommand{\LL}{{\mbox{\tiny LL}}} \newcommand{\TT}{{\mbox{\tiny TT}}} \newcommand{\sd}{\partial} \newcommand{\tilh}{\tilde{h}} \newcommand{\Ln}{\mbox{ln}} \newcommand{\scM}{\mathcal{M}} \newcommand{\scA}{\mathcal{A}} \newcommand{\mfb}{\mathfrak{b}} \newcommand{\mfa}{\mathfrak{a}} \newcommand{\mfc}{\mathfrak{c}} \newcommand{\mfd}{\mathfrak{d}} \newcommand{\mfe}{\mathfrak{e}} \newcommand{\mfg}{\mathfrak{g}} \newcommand{\frd}{f_{\RD}} \usepackage{tabularx} \usepackage{array} \usepackage{tabulary} \renewcommand{\arraystretch}{1.5} %%Defines new column type that is centered in tabularx \newcolumntype{Y}{>{\centering\arraybackslash}X} \newcolumntype{S}{>{\hsize=.2\hsize}X} \newcolumntype{M}{>{\centering\arraybackslash}S} \definecolor{mypink}{RGB}{232,132,161} \newcommand{\pink}[1]{\textcolor{mypink}{\textbf{#1}} } \usepackage{fancyhdr} \pagestyle{fancy} \fancyhf{} \rhead{Katie Chamberlain} \lhead{Black Hole Kicks} \cfoot{\thepage} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} \title{ Constraints on Black Hole Kicks \\ with Gravitational Wave Observations\\\textnormal{\footnotesize Interim Report One}\vspace{.1\baselineskip}} \begin{abstract} \emph{\textbf{Objective:}} We aim to develop a kicked gravitational waveform model in the frequency domain that can be used to place projected constraints on measurements of kick velocities with future ground- and space-based gravitational wave detectors. \end{abstract} \date{\today} \maketitle \thispagestyle{fancy} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Background} % What are black hole kicks? Black hole binaries are generically expected to contain black holes that are unequal in both mass and spin. These systems will display a high amount of asymmetry such that gravitational waves (GWs) produced during the evolution of the binary are radiated anisotropically. This anisotropic beaming of GWs causes the center of mass of the binary to move over time, changing as a function of the binary properties. During the final moments of the inspiral and merger, this anisotropy imparts some linear momentum to the black hole remnant that will then move in some direction with respect to the line of sight. This will cause the GWs emitted during late inspiral and merger-ringdown to be red- or blue-shifted according to the center of mass velocity. This velocity is called the recoil (or kick) velocity~\cite{Hannam:2013oca}. This process results in what is known as a ``black hole kick.'' % Why would we want to measure them? It was recently determined that black hole kicks could be directly detected using space-based and future ground-based GW detectors~\cite{Gerosa:2016vip}. In the case of space-based instruments, these measurements could be of particular importance to event rates for supermassive black hole binaries mergers if the remnant black holes of these mergers receive kicks that exceed the escape velocity of the host galaxy~\cite{Gerosa:2014gja}. \section{Approach} % How we are going to figure out how well we can measure black hole recoils? In order to determine how well we will be able to measure black hole kicks, we need to begin with a GW waveform that encapsulates the redshifting (blueshifting) of the GWs emitted from binaries during inspiral, merger, and ringdown. Ref.~\cite{Gerosa:2016vip} provides a model for the kick velocity as a function of time. However, we need the velocity as a function of frequency such that it can be incorporated into frequency-domain waveforms. In particular, we have an analytic phenomenological frequency-domain waveform, IMRPhenomD~\cite{Khan:2015jqa, Santamaria:2010yb,Husa:2015iqa}, that we need to modify to include black hole kicks. % What happens after we get the waveform? Once we have a kicked frequency-domain waveform, we can perform a Fisher Analysis to determine constraints on these kicks. However, this study requires derivatives of the frequency-domain waveform with respect to binary intrinsic and extrinsic parameters, as well as the kick velocity. In the code that was used for other similar analyses, these derivatives are done analytically to get rid of numerical error that is inherent in highly oscillatory numerical derivatives. Thus, it is of relative importance to obtain the velocity analytically as a function of frequency \textit{and} binary parameters. \section{Progress} % Discuss briefly how we explored the SPA. At the beginning of the project, we thought that obtaining the time evolution of the frequency of the GW in PhenomD would be relatively straightforward. It is easy to obtain $f(t)$ given a time-domain waveform by differentiating the argument with respect to time, but we do not have a waveform analytically in the time-domain. Thus we determined that it would be reasonable to complete an approximate analytic inverse Fourier transform (IFT) of our frequency-domain waveform. This can be done simply in some cases using the Stationary Phase Approximation (SPA)~\cite{Yunes:2009yz}. This approach appeared promising, but provided unphysical results in which time did not increase monotonically. We determined that this was due to a number of the assumptions that were made in order to use the SPA. In particular, we determined that after the inspiral regime of PhenomD, the amplitude of the waveform was oscillating too rapidly with respect to the phase for the SPA to be valid. % Discuss briefly how we explored the Method of steepest descents. We then explored whether or not our IFT could be completed using an asymptotic technique called the method of steepest descents~\cite{msd}. However, our integral was too complicated to be expanded in a convergent sum, so the method of steepest descent would not provide an easier way to obtain our time-domain waveform. % Discuss how we are currently attempting to do this integral. Currently, we are trying to integrate our waveform analytically using reasonable approximations to the waveform. This has proved successful for the merger-ringdown section of the PhenomD waveform\footnote{ The method by which we completed the IFT is presented in Appendix A.}, but we need to complete this same integration for four other waveform sections. We have also determined that we need to check that the inclusion of the velocity as a function of frequency into the frequency-domain PhenomD waveform gives the same results as introducing the kick into the time-domain waveform numerically and then taking the Fourier transform. If so, then it is reasonable to continue looking for a way to find the velocity as a function of frequency, but if not we may have no choice but to proceed numerically. \section{Challenges} % Discuss how built in mathematica functions are not well-equipped to handle what we are attempting. I have written a Mathematica script that not only computes the analytic PhenomD waveform, but also that uses the waveform to complete a Fisher analysis. However, in the event that we need to do the analysis numerically, it make take some time to implement the new process in Mathematica or a new language entirely. I have considered converting my Mathematica script to Python because it would complete numerical analyses much more quickly than Mathematica, however it may be useful to integrate the analytic derivatives taken by Mathematica and the numerical ones taken by Python. % Discuss how LAL's implementation of windowing and time-shifting may be challenging to replicate with our analytic waveform. We also need to consider how we will verify the kicked waveform in either the time or frequency domain. This can be done in the time domain by benchmarking our results with the LIGO Algorithm Library (LAL). However, it will be important to know how the waveform is first windowed and time-shifted in LAL in order to make accurate comparisons. Having never used (let alone ``simply'' downloaded) LAL, I will need to become familiar with it while we are testing our waveforms. \section{Future Directions} % Discuss the next step and how we are working to do it. In the future, we hope to find a relatively analytic way to implement a kick in our frequency-domain waveform. Once we do this, we will also need to redo the entire analysis for a PhenomD-type waveform model that introduces precession. Difficulties may arise if the precessing waveform differs drastically from PhenomD. We also hope to complete this analysis for a number of sources. Preferably, this analysis will be done for thousands of black hole binaries of various masses, spins, and redshifts, and such that we determine constraints for multiple detectors. Also, we have only briefly considered stacking when completing our parameter estimation study, but it seems highly nontrivial to do a stacking analysis for kicked black holes. \appendix \section{} Here, we present the method through which we are able to analytically complete the inverse Fourier transform of the merger ringdown portion of the waveform. Begin with a waveform of the form \be \tilh(f)=\scA(f)\, e^{-i\psi(f)}\ee The amplitude is \be \scA(f)\,=A\, f^{-7/6}\frac{ \, \gamma_{1}\, \gamma_{3}\, f_{\damp}\, e^{-\frac{\gamma_{2}\, (f-\frd)}{\gamma_{3}\ f_{\damp}}}}{m^{2}\,((f-f_{\RD})^{2} +(\gamma_{3}\, f_{\damp})^{2})} \ee where \be A=\bigg(\frac{\pi}{30}\bigg)^{1/2}\ \frac{\scM^{2}}{D_{L}}\ (\pi \scM)^{-7/6} \ee and $\gamma_{i}$ is a constant in frequency but is a function of system parameters. We can simplify this greatly if we reassign all constants. Let $\alpha\equiv A \gamma_{1}\gamma_{3}f_{\damp}/m^{2}$, $\beta\equiv \gamma_{3}\, f_{\damp}$, and $\mu\equiv\gamma_{2}(\gamma_{3}\, f_{\damp})^{-1}$. Then, \be \scA(f)\,=\frac{\alpha\, f^{-7/6}}{(f-f_{\RD})^{2}+\beta^{2}}\, e^{-\mu (f-\frd)}. \ee Now let's explore the phase, which is given as \be \psi(f)=\mfa+\mfb\, f -\frac{\mfc}{f}+\mfd \,f^{3/4}+\mfe \arctan\bigg(\frac{f-\mfg\, \frd}{\fdamp}\bigg) \ee where all gothic letters are constants in frequency and are found by matching coefficients to the inspiral and intermediate domains of the waveform. We are attempting to find \begin{align}\notag h(t)&=\frac{1}{2 \pi}\int \tilh(f)\, e^{2 \pi i f t}df\\&=\frac{1}{2 \pi}\int\scA(f)\, e^{-i\psi(f)}\, e^{2 \pi i f t}df\end{align} which is quite complicated for the amplitude and phase that are found above. Thus, it is necessary to make some simplifications. In the system that we have been exploring thus far (namely $m_{1}=10.25, \ m_{2}=10,\ \chi_{1}=0.001, \ \chi_{2}=0.002$), it seems that the phase is completely dominated by the term that is linear in frequency. Thus, it is straightforward to do a Taylor expansion to first order which will provide a good approximation. If we expand about $f_{0}$, then our phase becomes \begin{align}\notag \psi(f)&\simeq\bigg[\mfa-2\frac{\mfc}{f_{0}}+\frac{1}{4}\mfd \,f_{0}^{3/4}-\frac{\mfe f_{0} \fdamp}{\fdamp^{2}+(f_{0}-\mfg f_{\RD})^{2}}\\\notag &+\mfe \arctan\bigg(\frac{f_{0}-\mfg\, \frd}{\fdamp}\bigg)\bigg]+\bigg[\mfb+\frac{\mfc}{f_{0}^{2}}+\frac{3\ \mfd}{4 f_{0}^{1/4}}\\ &+\frac{\mfe\, \fdamp}{\fdamp^{2}+(f_{0}-f_{\RD})^{2}}\bigg]f. \end{align} Let us now define the constant term of the phase as \begin{align}\notag \psi_{1}\equiv\mfa-&2\frac{\mfc}{f_{0}}+\frac{1}{4}\mfd \,f_{0}^{3/4}-\frac{\mfe f_{0} \fdamp}{\fdamp^{2}+(f_{0}-\mfg f_{\RD})^{2}}\\ &+\mfe \arctan\bigg(\frac{f_{0}-\mfg\, \frd}{\fdamp}\bigg) \end{align} and the linear coefficient as \be \psi_{2}\equiv\mfb+\frac{\mfc}{f_{0}^{2}}+\frac{3\ \mfd}{4 f_{0}^{1/4}}+\frac{\mfe\, \fdamp}{\fdamp^{2}+(f_{0}-f_{\RD})^{2}} \ee such that the approximated phase can neatly be written as $\psi(f)\simeq\psi_{1}+f\, \psi_{2}$. Our inverse Fourier transform then becomes \begin{align} h(t)&=\frac{1}{2 \pi}\int\scA(f)\, e^{-i(\psi_{1}+f\, \psi_{2})}\, e^{2 \pi i f t}df\\\notag &=\frac{1}{2 \pi}e^{-i\psi_{1}}\int\scA(f)\, e^{-i\,f\, \psi_{2}}\, e^{2 \pi i f t}df\\\notag &=\frac{1}{2 \pi}e^{-i\psi_{1}}\int\scA(f)\, e^{-if(\psi_{2}-2 \pi t)}df \end{align} which includes an easily integrable exponential function. Now, we must concern ourselves with the amplitude portion of the transform. We can pull out the exponential that is constant in frequency such that \begin{align} &h(t)=\frac{1}{2 \pi}e^{-i\psi_{1}}\int\scA(f)\, e^{-if(\psi_{2}-2 \pi t)}df\\\notag &=\frac{1}{2 \pi}e^{-i\psi_{1}}\int\frac{\alpha\, f^{-7/6}}{(f-f_{\RD})^{2}+\beta^{2}}\, e^{-\mu (f-\frd)}\, e^{-if(\psi_{2}-2 \pi t)}df\\\notag &=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}}\int\frac{\, f^{-7/6}}{(f-f_{\RD})^{2}+\beta^{2}}\, e^{-\mu f}\, e^{-if(\psi_{2}-2 \pi t)}df\\\notag \end{align} which is looking more simple already. We can replace the amplitude term with the Taylor series approximation. (This looks like the best approximation when the expansion goes to 8th order in $f-f_{0}$. We can look at the errors that we get by decreasing or increasing the order when we are comparing how well this model fits to numerical models.) Then, \begin{align} &h(t)=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}}\int\frac{\, f^{-7/6}\, e^{-\mu f}\, e^{-if(\psi_{2}-2 \pi t)}df}{(f-f_{\RD})^{2}+\beta^{2}}\, \\\notag &=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}}\int\sum_{k} \sigma_{k} (f-f_{0})^{k}\, e^{-\mu f}\, e^{-if(\psi_{2}-2 \pi t)}df\\\notag &=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}}\sum_{k} \sigma_{k}\int (f-f_{0})^{k}\, e^{-\mu f}\, e^{-if(\psi_{2}-2 \pi t)}df\\\notag \end{align} where k is the order to which we evaluate the Taylor series. This is essentially the final form of our integral. Note that all of the following are constants that can be evaluated \textit{in the time domain} without having to do any manipulation in the frequency domain: ${\alpha, \psi_{1,2},\mu,\frd,\mbox{ and }\sigma_{k}}$. Then, this final integral is analytic for all $k$. In order to implement Feynman's trick here, we must convert the integral to one of the form $\int(f-f_{0})^{k}\,e^{\nu(f-f_{0})}df$. Thus, in order to simplify out integral at least a little, we can rewrite it as \begin{align} h(t)&=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}}\sum_{k} \sigma_{k}\int (f-f_{0})^{k}\, e^{-\mu f}\, e^{-if(\psi_{2}-2 \pi t)}df\\\notag &=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}}\sum_{k} \sigma_{k}\int (f-f_{0})^{k}\, e^{f(-\mu -i\psi_{2}+2 \pi i t)}df\\\notag &=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}}\sum_{k} \sigma_{k}\int (f-f_{0})^{k}\, e^{\nu f}df\\\notag &=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}}\sum_{k} \sigma_{k}\int (f-f_{0})^{k}\, e^{\nu f}e^{-\nu f_{0}+\nu f_{0}}df\\\notag &=\frac{\alpha}{2 \pi}e^{-i\psi_{1}}e^{\mu\,f_{\RD}+\nu f_{0}}\sum_{k} \sigma_{k}\int (f-f_{0})^{k}\, e^{\nu (f-f_{0})}df\\\notag \end{align} which is exactly of the form needed for Feynman's Trick. This is all implemented analytically in Mathematica. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \nocite{*} \bibliography{bibliography.bib} \end{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% %% %