\documentclass[url,11pt]{ligovirgodcc}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{longtable}
\usepackage{rotating}
\usepackage{color}
\usepackage{epsf}
\usepackage{fancyhdr}
\usepackage{ifthen}
\usepackage{setspace}
\usepackage[comma,numbers,sort&compress]{natbib}
\setlength{\parindent}{0.0cm}
\newcommand{\Ys}{{{}^{-s}Y}}
\newcommand{\Ytwo}{{{}^{-2}Y}}
\def\pasp{Publ. Astron. Soc. Pac.}
\def\prc{Phys. Rev. C.}
\def\prd{Phys. Rev. D.}
\def\aa{Astron. Astrophys.}
\def\aj{Astron. J.}
\def\araa{Ann. Rev. Astron. Astrophys.}
\def\iaucirc{Int. Astron. U. Circ.}
\def\aap{Astron. Astrophys.}
\def\aaps{Astron. Astrophys. Suppl.}
\def\apss{Astroph. \& Space Sc.}
\def\apj{Astrophys. J.}
\def\apjl{Astrophys. J. Lett.}
\def\apjs{Astrophys. J. Suppl. Ser. }
\def\mnras{Mon. Not. R. Astron. Soc.}
\def\physrep{Phys. Rep.}
\def\pasj{Publ. Astron. Soc. Jap.}
\def\nat{Nature}
\def\prl{Phys. Rev. Lett.}
\def\azh{Soviet Astron.}
\newcommand{\todo}[1]{{\color{blue}$\blacksquare$~\textsf{\small{[TODO: #1]}}\color{black}}}
\newcommand{\todored}[1]{{\color{red}$\blacksquare$~\textsf{[TODO: #1]}\color{black}}}
%\documentclass[12pt,letterpaper]{article}
%\setlength{\oddsidemargin}{-0.25in}
%\setlength{\evensidemargin}{-0.25in}
%\setlength{\topmargin}{-0.25in}
%\setlength{\textwidth}{7.0in}
%\setlength{\textheight}{8.5in}
%%%%%%%%%%%
\begin{document}
\pagestyle{fancy}
\headheight 14pt
\chead[{\large GWs from Accretion Disk Instabilities}]
{{\large GWs from Accretion Disk Instabilities}}
\rhead[]{}
\lhead[]{}
%
\title{\large Gravitational Wave Emission from Accretion Disk Instabilities -- Analytic Models}
%\title{\large LIGO Summer Undergraduate Research Report}
\ligodocdraft{\hfill}
\ligodoc{T1100093}{v2}
\virgodoc{}
\author{Luc\'ia Santamar\'ia (LIGO Lab, Caltech, {\tt luciasan@caltech.edu}),\\ Christian D. Ott (TAPIR/CaRT, Caltech, {\tt cott@tapir.caltech.edu})}
\date{August 31, 2013}
\ligoabstract{\noindent We derive the gravitational wave emission from
accretion disk instabilities in long gamma-ray bursts or black-hole
forming core-collapse supernovae based on simple analytic models of
van Putten and Piro \& Pfahl.}
\maketitle
%\tableofcontents
\pagebreak
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{s:intro}
We derive analytic estimates for the gravitational wave (GW)
emission from the following analytic accretion disk instability
models in the context of Collapsar-type gamma-ray bursts and black-hole (BH)
forming core-collapse supernovae:
\begin{itemize}
\item {\bf The suspended-accretion-driven disk instability proposed by van
Putten} \cite{mvp:04,mvp:01}.\\ In this instability, turbulence in a
thick accretion torus is driven by strong coupling between the torus
and the ergosphere of the central BH by MHD effects. In this
picture, the quadrupole components of the disk turbulence lead to GW
emission that (via the strong coupling) spin down the BH. The
strong coupling and the energetics of the associated GW emission
proposed by van Putten are generally regarded as too optimistic
(possibly by orders of magnitude). Nevertheless, we include this
model, since it makes useful, falsifiable predictions.
\emph{We advise that authors wishing to use the waveforms predicted
by our ad-hoc version of van Putten's model state clearly to
their readers that it is unlikely that a realistic accretion disk
instability will produce a gravitational wave signal of the kind
predicted by the model implemented here.}
\item {\bf Collapsar Disk Fragmentation by Piro \& Pfahl}
\cite{piro:2007}.\\ Gravitational instability causes fragmentation
of the massive accretion disk within a collapsing, rotating
star. This fragment spirals into the newly-formed BH at the center
of the accretion disk due to a combination of disk viscosity and GW
emission, producing a unique GW signature.
\end{itemize}
In the following, we work in physical units, including all relevant
factors of the gravitational constant $G$ and of the speed of light
$c$. Accompanying this technical report are two python scripts, {\tt
vanPuttengw.py} and {\tt pirogw.py}, that implement the models
described here and provide both GW polarizations as output.
\section{Suspended-Accretion Quadrupole Disk Instability}
Following along the lines of van Putten's ideas, we assume a spinning
Kerr BH with dimensionless Kerr spin parameter $a^\star = (c/G)
J_\mathrm{BH}/M_\mathrm{BH}^2$ ($0 \leq a^\star < 1$) with an
accretion disk/torus of mass $M_\mathrm{disk}$. The disk extends to
the radius of the innermost stable orbit \cite{bardeen:72},
%
\begin{align}
\label{eq:risco}
R_\mathrm{ISCO} &= \left(\frac{G}{c^2}\right) M_\mathrm{BH}\left( 3 + Z_2 \mp \left[(3-Z_1)(3+Z_1+2Z_2) \right]^{1/2}\right), \\
Z_1 & = 1 + (1 - {a^\star}^2)^{1/3}\left[ (1+a^\star)^{1/3} + (1-a^\star)^{1/3} \right], \\
Z_2 & = \left[3{a^\star}^2 + Z_1^2\right]^{1/2},
\end{align}
%
where the $\mp$ sign indicates prograde and retrograde orbits,
respectively. We will assume that the binary orbits in the same
direction as the BH spin and thus take the minus sign in
equation~\ref{eq:risco}.
Disk and BH are assumed to be coupled via strong magnetic fields and
MHD turbulence in the disk is assumed to be driven through this
coupling. The inner disk near the ISCO is expected to be neutrino
cooled and very thin. Further out, at $r_0 + R_\mathrm{ISCO}$ (where
$r_0 = 100$ km may be reasonable), the disk is a thick torus. We
assume that turbulence in the torus leads to two overdense regions
(which, in itself, is unlikely, since turbulent power will cascade to
small scales) with masses $M_1 = M_2 = \epsilon M_\mathrm{disk}$.
$\epsilon \approx 0.01 - 0.5$ (the latter is very unlikely). These two
`clumps' form a `binary' with separation
$2(r_0+R_\mathrm{ISCO})$ that efficiently emits GWs and would normally
lose $J$ by GW emission, leading to inspiral. The BH is located at the
center of our coordinate system and each of the clumps is located at a
distance $r_0+ R_\mathrm{ISCO}$ from the BH. Here, following van
Putten, we assume that the lost energy and angular momentum is
replenished by coupling to the central BH, so the BH loses $J$ and
spin energy. This leads to an incremental change of $R_\mathrm{ISCO}$
and, consequentially, of the binary separation
$2(r_0+R_\mathrm{ISCO})$.
\vskip.5cm
\subsection{Gravitational Wave Emission using the Newtonian Binary Approximation}
\label{sec:vanputtengw}
We assume, without loss of generality, that the BH angular momentum is
oriented in the $+z$-direction. Two overdense regions in the torus are
approximated as point masses $M_1$ and $M_2$, where we assume $m = M_1 = M_2$.
The two point masses orbit the BH at the origin in the
$xy$-plane with angular velocity $\Omega$. The orbital radius is $d
\equiv r_0+R_{\mathrm{ISCO}}$. We can compute the reduced mass
quadrupole momentum of the system, defined as
%
\begin{equation}
I_{ij} = m \sum_A (x_{Ai}x_{Aj} - \frac{1}{3}\delta_{ij}d^2),
\label{eq:quadbasic}
\end{equation}
where $A$ is an index denoting overdense regions.
Assuming that fragment 1 sits on the positive $x$ axis at $t=0$,
the time-dependent positions $(x_i,y_i)$ of the two fragments
are described by
\begin{equation}
\begin{array}{lcl}
x_1 &=& d \cos(\Omega t)\,\,,\\
y_1 &=& d \sin(\Omega t)\,\,,
\end{array} \hspace*{1cm}\begin{array}{lcl}
x_2 & =& - d \cos(\Omega t)\,\,,\\
y_2 & =& - d \sin(\Omega t)\,\,.
\end{array}
\end{equation}
The quadrupole moment components we are interested in
are $I_{xx}$, $I_{xy}$, and $I_{yy}$. Because of symmetry,
$I_{xy} = I_{yx}$. Using eq.~(\ref{eq:quadbasic}), we have
\begin{equation}
\begin{aligned}
I_{xx} & = m (x_1^2 + x_2^2 - \frac{2}{3}d^2)\,\,,\\
&= 2 m d^2 \cos^2{\Omega t} - \frac{2}{3}d^2\,\,,\\
&= m d^2 \cos{2\Omega t} + const.
\end{aligned}
\end{equation}
Similiarly, we obtain the remaining components:
\begin{equation}
I_{ij} =
m d^2 \left(\begin{array}{cc} \cos{2\Omega t} + const. &\sin{2 \Omega t} + const\\
\sin{2 \Omega t} + const. &-\cos{2 \Omega t} + const.
\end{array}\right),
\end{equation}
\begin{equation}
\ddot{I}_{ij} =
4 m d^2 \Omega^2 \left(\begin{array}{cc}
-\cos{2\Omega t}&-\sin{2 \Omega t}\\
-\sin{2 \Omega t}&\cos{2 \Omega t}
\end{array}\right),
\end{equation}
\begin{equation}
\dddot{I}_{ij} =
8 m d^2 \Omega^3 \left(\begin{array}{cc}
\sin{2\Omega t}&-\cos{2 \Omega t}\\
-\cos{2 \Omega t}&-\sin{2 \Omega t}
\end{array}\right).
\end{equation}
%
The gravitational wave signal emitted by the binary system and the change
in angular momentum and energy are given by
%
\begin{equation}
h^{TT}_{ij} = \frac{2}{D} \frac{G}{c^4} \ddot{I}^{TT}_{kl}\,\,,\qquad
\frac{dJ_i}{dt} = -\frac{2}{5}\frac{G}{c^5}\epsilon_{ijk}
\langle\ddot{I}_{jm}\dddot{I}_{mk}\rangle\,\,, \qquad
\frac{dE}{dt} = -\frac{1}{5}\frac{G}{c^5}\langle \dddot{I}_{jk}
\dddot{I}_{jk}\rangle\,\,.
\end{equation}
The change of orbital energy and angular moment straightforwardly
evaluate to
\begin{equation}
\label{eq:jdot}
\frac{dJ_{\mathrm{GW}}}{dt} = -\frac{128}{5} \frac{G}{c^5} m^2 d^4 \Omega^5 \,\,,
\end{equation}
%
\begin{equation}
\label{eq:edot}
\frac{dE_{\mathrm{GW}}}{dt} = P_{\mathrm{GW}} = -\frac{128}{5} \frac{G}{c^5} m^2 d^4
\Omega^6\,\, .
\end{equation}
Next, we must carry out the transverse-traceless projection of
$\ddot{I}$ to obtain the two GW polarizations $h+$ and
$h_\times$. These depend on the position (polar angle $\theta$
azimutal angle $\phi$ and distance $D$) the observer has with respect
to the source. This is done most easily by expressing the GW strain in
terms of spin-weighted tensor spherical harmonics,
\begin{equation}
h_+ - ih_\times = \frac{1}{D}\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^\ell H_{\ell m}(t)\,
\Ytwo_{\ell m}(\theta,\phi)\,.
\end{equation}
The expansion parameters $H_{lm}$ ($\ell = 2$ denotes the quadrupole)
are complex functions of the retarded source time $t$.
In order to express $H_{2m}$ in terms of $\ddot{I}_{ij}$, one first
expresses $h_+(\theta,\phi)$ and $h_\times(\theta,\phi)$ in terms of
$\ddot{I}_{kl}$, then convolves these with $^{-2}Y^*_{lm}$. The result
\cite{ajith:07} is
\begin{eqnarray}
\label{eq:Hlms1}
H^\mathrm{quad}_{20} &=& \sqrt{\frac{32\pi}{15}} \frac{G}{c^4}
\left(\ddot{I}_{zz}
- \frac{1}{2}(\ddot{I}_{xx}
+\ddot{I}_{yy})\right)\,\,,\\
\label{eq:Hlms2}
H^\mathrm{quad}_{2\pm1} &=& \sqrt{\frac{16\pi}{5}} \frac{G}{c^4}
\left(\mp \ddot{I}_{xz}
+ i \ddot{I}_{yz}\right)\,\,,\\
H^\mathrm{quad}_{2\pm2} &=& \sqrt{\frac{4\pi}{5}} \frac{G}{c^4}
\left(\ddot{I}_{xx}
- \ddot{I}_{yy} \mp
2 i \ddot{I}_{xy} \right)\,\,.
\label{eq:Hlms3}
\end{eqnarray}
For completeness, we give the definitions of the relevant
$^{-2}Y^*_{lm}$:
\begin{eqnarray}
\label{eq:7}
\Ytwo_{22} &=& \sqrt{\frac{5}{64\pi}}(1+\cos\theta)^2e^{2i\phi} \,,\\
\Ytwo_{21} &=& \sqrt{\frac{5}{16\pi}} \sin\theta( 1 + \cos\theta )e^{i\phi} \,,\\
\Ytwo_{20} &=& \sqrt{\frac{15}{32\pi}} \sin^2\theta \,,\\
\Ytwo_{2-1} &=& \sqrt{\frac{5}{16\pi}} \sin\theta( 1 - \cos\theta
)e^{-i\phi} \,,\\
\Ytwo_{2-2} &=& \sqrt{\frac{5}{64\pi}}(1-\cos\theta)^2e^{-2i\phi}\,.
\end{eqnarray}
For example, assuming the detector is located along the positive $z$
axis ($\theta = 0, \phi = 0$), we find:
\begin{eqnarray}
h_+ &=& \frac{1}{D} \frac{G}{c^4} (\ddot{I}_{xx} - \ddot{I}_{yy})\\
h_\times &=& \frac{2}{D} \frac{G}{c^4} \ddot{I}_{xy}\,\,.
\end{eqnarray}
\subsection{The coupled system of Ordinary Differential Equations}
We assume that the `binary' stays at a fixed radius. Angular
momentum $J$ lost to GW emission is provide from the BH spin. As a
consequence, the BH is spun down and $R_\mathrm{ISCO}$ changes. We
set
$\dot{J}_{\mathrm{BH}}=\dot{J}_{\mathrm{GW}}$. The change in the BH
mass is
\begin{equation}
\label{eq:mdot}
\dot{M}_{\mathrm{BH}}=\frac{P_{\mathrm{GW}}}{c^2},
\end{equation}
since the gravitational mass of the BH contains the contribution due
to its rotational energy. The change in $R_{\mathrm{ISCO}}$ can then
be computed by differentiating equation~\ref{eq:risco},
%
\begin{align}
\label{eq:riscodot}
\frac{dR_\mathrm{ISCO}}{dt}=\dot{R}_{\mathrm{ISCO}} &=
\left(\frac{G}{c^2}\right) \left[ \dot{M}_{\mathrm{BH}}
\left( 3 + Z_2 - \sqrt{(3-Z_1)(3+Z_1+2Z_2)} \right)
\right. \notag \\
& \qquad \qquad \qquad \left. + M_{\mathrm{BH}}
\left(\dot{Z}_2 +
\frac{(Z_1+Z_2)\dot{Z}_1- (3-Z_1)\dot{Z}_2}{\sqrt{(3-Z_1)(3+Z_1+2Z_2)}} \right)
\right],
\end{align}
%
where we can express $Z_1,Z_2$ and their derivatives as functions
of $J_{\mathrm{BH}}$, $M_{\mathrm{BH}}$, $\dot{J}_{\mathrm{BH}}$,
$\dot{M}_{\mathrm{BH}}$ only,
%
\begin{align}
\label{eq:zdot1}
\dot{Z}_1 &= \frac{c \left(M_{\mathrm{BH}}
\dot{J}_{\mathrm{BH}}-2 J_{\mathrm{BH}} \dot{M}_{\mathrm{BH}}\right)}
{3 G^3 M_{\mathrm{BH}}^7 \left( 1 - {a^\star}^2\right)^{4/3}}
\left[ 3 c^2 J_{\mathrm{BH}}^2 \left( \left(1 + a^\star \right)^{2/3} -
\left(1 - a^\star \right)^{2/3} \right) \right. \notag \\
& \left. - 2 c G J_{\mathrm{BH}} M_{\mathrm{BH}}^2
\left( \left(1 + a^\star \right)^{2/3} + \left(1 - a^\star \right)^{2/3}\right)
+ G^2 M_{\mathrm{BH}}^4
\left( \left(1 - a^\star \right)^{2/3} - \left(1 + a^\star \right)^{2/3}\right)
\right]\,\,,
\end{align}
%
\begin{equation}
\label{eq:zdot2}
\dot{Z}_2=\frac{3 c^2 J_{\mathrm{BH}} M_{\mathrm{BH}} \dot{J}_{\mathrm{BH}}
- 6 c^2 J_{\mathrm{BH}}^2 \dot{M}_{\mathrm{BH}}
+ G^2 M_{\mathrm{BH}}^5 Z_1 \dot{Z}_1}
{G^2 M_{\mathrm{BH}}^5 \sqrt{\frac{3 c^2 J_{\mathrm{BH}}^2}{G^2 M_{\mathrm{BH}}^4}
+ Z_1^2}}\,\,.
\end{equation}
These expressions allows us to calculate the change in
$R_\mathrm{ISCO}$ when the mass and angular momentum of the central BH
change.
\subsection{Application of the coupled system}
At every time $t$, $\ddot{I}_{ij}$ and, thus, $h^{TT}_{ij}$ depend on
(\emph{i}) the orbital radius of the `binary', $d = r_0+
R_\mathrm{ISCO}$, (\emph{ii}) the mass $m$ of the chunks (assumed to
be constant), and (\emph{iii}) on the angular velocity, which,
according to Kepler's law, we set to
%
\begin{equation}
\label{eq:omega}
\Omega = \sqrt{\frac{G M}{d^3}}\,\,.
\end{equation}
%
Let us assume that the mass of the chunks forming the `binary' is
negligible with respect to the mass of the central black hole, in
which case $\Omega = \sqrt{G M_{\mathrm{BH}}/(r_0+
R_{\mathrm{ISCO}})^3}.$ The coupled system of ODEs is then formed by
equations~\ref{eq:jdot}, \ref{eq:mdot} and \ref{eq:riscodot}, with $d
= r_0 + R_\mathrm{ISCO}$ and $\Omega$ as in
equation~\ref{eq:omega}. $P_{\mathrm{GW}}$ in equation~\ref{eq:mdot}
is given by equation~\ref{eq:edot}. This system describes the
evolution of $R_{\mathrm{ISCO}}, J_{\mathrm{BH}}$ and
$M_{\mathrm{BH}}$. We integrate the coupled system of ODEs with a
fourth-order Runge-Kutta integrator.
Note that once all spin has been extracted from the hole, the `binary'
will inspiral. However, we stop our integration when the BH is
completelly spun down and do not calculate the subsequent chirp.
\subsection{Astrophysically Meaningful Parameters}
\label{subsec:MvPpars}
\begin{tabular}{rl}
Mass & BH of mass $M_{\mathrm{BH}} = 5-10\, M_\odot$ \\
Initial BH spin & $a^\star = 0.3-0.95$ \\
Fragment mass & Assume $M_{\mathrm{disk}} = 1.5 M_{\odot}$,
$m_{\mathrm{chunk}} = \epsilon M_{\mathrm{disk}}$ with $\epsilon =
0.01-0.2$ \\
Const. separation of torus from ISCO & $r_0 = 100$ km \\
End integration & When $J_{\mathrm{BH}} = 0$ or pre-specified run time is reached\\
\end{tabular}
\subsection{Usage of the python script}
{\it Code usage:} {\tt \$ python vanPutten.py}
The code section below {\tt \#physical parameters} allows to specify the
physical parameters
of the system as defined in subsection~\ref{subsec:MvPpars}.
The colatitude and azimuth of the system can be specified as well.
The code
section under {\tt \#parameters} allows to change the total run
time (in seconds) of the integration (provided that the BH is not
completely spun down, in which case the integration stops) as well as
the sampling time {\tt dt} of the output.
The script produces two output files:
\begin{itemize}
\item {\tt pmvp.dat} is a diagnosis and debug output file, containing
the following variables: \\
$\boxed{\mathrm{time} \qquad R_\mathrm{ISCO} \qquad J_{\mathrm{BH}} \qquad
M_{\mathrm{BH}} \qquad a^\star \qquad E_{\mathrm{rad}} \qquad h_+
\qquad h_\times}$
\item {\tt M*a*eps*.dat}, where the {\tt *}'s denote the values of the physical parameters
$M_{\mathrm{BH}}$, $a^\star$ and $\epsilon$ given as input,
is the production output file, containing:\\
$\boxed{\mathrm{time} \qquad h_+ \qquad h_\times }$
\end{itemize}
An example of the antichirp-like signal obtained for the van Putten
model is shown in figure~\ref{fig:vp}.
\begin{figure}
\includegraphics[width=\linewidth]{vPwf.pdf}
\caption{Waveform computed following the van Putten model. The
parameters of the system are $M_{\mathrm{BH}}= 10 M_\odot$, $a^\star
= 0.95$, $\epsilon = 0.2$. The strain corresponds to a face-on,
optimally-oriented system situated at 10 kpc. The main plot shows
the absolute magnitude of the strain $|h| = \sqrt{h_+^2 +
h_\times^2}$ and the inset plot shows the two polarizations
zooming into the first 0.1 s of evolution.}
\label{fig:vp}
\end{figure}
\section{Torus Fragmentation Instability and Inspiral of a Single Overdense Blob}
When the core-collapse supernova mechanism fails to re-energize the
stalled shock (see, e.g., \cite{janka:07}), the protoneutron star
collapses to a BH on a timescale of $\sim$$500\,\mathrm{ms} -
\mathrm{few}\,\,\mathrm{s}$~\cite{oconnor:11a}. Provided that the star
has sufficient angular momentum, its collapse will eventually be
halted and a massive accretion disk/torus will form around the nascent
stellar-mass BH. This collapsar scenario is the leading model for long
duration gamma-ray bursts and may also power an ``engine driven''
supernova \cite{wb:06}.
The inner regions of the accretion disk are geometrically thin due to
neutrino cooling, but the outer parts are unable to cool efficiently
causing them to be thick. This potentially leads to fragmentation at
large radii. We investigate the expected gravitational radiation from
such a system, inspired by the discussion by Piro \& Pfahl~\cite{piro:2007}.
\vskip.2cm
We assume a central BH of mass $M_{\mathrm{BH}}$ surrounded by a
Keplerian accretion disk with orbital frequency
$\Omega_\mathrm{disk}=(G M_{\mathrm{BH}}/r^3)^{1/2}$ and vertical
scale height $H$. The accretion rate is $\dot{M} = 3 \pi \nu \Sigma$,
where $\Sigma$ is the disk's surface density, $\nu = \alpha c_S H$ the
usual $\alpha$-viscosity prescription and $c_s$ the isothermal speed
of sound. We assume that $H/r = \eta$ is a fixed parameter.
\subsection{Gravitationaly instability and fragmentation}
Gravitational instability arises when
\begin{equation}
Q \equiv \frac{\Omega_\mathrm{disk} c_s}{\pi G \Sigma} <
Q_{\textrm{crit}} \simeq 1\,\,.
\end{equation}
Even once a disk becomes gravitationally unstable, numerical
simulations demonstrate that fragmentation is only promoted if there
is a sufficiently rapid cooling mechanism available to allow collapse
to high densities. Piro \& Pfahl \cite{piro:2007} argue that
photodisintegration of $^4 \textrm{He}$ naturally leads to
fragmentation because this endothermic reaction removes energy at a
rate that far exceeds the viscous dissipation in the disk. In
addition, the photodisintegration is very temperature sensitive and
sets a clear radius were it first becomes active at around $100$
gravitational radii $(R_g = GM/(rc^2))$. As explained in
\cite{piro:2007}, we can therefore identify the quantity $(QH)^2 \Sigma$
with the mass of the bound fragment, which is estimated to be
\begin{equation}
\label{eq:massfrag}
M_f \approx 0.2 \left(\frac{\eta}{0.5} \right)^3 \frac{M_\textrm{BH}}{3}\,\,.
\end{equation}
\subsection{Migration and associated gravitational waves}
Due to it's relatively large mass in comparison to surrounding disk
material, the fragment accretes material within a tidal radius of
itself and opens a gap in the accretion disk. It then migrates inward
toward the central BH due to tidal interactions with the edges of the
viscously accreting disk in a well-known process called Type II
migration. (This is opposed to Type I migration, which takes places
when a fragment is too small to open a gap and instead migrates inward
via the generation of density waves in the disk.) This viscous
migration happens on the approximate viscous timescale given by Piro \&
Pfahl~\cite{piro:2007},
\begin{equation}
t_\nu \approx \frac{1}{\alpha\eta^2\Omega}\,\,.
\end{equation}
Note that since the mass of the fragment is not necessarily
negligible, we use here $\Omega = \sqrt{G (M_\mathrm{BH}+M_f)
r^{-3}}$, whereas $\Omega_\mathrm{disk} = \sqrt{GM_\mathrm{BH}
r^{-3}}$.
Since the fragment-BH system is effectively a binary, there is
associated GW emission that can also contribute to the inward
migration of the fragment. The timescale for this can be estimated to
first order via the quadrupole formula to be
\begin{equation}
t_\textrm{GW} = \frac{5}{64 \Omega} \left( \frac{G
\mathcal{M}\Omega}{c^3} \right)^{-5/3},
\end{equation}
where $\mathcal{M} = \mu^{3/5} (M_\mathrm{BH} + M_f)^{2/5})$ is the
chirp mass and $\mu = M_\mathrm{BH} M_f / (M_\mathrm{BH} + M_f)$ is
the reduced mass of the system. The evolution of the orbit under the
influence of both GW and viscous effects can be computed by simply
solving the differential equation
\begin{equation}
\label{eq:drdt}
\frac{dr}{dt} = -r \left(\frac{1}{t_\textrm{GW}} + \frac{1}{t_v} \right)\,\,,
\end{equation}
Generally the viscous processes dominate in causing migration at large
radii, while GW emission dominates at small radii. Together this leads
to a unique GW signature from the fragment-BH system. This may allow
the physics of the disk to be probed via the GW signal by identifying
when the inspiral switches from being predominantly controlled by
viscosity to being controlled by GWs.
\vskip.2cm
We integrate Eq.~(\ref{eq:drdt}) with a standard fourth-order
Runge-Kutta integrator. The chirp-like gravitational wave
emission expected from the system via the quadrupole approximation
as discussed in the next subsection.
\subsection{Gravitational Wave Signal}
As in \S\ref{sec:vanputtengw}, we consider two point masses,
$M_\mathrm{BH}$ and $M_f$, separated by distance $r$. The total mass
is $M = M_\mathrm{BH} + M_f$ and the reduced mass is $\mu =
M_\mathrm{BH}M_f / M$. The motion (arbitrarily assumed to be in the
$x-y$ plane) of the two masses about their center of mass (at the
origin) is described by the following equations:
\vspace*{-0.75cm}
{\setstretch{2.25}
\begin{equation}
\begin{array}{lcl}
x_\mathrm{BH} &=& \dfrac{M_\mathrm{BH}}{M} r \cos(\Omega t)\,\,,\\
y_\mathrm{BH} &=& \dfrac{M_\mathrm{BH}}{M} r \sin(\Omega t)\,\,,
\end{array} \hspace*{1cm}\begin{array}{lcl}
x_f & =& - \dfrac{M_f}{M} r \cos(\Omega t)\,\,,\\
y_f & =& - \dfrac{M_f}{M} r \sin(\Omega t)\,\,.
\end{array}
\end{equation}
\setstretch{1.00} } Here, $\Omega = \sqrt{G (M_\mathrm{BH}+M_f)
r^{-3}}$. Now, using Eq.~(\ref{eq:quadbasic}), we have
\begin{equation}
I_{ij} =
\frac{\mu}{2} r^2 \left(\begin{array}{cc} \cos{2\Omega t} + const. &\sin{2 \Omega t} + const\\
\sin{2 \Omega t} + const. &-\cos{2 \Omega t} + const.
\end{array}\right),
\end{equation}
\begin{equation}
\ddot{I}_{ij} =
2 \mu r^2 \Omega^2 \left(\begin{array}{cc} -\cos{2\Omega t} &-\sin{2 \Omega t} \\
-\sin{2 \Omega t} &\cos{2 \Omega t}
\end{array}\right)\,\,,
\end{equation}
\begin{equation}
\dddot{I}_{ij} =
4 \mu r^2 \Omega^3 \left(\begin{array}{cc} \sin{2\Omega t} &-\cos{2 \Omega t} \\
-\cos{2 \Omega t} &-\sin{2 \Omega t}
\end{array}\right)\,\,.
\end{equation}
The source-angle-dependent $h_+$ and $h_\times$ GW polarizations can then be
obtained via Eqs.~(\ref{eq:Hlms1}--\ref{eq:Hlms3}).
\subsection{Astrophysically Meaningful Parameters}
\label{subsec:piropars}
\begin{tabular}{rl}
Mass & BH of mass $M_{\mathrm{BH}} = 3-10 M_\odot$ \\
Viscosity & Standard value of $\alpha = 0.1$ \cite{piro:2007}\\
Geometrical parameter & $\eta = 0.3 -0.6$ \cite{piro:2007}\\
Mass of the bound fragment & The approximate factor 0.2 in
Eq.~(\ref{eq:massfrag}) can be varied from $0.2-0.5$.\\
Starting Radius & Start at $r = 100 r_g$, where $r_g =
\frac{GM}{c^2}$ is the gravitational radius \\
& End integration close to the ISCO \\
\end{tabular}
\subsection{Usage of the python script}
{\it Code usage:} {\tt \$ python pirogw.py}
The code section below {\tt \#physical parameters} allows to specify the
physical parameters
of the system as defined in subsection~\ref{subsec:piropars}.
The colatitude and azimuth of the system can be specified as well.
The code
section under {\tt \#parameters} allows to change the total run
time (in seconds) of the integration (provided that the orbital radius
is larger than $2.5 R_{\mathrm{ISCO}}$ where we take a multiple of the ISCO radius
of a non-spinning black hole $R_{\mathrm{ISCO}} = 6 G
M_{\mathrm{BH}}/c^2$ as a limit for the integration) as well as
the sampling time {\tt dt} of the output.
The script produces two output files:
\begin{itemize}
\item {\tt piro.dat} is a diagnosis and debug output file, containing
the following variables: \\
$\boxed{\mathrm{time} \qquad r\, \mathrm{(cm)} \qquad r\, (r_S) \qquad
h_+ \qquad h_\times \qquad \Omega \qquad t_{GW} \qquad t_v }$
\item {\tt piroM*eta*fac*.dat}, where {\tt *} denote the values of the physical parameters
$M_{\mathrm{BH}}$, $\eta$ and the factor in the RHS of
equation~\ref{eq:massfrag} given as input, is the production output file, containing:\\
$\boxed{\mathrm{time} \qquad h_+ \qquad h_\times }$
\end{itemize}
%%% \subsection{Gravitationaly instability and fragmentation}
%%%
%%% Gravitational instability arises when
%%%
%%% \begin{equation}
%%% Q \equiv \frac{\Omega c_S}{\pi G \Sigma} < Q_{\textrm{crit}} \simeq 1\,\,,
%%% \end{equation}
%%%
%%% Under given circunstances, numerical
%%% simulations have shown that gravitational instability leads to
%%% fragmentation. As explained in~\cite{piro:2007}, we can identify the quantity
%%% $(QH)^2\Sigma $ with the mass of a bound clump if cooling is rapid
%%% enough to permit collapse to high densities.
%%%
%%% Various possible cooling mechanisms in the disk can be studied, among
%%% them radiative diffusion, neutrino cooling processes, and
%%% photodisintegration. Photodisintegration of $^4 \textrm{He}$ will
%%% absorb sufficient energy quickly enough to allow fragmentation. Thus,
%%% it may be a very effective coolant, permitting fragmentation over
%%% a small range of radii near the location where $Q = Q_{\textrm{crit}}$
%%% in the perturbed disk. This allows us to estimate the mass of the bound
%%% fragment as
%%%
%%% \begin{equation}
%%% \label{eq:massfrag}
%%% M_f \approx 0.2 \left(\frac{\eta}{0.5} \right)^{\!3} \frac{M_\textrm{BH}}{3}\,\,.
%%% \end{equation}
%%%
%%% \subsection{Migration and associated gravitational waves}
%%%
%%% If the fragment is massive enough to form a bound object and open a
%%% gap in the accretion disk,
%%% it migrates inwards, forming an inspiralling binary with the central
%%% BH. The loss of angular momentum is due to dissipation within the disk
%%% and to the emission of gravitational waves. The viscous migration happens
%%% on a the viscous timescale,
%%%
%%% \begin{equation}
%%% t_\nu \approx \frac{1}{\alpha\nu^2\Omega}\,\,.
%%% \end{equation}
%%%
%%% Loss of angular momentum in the form of gravitational waves causes
%%% inspiral on a time that can be estimated to first order via the
%%% quadrupole formula
%%%
%%% \begin{equation}
%%% t_\textrm{GW} = \frac{5}{64 \Omega} \left( \frac{G
%%% \mathcal{M}\Omega}{c^3} \right)^{-5/3},
%%% \end{equation}
%%%
%%% where $\mathcal{M}$ is the chirp mass of the system formed by BH plus
%%% fragment. The evolution of the orbit can be computed by simply solving
%%% the differential equation
%%%
%%% \begin{equation}
%%% \label{eq:drdt}
%%% \frac{dr}{dt} = -r \left(\frac{1}{t_\textrm{GW}} + \frac{1}{t_v} \right)
%%% \end{equation}
%%%
%%% The integration of equation~\ref{eq:drdt} allows us to compute the chirp-like
%%% gravitational wave emission expected from this system.
%%% We make use of equation~\ref{eq:hTT} where $D$ is the distance between
%%% the central BH and the bound fragment.
%%% The ODE integration is implemented with a fourth-order Runge-Kutta integrator.
%%%
%%% \subsection{Astrophysically Meaningful Parameters}
%%% \label{subsec:piropars}
%%% \begin{tabular}{rl}
%%% Mass & BH of mass $M_{\mathrm{BH}} = 3-10 M_\odot$ \\
%%% Viscosity & Standard value of $\alpha = 0.1$ \cite{piro:2007}\\
%%% Geometrical parameter & $\eta = 0.3 -0.6$ \cite{piro:2007}\\
%%% Mass of the bound fragment & The approximate factor 0.2 in
%%% equation~\ref{eq:massfrag} can be varied from $0.2-0.5$.\\
%%% Starting Radius & Start at $r = 100 r_S$, where $r_S =
%%% \frac{GM}{c^2}$ is the gravitational radius. \\
%%% & End integration close to the ISCO. \\
%%% \end{tabular}
%%%
%%% \subsection{Usage of the python script}
%%%
%%% {\it Code usage:} {\tt \$ python pirogw.py}
%%%
%%% The code section below {\tt \#physical parameters} allows to specify the
%%% physical parameters
%%% of the system as defined in subsection~\ref{subsec:piropars}.
%%% The colatitude and azimuth of the system can be specified as well.
%%% The code
%%% section under {\tt \#parameters} allows to change the total run
%%% time (in seconds) of the integration (provided that the orbital radius
%%% is larger than $2.5 R_{\mathrm{ISCO}}$ where we take a multiple of the ISCO radius
%%% of a non-spinning black hole $R_{\mathrm{ISCO}} = 6 G
%%% M_{\mathrm{BH}}/c^2$ as a limit for the integration) as well as
%%% the sampling time {\tt dt} of the output.
%%%
%%% The script produces two output files:
%%% \begin{itemize}
%%% \item {\tt piro.dat} is a diagnosis and debug output file, containing
%%% the following variables: \\
%%% $\boxed{\mathrm{time} \qquad r\, \mathrm{(cm)} \qquad r\, (r_S) \qquad
%%% h_+ \qquad h_\times \qquad \Omega \qquad t_{GW} \qquad t_v }$
%%%
%%% \item {\tt piroM*eta*fac*.dat}, where {\tt *} denote the values of the physical parameters
%%% $M_{\mathrm{BH}}$, $\eta$ and the factor in the RHS of
%%% equation~\ref{eq:massfrag} given as input, is the production output file, containing:\\
%%% $\boxed{\mathrm{time} \qquad h_+ \qquad h_\times }$
%%%
%%% \end{itemize}
An example of the chirp-like signal obtained for the Piro \& Pfahl
model is shown in figure~\ref{fig:piro}.
\begin{figure}
\includegraphics[width=\linewidth]{pirowf.pdf}
\caption{Waveform computed following the Piro \& Pfahl model. The parameters of the system are
$M_{\mathrm{BH}}= 8 M_\odot$, $\eta = 0.3$, factor for the mass of the
bound fragment = 0.2. The strain corresponds to a face-on,
optimally-oriented system situated at 10 kpc.}
\label{fig:piro}
\end{figure}
\section{Changes to this Document}
\noindent {\bf T1100093-v2, update by C.~D.~Ott}
\begin{itemize}
\item Added more details on the quadrupole approximation used
in section 2.1.
\item Added a caveat regarding the van Putten waveforms.
\item Removed duplicate text at the end of section 3.
\item Updated most of the text in section 3 based on input from Tony Piro.
\item Improved description of how the gravitational wave signal
is calculated for the Piro \& Pfahl model.
\item Updated text and {\tt pirogw.py} script (now version 2) to
correct a slight inconsistency in calculating the orbital angular
velocity. Previously only the black hole mass was taken into
account. Now the total system mass $M_\mathrm{BH} + M_f$ is used.
\end{itemize}
\section*{Acknowledgements}
We thank Peter Kalmus for helpful discussions, Tony Piro for helpful
suggestions, and Eric Thrane and Sarah Gossan for checking the
calculations.
%\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Bibliography}
\label{s:bib}
\bibliographystyle{unsrtnat.cott}
\bibliography{disk_instabilities}
\end{document}