# Polarimeter optical spectrum analyzer

Eyal Buks  Andrew and Erna Viterbi Department of Electrical Engineering,
Technion, Haifa 32000 Israel

###### Abstract

A coherent optical spectrum analyzer is integrated with a rotating quarter
wave plate polarimeter. The combined polarimeter optical spectrum analyzer
(POSA) allows the extraction of the state of polarization with high spectral
resolution. POSA is used in this work to study two optical systems. The first
is an optical modulator based on a ferrimagnetic sphere resonator. POSA is
employed to explore the underlying magneto-optical mechanism responsible for
modulation sideband asymmetry. The second system under study is a cryogenic
fiber loop laser, which produces an unequally spaced optical comb.
Polarization measurements provide insights on the nonlinear processes
responsible for comb creation. Characterizations extracted from POSA data
provide guidelines for performance optimization of applications based on these
systems under study.

##  I Introduction

Optical state of polarization (SOP) can be measured using a variety of
techniques [1]. For some applications, the dependency of SOP on optical
wavelength $\lambda$ has to be determined. This dependency on $\lambda$ can be
obtained by combining a wavelength tunable band-pass optical filter with the
polarimeter being used to measure the SOP. This relatively simple
configuration has a spectral resolution that is limited by the ratio between
the filter’s free spectral range (FSR) and its finesse. Lowering the spectral
resolution below about $1\operatorname{GHz}$ over a wide tuning range in the
optical band is practically challenging. On the other hand, for some
applications, a much higher spectral resolution is needed.

A coherent optical spectrum analyzer (OSA) [2, 3, 4] is an instrument based on
heterodyne detection of the optical signal under study [5, 6, 7]. Commonly, a
coherent OSA has significantly higher spectral resolution compared to grating-
based OSA. In its standard configuration, SOP cannot be extracted from
coherent OSA data. In this work a fiber-based system is proposed, which
integrates a coherent OSA with a polarimeter (PM) that is based on a rotating
quarter wave plate (RQWP). The integrated polarimeter optical spectrum
analyzer (POSA) allows determining both SOP and degree of polarization (DOP)
with high spectral resolution.

The SOP can be described as a point in the Poincaré unit sphere [see Fig.
1(a)]. The colinear vertical, horizontal, diagonal and anti-diagonal SOP are
denoted by $\left|V\right\rangle$, $\left|H\right\rangle$,
$\left|D\right\rangle=2^{-1/2}\left(\left|H\right\rangle+\left|V\right\rangle\right)$
and
$\left|A\right\rangle=2^{-1/2}\left(\left|H\right\rangle-\left|V\right\rangle\right)$
, respectively, whereas the circular right-hand and left-hand SOP are denoted
by
$\left|R\right\rangle=2^{-1/2}\left(\left|H\right\rangle-i\left|V\right\rangle\right)$
and
$\left|L\right\rangle=2^{-1/2}\left(\left|H\right\rangle+i\left|V\right\rangle\right)$,
respectively. The unit vectors in the Poincaré sphere corresponding to the SOP
$\left|V\right\rangle$, $\left|H\right\rangle$, $\left|D\right\rangle$,
$\left|A\right\rangle$, $\left|R\right\rangle$ and $\left|L\right\rangle$, are
$\mathbf{\hat{x}}_{3}$, $-\mathbf{\hat{x}}_{3}$, $\mathbf{\hat{x}}_{1}$,
$-\mathbf{\hat{x}}_{1}$, $-\mathbf{\hat{x}}_{2}$ and $\mathbf{\hat{x}}_{2}$,
respectively.

In this work POSA is employed for studying two systems. The first is a
ferrimagnetic sphere resonator (FMSR) made of yttrium iron garnet (YIG) (see
sections III and IV). It has been recently demonstrated that optical single
sideband modulation (SSM) can be implemented in the telecom band using a FMSR
[8, 9]. Section V is devoted to POSA measurements of a cryogenic fiber loop
laser, which is operated in a region where an unequally spaced optical comb
(USOC) is formed [10]. For both systems, theoretical interpretation of POSA
measurement results is discussed.

##  II POSA

Figure 1: SOP. (a) The Poincaré sphere. (b) POSA setup.

POSA setup is schematically shown in Fig. 1(b). The input section contains two
polarization controllers [PC1 and PC2 in Fig. 1(b)], and a RQWP (which is
based on a rotating stage placed between two fiber collimators). Heterodyne
detection in the coherent OSA section is performed using a wavelength tunable
laser [TL2 in Fig. 1(b)], two polarized beam splitters [PBS1 and PBS2 in Fig.
1(b)], two 50:50 optical couplers [OC1 and OC2 in Fig. 1(b)], and two
differential photodetectors [DPD1 and DPD2 in Fig. 1(b)].

In some other RQWP-based polarimeters, the RQWP is directly attached to a PBS
and a photodetector [1]. In contrast, a single mode fiber is used in our setup
to connect the RQWP to the coherent OSA section [see Fig. 1(b)]. This inter-
section connection gives rise to a unitary transformation denoted by $B_{2}$.
The transformation $B_{2}$ can be manipulated using PC2, however, it is a
priori unknown. Consequently, the extraction process of SOP from POSA data,
which is explained below, is more complicated than the process that is
commonly employed in other RQWP-based polarimeters, for which $B_{2}$
represents the identity transformation [1].

Let $\rho$ be the density matrix of a given SOP, and let $P_{\mathrm{m}}$ be
the matrix representation of a given projection operator (associated with a
given polarization filter). The $2\times 2$ matrices $\rho$ and
$P_{\mathrm{m}}$ are expressed as
$\rho=\left(1/2\right)\left(1+\gamma\mathbf{\hat{n}}\cdot\bm{\sigma}\right)$
and
$P_{\mathrm{m}}=\left(1/2\right)\left(1+\mathbf{\hat{n}}_{\mathrm{m}}\cdot\bm{\sigma}\right)$,
where both
$\mathbf{\hat{n}}=\gamma^{-1}\left(P_{1},P_{2},P_{3}\right)=\left(\sin\theta_{\mathrm{n}}\cos\varphi_{\mathrm{n}},\cos\theta_{\mathrm{\
n}},\sin\theta_{\mathrm{n}}\sin\varphi_{\mathrm{n}}\right)$ and $\mathbf{\
\hat{n}}_{\mathrm{m}}=\left(m_{1},m_{2},m_{3}\right)=\left(\sin\theta\cos\varphi,\cos\theta,\sin\theta\sin\varphi\right)$
are real unit vectors (the over-hat symbol is used to denote unit vectors),
$0\leq\gamma\leq 1$ is the DOP, and
$\bm{\sigma}=\left(\sigma_{1},\sigma_{2},\sigma_{3}\right)$ is the Pauli
matrix vector

$\sigma_{1}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),\;\sigma_{2}=\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right),\;\sigma_{3}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)\;.$ (1)

Using the relation
$\left(\bm{\sigma}\cdot\mathbf{a}\right)\left(\bm{\sigma}\cdot\mathbf{b}\right)=\mathbf{a}\cdot\mathbf{b}+i\bm{\sigma}\cdot\left(\mathbf{a}\times\mathbf{b}\right)$,
one finds that the probability $p$ to find the SOP pointing in the
$\mathbf{\hat{n}}_{\mathrm{m}}$ direction is given by

$p=\operatorname{Tr}\left(\rho
P_{\mathrm{m}}\right)=\frac{1+\gamma\mathbf{\hat{n}}\cdot\mathbf{\hat{n}}_{\mathrm{m}}}{2}\;.$
(2)

While the unit vector $\mathbf{\hat{n}}_{\mathrm{m}}$ is associated with the
polarization filter of OSA trace 1 (DPD1), OSA trace 2 (DPD2) represents the
orthogonal SOP corresponding to the unit vector
$-\mathbf{\hat{n}}_{\mathrm{m}}$. For a given input signal having intensity
$I$, the OSA trace 1 (2) signal intensity $I_{+}$ ($I_{-}$) is given by
$I_{\pm}/I=\left(1/2\right)\left(1\pm\gamma\mathbf{\hat{n}}\cdot\mathbf{\hat{n}}_{\mathrm{m}}\right)$
[see Eq. (2)], and thus the following holds

$\gamma\mathbf{\hat{n}}\cdot\mathbf{\hat{n}}_{\mathrm{m}}=\frac{I_{+}-I_{-}}{I_{+}+I_{-}}\;.$
(3)

Any loss-less linear SOP transformation can be described using a unitary Jones
matrix $B\left(\mathbf{\hat{u}},\phi\right)$ given by [11]

$B\left(\mathbf{\hat{u}},\phi\right)\dot{=}\exp\left(-\frac{i\bm{\
\sigma}\cdot\mathbf{\hat{u}}\phi}{2}\right)=\mathbf{1}\cos\frac{\phi}{2}-i\bm{\sigma}\cdot\mathbf{\hat{u}}\sin\frac{\phi}{2}\;,$
(4)

where $\mathbf{\hat{u}}$ is a unit vector and $\phi$ is a rotation angle. By
expressing the unit vector $\mathbf{\hat{n}}$ as
$\mathbf{\hat{n}}=\mathbf{n}_{\parallel}+\mathbf{n}_{\perp}$, where
$\mathbf{n}_{\parallel}=\mathbf{\
\left(\hat{u}\cdot\mathbf{\hat{n}}\right)\hat{u}}$ (parallel component of
$\mathbf{\hat{n}}$ in the $\mathbf{\hat{u}}$ direction) and
$\mathbf{n}_{\perp}=\mathbf{\hat{u}}\times\left(\mathbf{\hat{n}}\times\mathbf{\hat{u}}\right)$
(perpendicular component), one finds that

$B\left(\mathbf{\hat{u}},\phi\right)\left(\bm{\sigma}\cdot\mathbf{\hat{n}}\right)B^{{\dagger}}\left(\mathbf{\hat{u}},\phi\right)=\bm{\
\sigma}\cdot\mathbf{\hat{n}}_{\mathrm{T}}\;,$ (5)

where

$\mathbf{\hat{n}}_{\mathrm{T}}=\mathbf{n}_{\parallel}+\mathbf{n}_{\perp}\cos\phi+\left(\mathbf{\hat{u}}\times\mathbf{n}_{\perp}\right)\sin\phi\;.$
(6)

The transformation from $\mathbf{\hat{n}}$ to $\mathbf{\hat{n}}_{\mathrm{T}}$
(6) is a rotation about the $\mathbf{\hat{u}}$ axis with angle $\phi$.

The unit vector $\mathbf{\hat{u}}$ corresponding to the RQWP is given by
$\mathbf{\hat{u}}=-\sin\left(2\alpha\right)\mathbf{\hat{x}}_{1}+\cos\left(2\alpha\right)\mathbf{\hat{x}}_{3}$,
where $\alpha$ is the RQWP axis angle, and the rotation angle $\phi$ is given
by $\phi=\pi/2$, thus [see Eqs. (5) and (6)]

$\displaystyle\gamma\mathbf{\hat{n}}_{\mathrm{T}}\cdot\mathbf{\hat{n}}_{\mathrm{m}}$
$\displaystyle=\gamma\left(\mathbf{\left(\hat{u}\cdot\mathbf{\hat{n}}\right)\hat{u}}+\left(\mathbf{\hat{u}}\times\mathbf{n}\right)\right)\cdot\mathbf{\hat{n}}_{\mathrm{m}}$
$\displaystyle=a_{0}+a_{1}\cos\left(2\alpha\right)+b_{1}\sin\left(2\alpha\right)$
$\displaystyle+a_{2}\cos\left(4\alpha\right)+b_{2}\sin\left(4\alpha\right)\;,$
(7)

where

$\displaystyle a_{0}$ $\displaystyle=\frac{P_{1}m_{1}+P_{3}m_{3}}{2}\;,$ (8)
$\displaystyle a_{1}$ $\displaystyle=P_{1}m_{2}-P_{2}m_{1}\;,$ (9)
$\displaystyle b_{1}$ $\displaystyle=P_{3}m_{2}-P_{2}m_{3}\;,$ (10)
$\displaystyle a_{2}$ $\displaystyle=\frac{P_{3}m_{3}-P_{1}m_{1}}{2}\;,$ (11)
$\displaystyle b_{2}$ $\displaystyle=-\frac{P_{1}m_{3}+P_{3}m_{1}}{2}\;.$ (12)

The following holds [see Eqs. (11) and (12)]

$a_{2}^{2}+b_{2}^{2}=\frac{\gamma^{2}\sin^{2}\theta\sin^{2}\theta_{\mathrm{n}}}{4}\;.$
(13)

The Poincaré vector $\mathbf{P}=\left(P_{1},P_{2},P_{3}\right)$ can be
extracted from Eqs. (8)-(12), and the measured values of $a_{0}$, $a_{1}$,
$b_{1}$, $a_{2}$ and $b_{2}$, provided that the transformation $B_{2}$, which
determines the unit vector $\mathbf{\
\hat{n}}_{\mathrm{m}}=\left(m_{1},m_{2},m_{3}\right)$, is given. POSA
calibration (which is needed because $B_{2}$ is a priori unknown) is done by
varying the angles $\theta_{\mathrm{n}}$ and $\varphi_{\mathrm{n}}$ associated
with the SOP transformation $B_{1}$ controlled by PC1, while keeping $\theta$
and $\varphi$ unchanged ($\theta$ and $\varphi$ are determined by the SOP
transformation $B_{2}$ from the RQWP to the coherent OSA section, which can be
manipulated using PC2). This process allows determining the angle
$\theta_{\mathrm{n}}$ using the relation [see Eq. (13)]

$\sin^{2}\theta_{\mathrm{n}}=\frac{a_{2}^{2}+b_{2}^{2}}{\max_{\theta_{\mathrm{n}},\varphi_{\mathrm{n}}}\left(a_{2}^{2}+b_{2}^{2}\right)}\;.$
(14)

The unit vector
$\mathbf{\hat{n}}_{\mathrm{m}}=\left(\sin\theta\cos\varphi,\cos\theta,\sin\theta\sin\varphi\right)$
can be manipulated by PC2 [see Fig. 1(b)]. As can be seen from Eq. (13), the
angle $\theta$ can be tuned such that $\sin^{2}\theta=1$, by maximizing
$a_{2}^{2}+b_{2}^{2}$ using PC2. For this case
$a_{1}^{2}+b_{1}^{2}+4\left(a_{2}^{2}+b_{2}^{2}\right)=\gamma^{2}$. This PC2
calibration, which greatly simplifies SOP extraction from POSA data, has been
performed prior to all measurements with both a driven FMSR (sections III and
IV) and with a cryogenic fiber loop laser (section V).

##  III FMSR modulator

FMSR are widely employed as magnetically tunable microwave filters having high
quality factor. The FMSR ellipsoidal shape allows uniform magnetization.
Magneto-optical (MO) coupling [12, 13, 14, 15, 16, 17, 18, 19, 20] between
FMSR optical [21, 22, 23] and magnetic Walker [24] modes can be used for
implementing optical modulation [8, 25, 26, 27, 28, 29, 30, 31, 32, 9, 33, 34,
35]. The sequence of sidebands that are generated in a FMSR-based modulator
can be determined from Brillouin scattering selection rules [8, 25, 26, 27,
29, 30, 31, 36, 37, 38, 39, 40], and angular momentum conservation in photon-
magnon scattering [41, 42, 43, 44, 45, 46, 47]. Contrary to some other
modulation techniques (such as amplitude, phase and frequency modulations)
[48, 49, 50], symmetry between Stokes and anti-Stokes sidebands can be broken
in FMSR-based modulation [8, 9, 35]. In particular, the method of single
sideband modulation (SSM), which allows reducing both transmission power and
bandwidth, can be implemented using FMSR modulation.

Figure 2: Experimental setup. Optical components and fibers are red colored,
whereas blue color is used to label MW components and coaxial cables. (a)
Optical fibers are installed on both sides of the FMSR for transmission of
light through the sphere (see Ref. [35] for more details on the experimental
setup). (b) The cryogenic fiber loop laser (see Ref. [51] for more details).

The FMSR modulator is schematically shown in Fig. 2(a). An FMSR made of YIG
having radius of $R_{\mathrm{s}}=125\operatorname{\mu m}$ is held by two
ceramic ferrules (CF). The two CFs, which are held by a concentric sleeve,
provide transverse alignment for both input and output single mode optical
fibers. All optical measurements are performed in the telecom band, in which
YIG has refractive index of $2.19$ and absorption coefficient of
$\left(0.5\mathrm{m}\right)^{-1}$ [52]. The intensity and SOP of light
illuminating the FMSR are controlled by an optical attenuator (Att) and a PC,
respectively. The blue-colored microwave (MW) components shown in Fig. 2(a)
allow both driving and detection of FMSR magnetic resonances. The FMSR is
inductively coupled to a microwave loop antenna (MWA). Magnetic resonances are
identified using a vector network analyzer (VNA). A signal generator (SG) is
employed for driving, and the response is monitored using a radio frequency
spectrum analyzer (RFSA). A circulator (C) and a splitter (S) are employed to
direct the input and output microwave signals [see Fig. 2(a)].

The angular frequency of the FMSR Kittel mode $\omega_{\mathrm{m}}$ is
approximately given by
$\omega_{\mathrm{m}}=\mu_{0}\gamma_{\mathrm{e}}H_{\mathrm{s}}$, where
$H_{\mathrm{s}}$ is the static magnetic field, $\mu_{0}$ is the free space
permeability, and $\gamma_{\mathrm{e}}/2\pi=28\mathrm{GHz}\mathrm{T}^{-1}$ is
the gyromagnetic ratio [53, 54, 20, 55]. The applied static magnetic field
$\mathbf{H}_{\mathrm{s}}$ is controlled by adjusting the relative position of
a magnetized Neodymium using a motorized stage. The static magnetic field
$\mathbf{H}_{\mathrm{s}}$ is normal to the unit vector $\mathbf{\hat{k}}$
pointing in the light propagation direction, and the MWA driving magnetic
field is nearly parallel to $\mathbf{\hat{k}}$. The FMSR is installed such
that its crystallographic direction $\left[111\right]$ is aligned to
$\mathbf{\hat{k}}$.

The plot in Fig. 3(2) shows POSA DPD1 signal intensity [see Fig. 1(b)] as a
function of TL2 optical wavelength $\lambda$. For this plot the TL1 [see Fig.
2(a)] power is $P_{\mathrm{L}}=6\operatorname{mW}$ and wavelength is
$\lambda_{\mathrm{L}}=1537.7\operatorname{nm}$. The MWA driving, which has
power of $12\operatorname{dBm}$, is frequency tuned to resonance at
$\omega_{\mathrm{m}}/\left(2\pi\right)=4.1755\operatorname{GHz}$. The driving-
induced sidebands around the central peak [see Fig. 3(2)] have wavelengths
$\lambda_{\mathrm{L}}\left(1\pm\lambda_{\mathrm{L}}\omega_{\mathrm{m}}/\left(2\pi
c\right)\right)=\lambda_{\mathrm{L}}\pm 32.9\operatorname{pm}$, where $c$ is
the vacuum speed of light. Sidebands’ intensities strongly depend on the input
light SOP. This dependency is explored in the next section, which is devoted
to the FMSR transverse dielectric tensor.

Figure 3: FMSR. Normalized sideband intensities
$I_{\mathrm{L}}/I_{\mathrm{P}}$ and $I_{\mathrm{R}}/I_{\mathrm{P}}$ are shown
in (1) and (3), respectively ($I_{\mathrm{P}}$ is the central peak intensity).
The letters a, b and c in the Poincaré plots’ labels refer to values of
$\upsilon$ indicated in (1) and (3). The numbers 1, 2 and 3 in the Poincaré
plots’ labels refer to the left sideband, central peak and right sideband,
respectively. The Poincaré vectors are red colored.

##  IV FMSR transverse dielectric tensor

In a Cartesian coordinate system, for which the $[100]$, $[010]$, and $[001]$
crystallographic directions are pointing in the $x$, $y$ and $z$ axes,
respectively, the dielectric tensor is represented by a $3\times 3$ matrix
denoted by $\varepsilon$ [56, 57, 58]. The Onsager reciprocal relations, which
originate from time-reversal symmetry, read
$\epsilon_{nm}\left(\mathbf{M}\right)=\epsilon_{mn}\left(-\mathbf{M}\right)$,
where $\mathbf{M}=\left(M_{1},M_{2},M_{3}\right)\mathbf{\ }$is the
magnetization vector, and $n,m\in\left\\{1,2,3\right\\}$. The MO
Stoner–Wohlfarth energy density is given by
$u=\left(1/4\right)\mathbf{E}^{\ast}\varepsilon\mathbf{E}^{\mathrm{T}}$, where
$\mathbf{E}=\left(E_{1},E_{2},E_{3}\right)$ is the electric field vector
phasor.

To second order in the magnetization $\mathbf{M}$, the dielectric tensor can
be expressed as
$\varepsilon=\varepsilon_{\mathrm{r}}+\varepsilon_{1}+\varepsilon_{2}$, where
$\varepsilon_{\mathrm{r}}$ is the relative permittivity. The first order (in
$\mathbf{M}$) contribution to the dielectric tensor is given by
$\varepsilon_{1}=iQ_{\mathrm{s}}\epsilon_{ijk}M_{k}$, where $\epsilon_{ijk}$
is the Levi-Civita symbol. For YIG in saturated magnetization
$Q_{\mathrm{s}}\simeq 10^{-4}$ [59]. Crystal cubic symmetry implies that the
second order contribution $\varepsilon_{2}$ can be expressed as
$\varepsilon_{2}=g\operatorname{diag}\left(M_{1}^{2},M_{2}^{2},M_{3}^{2}\right)+g_{12}\mathbf{M}\mathbf{M}^{\mathrm{T}}+2g_{44}\mathbf{M}^{\mathrm{T}}\mathbf{M}$,
where $g$, $g_{12}$, and $g_{44}$ are constants [60, 61] (in the notation used
in Ref. [60] $g=g_{11}-g_{12}-2g_{44}$).

Consider the coordinate transformation
$\mathbf{M}^{\prime}=\mathbf{M}\mathcal{R}$, where $\mathcal{R}$ is a $3\times
3$ rotation matrix. The energy density can be expressed as
$u=\left(1/4\right)\mathbf{E}^{\prime\ast}\varepsilon^{\prime}\mathbf{E}^{\prime\mathrm{T}}$,
where $\mathbf{E}^{\prime}=\mathbf{E}\mathcal{R}$ is the transformed value of
the electric field $\mathbf{E}$, and where
$\varepsilon^{\prime}=\mathcal{R}^{-1}\varepsilon\mathcal{R}$ is the
transformed value of the dielectric tensor $\varepsilon$ (note that
$\mathcal{R}^{\mathrm{T}}=\mathcal{R}^{-1}$ and that $\mathcal{R}$ is real).
The corresponding transformed first (second) order contribution to the
dielectric tensor is denoted by $\varepsilon_{1}^{\prime}$
($\varepsilon_{2}^{\prime}$).

The following holds
$\mathbf{E}^{\ast}\varepsilon_{1}\mathbf{E}^{\mathrm{T}}=-iQ_{\mathrm{s}}\mathbf{E}^{\ast}\cdot\left(\mathbf{M}\times\mathbf{E}\right)$,
thus using the matrix identity
$\left(\mathcal{R}\mathbf{X}\right)\times\left(\mathcal{R}\mathbf{Y}\right)=$
$\mathcal{R}\left(\mathbf{X}\times\mathbf{Y}\right)$, where $\mathbf{X}$ and
$\mathbf{Y}$ are vectors and $\mathcal{R}$ is a rotation matrix, one finds
that
$\varepsilon_{1}^{\prime}\left(\mathbf{M}^{\prime}\right)=\varepsilon_{1}\left(\mathbf{M}\right)$,
i.e. $\varepsilon_{1}^{\prime}=iQ_{\mathrm{s}}\epsilon_{ijk}M_{k}^{\prime}$.
The transformed second order contribution $\varepsilon_{2}^{\prime}$ is given
by
$\varepsilon_{2}^{\prime}=gE_{\mathrm{g}}^{\prime}+g_{12}\mathbf{M}^{\prime}\left(\mathbf{M}^{\prime}\right)^{\mathrm{T}}+2g_{44}\left(\mathbf{M}^{\prime}\right)^{\mathrm{T}}\mathbf{M}^{\prime}$,
where

$E_{\mathrm{g}}^{\prime}=\sum_{n=1}^{3}\mathcal{R}^{-1}\eta_{n,n}\mathcal{R}\left(\mathbf{M}^{\prime}\right)^{\mathrm{T}}\mathbf{M}^{\prime}\mathcal{R}^{-1}\eta_{n,n}\mathcal{R}\;,$
(15)

and where $\eta_{n,m}$ is a $3\times 3$ matrix, whose entries are given by
$\left(\eta_{n,m}\right)_{n^{\prime},m^{\prime}}=\delta_{n,n^{\prime}}\delta_{m,m^{\prime}}$.
Note that only the term $E_{\mathrm{g}}^{\prime}$ given by Eq. (15) gives rise
to dependency of $\varepsilon_{2}^{\prime}$ on crystallographic directions.

The rotation matrix $\mathcal{R}$, which is given by

$\mathcal{R}=\left(\begin{array}[]{ccc}\frac{\cos\phi}{\sqrt{2}}-\frac{\sin\phi}{\sqrt{6}}&-\frac{\cos\phi}{\sqrt{6}}-\frac{\sin\phi}{\sqrt{2}}&\frac{1}{\sqrt{3}}\\\
\frac{2\sin\phi}{\sqrt{6}}&\frac{2\cos\phi}{\sqrt{6}}&\frac{1}{\sqrt{3}}\\\
-\frac{\cos\phi}{\sqrt{2}}-\frac{\sin\phi}{\sqrt{6}}&-\frac{\cos\phi}{\sqrt{6}}+\frac{\sin\phi}{\sqrt{2}}&\frac{1}{\sqrt{3}}\end{array}\right)\;,$
(16)

maps the $z$ direction to $\left[111\right]$ (which is parallel to light
propagation direction $\mathbf{\hat{k}}$), and the $xy$ plane to
$\left(111\right)$, with a variable angle $\phi$. Consider the case where the
static magnetic field is applied parallel to the $x^{\prime}$ direction. It is
assumed that $M_{1}^{\prime}\simeq M_{0}$,
$M_{2}^{\prime}/M_{0}\equiv\mu_{2}\ll 1$ and
$M_{3}^{\prime}/M_{0}\equiv\mu_{3}\ll 1$, where $M_{0}$ is the saturation
magnetization.

In the truncation approximation, the transverse dielectric tensor is taken to
be the $2\times 2$ upper left block of the $3\times 3$ tensor
$\varepsilon^{\prime}$. Moreover, terms independent on $Q_{\mathrm{s}}$,
$\mu_{2}$ and $\mu_{3}$ are disregarded. In this approximation
$\varepsilon^{\prime}-\varepsilon_{\mathrm{r}}$ is represented by a $2\times
2$ matrix given by
$M_{0}^{2}\left(\mu_{3}\mathbf{a}^{\prime}+\mu_{2}\mathbf{a}^{\prime\prime}\right)\cdot\bm{\sigma}$,
where $\mathbf{a}^{\prime}$ ($\mathbf{a}^{\prime\prime}$) is the real
(imaginary) part of the vector
$\mathbf{a}=\mathbf{a}^{\prime}+i\mathbf{a}^{\prime\prime}$, which is given by
(to first order in $Q_{\mathrm{s}}$)

$\mathbf{a}=-\frac{\sqrt{2}g}{3}\left(\cos
3\phi-\frac{3i\left(\frac{1}{3}+\frac{2g_{44}}{g}\right)}{\sqrt{2}},\frac{3Q_{\mathrm{s}}}{\sqrt{2}M_{0}g},\sin
3\phi\right)\;,$ (17)

and where $\bm{\sigma}$ is the Pauli matrix vector [see Eq. (1)].
Alternatively

$\varepsilon^{\prime}-\varepsilon_{\mathrm{r}}=\frac{M_{0}^{2}}{2}\left(\mu_{+}\mathbf{a}^{\ast}+\mu_{-}\mathbf{a}\right)\cdot\bm{\sigma}\;,$
(18)

where $\mu_{\pm}=\mu_{3}\pm i\mu_{2}$.

The Hermitian matrices
$\left(\mathbf{a}^{\prime\prime}\times\mathbf{a}^{\prime}\right)\cdot\bm{\sigma}$,
$S^{{\dagger}}S$ and $SS^{{\dagger}}$, where $S=\mathbf{a}\cdot\bm{\sigma}$
[see Eq. (18)], share the same eigenvectors, which are denoted by
$\left|\mathbf{a}_{+}\right\rangle$ and $\left|\mathbf{a}_{-}\right\rangle$.
The transformed matrix $\tilde{S}=U^{-1}SU$, where $U$ is the $2\times 2$
unitary transformation matrix corresponding to the basis
$\left\\{\left|\mathbf{a}_{+}\right\rangle,\left|\mathbf{a}_{-}\right\rangle\right\\}$,
has a hollow form [see Eq. (6.150) of Ref. [62]]. The matrix representation of
$\varepsilon^{\prime}-\varepsilon_{\mathrm{r}}$ in that basis, which is
denoted by $\varepsilon_{\mathrm{h}}$, is given by

$\varepsilon_{\mathrm{h}}=\frac{M_{0}^{2}}{2}\left(\begin{array}[]{cc}0&\mu_{-}\zeta_{1}+\mu_{+}\zeta_{2}^{\ast}\\\
\mu_{-}\zeta_{2}+\mu_{+}\zeta_{1}^{\ast}&0\end{array}\right)\;,$ (19)

where $\zeta_{1}\zeta_{2}=\mathbf{a}\cdot\mathbf{a}$,
$\left|\zeta_{1}\right|^{2}+\left|\zeta_{2}\right|^{2}=2\mathbf{a}\cdot\mathbf{a}^{\ast}$
and
$\left|\zeta_{2}\right|^{2}-\left|\zeta_{1}\right|^{2}=4\left|\mathbf{a}^{\prime\prime}\times\mathbf{a}^{\prime}\right|$.

Both scattering selection rules and sideband asymmetry (between Stokes and
anti-Stokes) can be characterized in terms of the ratio $\zeta_{\mathbf{a}}$,
which is defined by [see Eq. (19), and note that
$\left|\zeta_{\mathbf{a}}\right|\leq 1$]

$\zeta_{\mathbf{a}}=\frac{\left|\zeta_{2}\right|^{2}-\left|\zeta_{1}\right|^{2}}{\left|\zeta_{1}\right|^{2}+\left|\zeta_{2}\right|^{2}}=\frac{2\left|\mathbf{a}^{\prime\prime}\times\mathbf{a}^{\prime}\right|}{\left|\mathbf{a}^{\prime}\right|^{2}+\left|\mathbf{a}^{\prime\prime}\right|^{2}}\;.$
(20)

Note that sideband asymmetry is possible only when
$\left|\mathbf{a}^{\prime\prime}\times\mathbf{a}^{\prime}\right|>0$ [see Eqs.
(19) and (20)]. For radio and optical communication, it is well known that
sideband asymmetry is impossible when only the amplitude, or only the
frequency, or only the phase, is modulated by a monochromatic signal.
Moreover,with simultaneously applied two different modulation methods,
sideband asymmetry becomes possible only when one modulation method is phase-
shifted with respect to the other. These properties are all reflected by the
condition $\left|\mathbf{a}^{\prime\prime}\times\mathbf{a}^{\prime}\right|>0$.

For YIG $\left(Q_{\mathrm{s}}/\left(M_{0}g\right)\right)^{2}=44.2$ and
$2g_{44}/g=-4.0$ [20], hence $\zeta_{\mathbf{a}}\simeq 0.84$ [see Eqs. (17)
and (20), and note that the term $\sin^{2}3\phi\ $has a relatively small
effect]. The largest sideband asymmetry is obtained for input states
$\left|\mathbf{a}_{+}\right\rangle$ and $\left|\mathbf{a}_{-}\right\rangle$
[eigenvectors of
$\left(\mathbf{a}^{\prime\prime}\times\mathbf{a}^{\prime}\right)\cdot\bm{\sigma}$].
As can be seen from Eq. (17), for YIG the vector
$\mathbf{a}^{\prime\prime}\times\mathbf{a}^{\prime}$ is nearly parallel to the
unit vector $\mathbf{\hat{z}}$.

Intensities of both left (shorter wavelength, anti-Stokes) and right (longer
wavelength, Stokes) sidebands can be manipulated using the PC shown in Fig.
2(a), which controls the input SOP. The vector of three voltages controlling
this PC are denoted by
${\mathbf{V}}_{\mathrm{PC}}=\left(V_{\mathrm{PC},1},V_{\mathrm{PC},2},V_{\mathrm{PC},3}\right)$.
A searching procedure is implemented to determine the value
${\mathbf{V}}_{\mathrm{PC,L}}$ (${\mathbf{V}}_{\mathrm{PC,R}}$), for which the
intensity of the left (right) sideband, which is denoted by $I_{\mathrm{L}}$
($I_{\mathrm{R}}$), is maximized. The measurements shown in Fig. 3 are
performed with values of ${\mathbf{V}}_{\mathrm{PC}}$ given by
${\mathbf{V}}_{\mathrm{PC}}=\upsilon{\mathbf{V}}_{\mathrm{PC,L}}+\left(1-\upsilon\right){\mathbf{V}}_{\mathrm{PC,R}}$,
where $0\leq\upsilon\leq 1$. The Poincaré vector
$\mathbf{P}=\left(P_{1},P_{2},P_{3}\right)$ can be extracted from POSA data
using the method above-discussed in section II. Poincaré vectors $\mathbf{P}$
are shown in Fig. 3 for three different values of the parameter $\upsilon$,
which are indicated in Fig. 3(1) and (3) by the letters a, b and c. The
numbers 1, 2 and 3 in the Poincaré plots’ labels in Fig. 3 refer to the left
sideband, central peak and right sideband, respectively.

A value for the asymmetry parameter given by $\zeta_{\mathbf{a}}=0.92$ is
extracted from the data shown in Fig. 3(1) and (3). This result suggests that
the above-discussed theoretical analysis underestimates $\zeta_{\mathbf{a}}$
(recall that the calculated value is $0.84$) [44]. Disagreement in SOP between
theory and experiment is quantified by the parameter
$\delta_{\mathrm{P}}=\left\|\mathbf{P}_{\mathrm{E}}-\mathbf{P}_{\mathrm{T}}\right\|/\left\|\mathbf{P}_{\mathrm{T}}\right\|$,
where $\mathbf{P}_{\mathrm{E}}$ denotes a measured Poincaré vector, and
$\mathbf{P}_{\mathrm{T}}$ is theoretically derived using Eq. (19). For the
data shown in Fig. 3 $\delta_{\mathrm{P}}\lesssim 0.13$. The disagreement is
partially attributed to magneto-dichroism (i.e. polarization-dependent
absorption), which is theoretically disregarded. Further study is needed to
explore the effect of magneto-dichroism and other possible mechanisms
contributing to the observed disagreement.

##  V USOC

USOC formation in a cryogenic fiber loop laser has been recently demonstrated
[10, 51]. USOC can be employed as a tunable multi-mode lasing source [63].
Multimode lasing has a variety of applications in the fields of sensing,
spectroscopy, signal processing and communication. Multimode lasing in the
telecom band has been demonstrated by integrating Erbium doped fibers (EDF)
cooled by liquid nitrogen into a fiber ring laser [64, 65, 66, 67]. It has
been recently proposed that EDF operating at low temperatures can be used for
storing quantum information [68, 69, 70, 71, 72]. For some USOC-based
applications, DOP specification is important. As is discussed below, the POSA
high spectral resolution allows such a specification for individual USOC
peaks.

A sketch of the cryogenic fiber loop laser used for studying USOC is shown in
Fig. 2(b). The loop is made of an undoped single mode fiber (Corning 28), and
a $5\operatorname{m}$ long EDF, which is used for optical amplification. The
EDF has absorption of $30$ dB $\mathrm{m}^{-1}$ at $1530\mathrm{nm}$, and mode
field diameter of $6.5\mu\mathrm{m}$ at $1550\mathrm{nm}$. The cold EDF is
integrated with a room temperature fiber loop using a wavelength-division
multiplexing (WDM) device. The EDF is pumped using a $980\operatorname{nm}$
laser diode (LD) biased with current denoted by $I_{\mathrm{D}}$. A 10:90 OC,
and two isolators [labeled by arrows in the sketch shown in Fig. 2(b)], are
integrated into the fiber loop. The output port of the 10:90 OC is connected
to the POSA setup shown in Fig. 1(b).

The plot in Fig. 4(a) shows the total POSA signal intensity $I_{+}+I_{-}$ as a
function of the TL2 wavelength $\lambda$ in the band where USOC is formed
($1539.7-1541.5\operatorname{nm}$). The underlying mechanism responsible for
USOC formation has remained mainly unknown [10, 51]. POSA is employed to
determine the DOP of each USOC peak. The result, which is shown in Fig. 4(b),
reveals that the first (shortest wavelength) USOC peak has the largest DOP.
Note that commonly in fiber lasers DOP is below about $0.15$, unless a
polarization-maintaining fiber is used [73]. Further study is needed to
explore how higher USOC DOP, which is needed for some applications, can be
achieved.

Figure 4: USOC intensity (a) and DOP (b). Diode current is
$0.2\operatorname{A}$, and EDF temperature is $3.2\operatorname{K}$.

##  VI Discussion and summary

A variety of POSA configurations have been developed and optimized for
specific applications [74]. Commonly, wavelength separation in POSA
instruments is achieved by integrating either a diffraction grating or a
tunable optical filter. For some applications, in fields including astronomy
and mineralogy analysis [75], these methods provide sufficiently high spectral
resolution. However, a much higher spectral resolution can be achieved by
implementing the method of coherent heterodyne detection [2, 3, 4].

Our coherent POSA setup operates in the entire C band
($1530-1565\operatorname{nm}$), and it has dynamic range of about 70 dB. It
allows SOP measurement with spectral resolution of $5\operatorname{MHz}$
(corresponding wavelength resolution of $0.04\operatorname{pm}$ in the telecom
band). SOP can be alternatively measured by combining a polarimeter with a
tunable optical filter. The filter-based method is simpler to implement,
however, its spectral resolution is limited by the ratio between the filter’s
FSR and finesse. Filter-based SOP measurements of a driven FMSR have been
reported in [44, 8]. In that study, the FMSR resonance frequency
$\omega_{\mathrm{m}}/\left(2\pi\right)$ was tuned close to
$7\operatorname{GHz}$, and a scanning Fabry–Pérot etalon having FSR of
$10\operatorname{GHz}$ (corresponding wavelength spacing of
$80\operatorname{pm}$ in the telecom band) was employed as a tunable filter
[44, 8]. Accurate SOP measurements of a driven FMSR can be obtained provided
that the filter is capable of resolving individual sidebands. This task is
challenging for a driven FMSR, because the power carried by both Stokes and
anti-Stokes optical sidebands is commonly at least 3 orders of magnitude
smaller than the power carried by the pump tone [see Fig. 3(2)]. Thus, for SOP
measurements of a driven FMSR, the much higher spectral resolution offered by
the coherent POSA significantly improves accuracy. For the USOC, as can be
seen from Fig. 4(a), the spacing between peaks is far too small to practically
allow any accurate SOP measurements using the filter-based method.

Future work will be devoted to optimize some applications, that are based on
the optical systems under-study. For the FMSR modulator, ways to further
increase the asymmetry parameter $\zeta_{\mathbf{a}}$ will be explored to
allow more efficient optical communication based on SSM. Moreover, ways to
enhance USOC DOP will be investigated in order to open the way for some novel
applications that require a multi-mode lasing source having high level of
intermode coherency.

## References

  * [1] Dennis H Goldstein,  Polarized light,  CRC press, 2017.
  * [2] X Steve Yao, Bo Zhang, Xiaojun Chen, and Alan E Willner,  “Real-time optical spectrum analysis of a light source using a polarimeter”,  Optics Express, vol. 16, no. 22, pp. 17854–17863, 2008.
  * [3] Eunha Kim, Digant Dave, and Thomas E Milner,  “Fiber-optic spectral polarimeter using a broadband swept laser source”,  Optics communications, vol. 249, no. 1-3, pp. 351–356, 2005.
  * [4] HK Kim, SK Kim, HG Park, and Byoung Yoon Kim,  “Polarimetric fiber laser sensors”,  Optics letters, vol. 18, no. 4, pp. 317–319, 1993.
  * [5] Douglas M Baney, Bogdan Szafraniec, and Ali Motamedi,  “Coherent optical spectrum analyzer”,  Ieee Photonics Technology Letters, vol. 14, no. 3, pp. 355–357, 2002.
  * [6] Kunpeng Feng, Jiwen Cui, Hong Dang, Weidong Wu, Xun Sun, Xuelin Jiang, and Jiubin Tan,  “An optoelectronic equivalent narrowband filter for high resolution optical spectrum analysis”,  Sensors, vol. 17, no. 2, pp. 348, 2017.
  * [7] Hong Dang, Huanhuan Liu, Linqi Cheng, Yuqi Tian, Jinna Chen, Kunpeng Feng, Jiwen Cui, and Perry Ping Shum,  “Deconvolutional suppression of resolution degradation in coherent optical spectrum analyzer”,  Journal of Lightwave Technology, vol. 13, pp. 4430, 2023.
  * [8] JA Haigh, Andreas Nunnenkamp, AJ Ramsay, and AJ Ferguson,  “Triple-resonant brillouin light scattering in magneto-optical cavities”,  Physical review letters, vol. 117, no. 13, pp. 133602, 2016.
  * [9] Cheng-Zhe Chai, Zhen Shen, Yan-Lei Zhang, Hao-Qi Zhao, Guang-Can Guo, Chang-Ling Zou, and Chun-Hua Dong,  “Single-sideband microwave-to-optical conversion in high-q ferrimagnetic microspheres”,  Photonics Research, vol. 10, no. 3, pp. 820–827, 2022.
  * [10] Eyal Buks,  “Low temperature spectrum of a fiber loop laser”,  Physics Letters A, vol. 458, pp. 128591, 2023.
  * [11] Richard J Potton,  “Reciprocity in optics”,  Reports on Progress in Physics, vol. 67, no. 5, pp. 717, 2004.
  * [12] Babak Zare Rameshti, Silvia Viola Kusminskiy, James A Haigh, Koji Usami, Dany Lachance-Quirion, Yasunobu Nakamura, Can-Ming Hu, Hong X Tang, Gerrit EW Bauer, and Yaroslav M Blanter,  “Cavity magnonics”,  Physics Reports, vol. 979, pp. 1–61, 2022.
  * [13] Silvia Viola Kusminskiy,  “Cavity optomagnonics”,  in Optomagnonic Structures: Novel Architectures for Simultaneous Control of Light and Spin Waves, pp. 299–353. World Scientific, 2021.
  * [14] Na Zhu, Xufeng Zhang, Xu Han, Chang-Ling Zou, and Hong X Tang,  “Inverse faraday effect in an optomagnonic waveguide”,  arXiv:2012.11119, 2020.
  * [15] Dominik M Juraschek, Derek S Wang, and Prineha Narang,  “Sum-frequency excitation of coherent magnons”,  Physical Review B, vol. 103, no. 9, pp. 094407, 2021.
  * [16] VASV Bittencourt, I Liberal, and S Viola Kusminskiy,  “Light propagation and magnon-photon coupling in optically dispersive magnetic media”,  Physical Review B, vol. 105, no. 1, pp. 014409, 2022.
  * [17] Xufeng Zhang, Na Zhu, Chang-Ling Zou, and Hong X Tang,  “Optomagnonic whispering gallery microresonators”,  Physical review letters, vol. 117, no. 12, pp. 123605, 2016.
  * [18] JP Van der Ziel, Peter S Pershan, and LD Malmstrom,  “Optically-induced magnetization resulting from the inverse faraday effect”,  Physical review letters, vol. 15, no. 5, pp. 190, 1965.
  * [19] Shasha Zheng, Zhenyu Wang, Yipu Wang, Fengxiao Sun, Qiongyi He, Peng Yan, and HY Yuan,  “Tutorial: Nonlinear magnonics”,  arXiv:2303.16313, 2023.
  * [20] Daniel D Stancil and Anil Prabhakar,  Spin waves,  Springer, 2009.
  * [21] Anatolii N Oraevsky,  “Whispering-gallery waves”,  Quantum electronics, vol. 32, no. 5, pp. 377, 2002.
  * [22] Gerhard Schunk, Josef U Fürst, Michael Förtsch, Dmitry V Strekalov, Ulrich Vogl, Florian Sedlmeir, Harald GL Schwefel, Gerd Leuchs, and Christoph Marquardt,  “Identifying modes of large whispering-gallery mode resonators from the spectrum and emission pattern”,  Optics express, vol. 22, no. 25, pp. 30795–30806, 2014.
  * [23] Michael L Gorodetsky and Aleksey E Fomin,  “Geometrical theory of whispering-gallery modes”,  IEEE Journal of Selected Topics in Quantum Electronics, vol. 12, no. 1, pp. 33–39, 2006.
  * [24] Laurence R Walker,  “Magnetostatic modes in ferromagnetic resonance”,  Physical Review, vol. 105, no. 2, pp. 390, 1957.
  * [25] A Osada, A Gloppe, Y Nakamura, and K Usami,  “Orbital angular momentum conservation in brillouin light scattering within a ferromagnetic sphere”,  New Journal of Physics, vol. 20, no. 10, pp. 103018, 2018.
  * [26] A Osada, R Hisatomi, A Noguchi, Y Tabuchi, R Yamazaki, K Usami, M Sadgrove, R Yalla, M Nomura, and Y Nakamura,  “Cavity optomagnonics with spin-orbit coupled photons”,  Physical review letters, vol. 116, no. 22, pp. 223601, 2016.
  * [27] Sanchar Sharma, Yaroslav M Blanter, and Gerrit EW Bauer,  “Light scattering by magnons in whispering gallery mode cavities”,  Physical Review B, vol. 96, no. 9, pp. 094412, 2017.
  * [28] Sanchar Sharma, Yaroslav M Blanter, and Gerrit EW Bauer,  “Optical cooling of magnons”,  Physical review letters, vol. 121, no. 8, pp. 087205, 2018.
  * [29] Evangelos Almpanis,  “Dielectric magnetic microparticles as photomagnonic cavities: Enhancing the modulation of near-infrared light by spin waves”,  Physical Review B, vol. 97, no. 18, pp. 184406, 2018.
  * [30] R Zivieri, P Vavassori, L Giovannini, F Nizzoli, Eric E Fullerton, M Grimsditch, and V Metlushko,  “Stokes–anti-stokes brillouin intensity asymmetry of spin-wave modes in ferromagnetic films and multilayers”,  Physical Review B, vol. 65, no. 16, pp. 165406, 2002.
  * [31] BERNARD Desormiere and HENRI Le Gall,  “Interaction studies of a laser light with spin waves and magnetoelastic waves propagating in a yig bar”,  IEEE Transactions on Magnetics, vol. 8, no. 3, pp. 379–381, 1972\.
  * [32] Zeng-Xing Liu, Bao Wang, Hao Xiong, and Ying Wu,  “Magnon-induced high-order sideband generation”,  Optics Letters, vol. 43, no. 15, pp. 3698–3701, 2018.
  * [33] Na Zhu, Xufeng Zhang, Xu Han, Chang-Ling Zou, Changchun Zhong, Chiao-Hsuan Wang, Liang Jiang, and Hong X Tang,  “Waveguide cavity optomagnonics for microwave-to-optics conversion”,  Optica, vol. 7, no. 10, pp. 1291–1297, 2020.
  * [34] Jie Li, Yi-Pu Wang, Wei-Jiang Wu, Shi-Yao Zhu, and JQ You,  “Quantum network with magnonic and mechanical nodes”,  PRX Quantum, vol. 2, no. 4, pp. 040344, 2021.
  * [35] Banoj Kumar Nayak and Eyal Buks,  “Polarization-selective magneto-optical modulation”,  Journal of Applied Physics, vol. 132, no. 19, pp. 193905, 2022.
  * [36] Zeng-Xing Liu and Hao Xiong,  “Magnon laser based on brillouin light scattering”,  Optics Letters, vol. 45, no. 19, pp. 5452–5455, 2020.
  * [37] JR Sandercock and W Wettling,  “Light scattering from thermal acoustic magnons in yttrium iron garnet”,  Solid State Communications, vol. 13, no. 10, pp. 1729–1732, 1973\.
  * [38] HL Hu and FR Morgenthaler,  “Strong infrared-light scattering from coherent spin waves in yttrium iron garnet”,  Applied Physics Letters, vol. 18, no. 7, pp. 307–310, 1971.
  * [39] E Ghasemian,  “Dissipative dynamics of optomagnonic nonclassical features via anti-stokes optical pulses: squeezing, blockade, anti-correlation, and entanglement”,  Scientific Reports, vol. 13, no. 1, pp. 12757, 2023.
  * [40] Silvia Viola Kusminskiy, Hong X Tang, and Florian Marquardt,  “Coupled spin-light dynamics in cavity optomagnonics”,  Physical Review A, vol. 94, no. 3, pp. 033821, 2016.
  * [41] W Wettling, MG Cottam, and JR Sandercock,  “The relation between one-magnon light scattering and the complex magneto-optic effects in yig”,  Journal of Physics C: Solid State Physics, vol. 8, no. 2, pp. 211, 1975.
  * [42] Michael G Cottam and David J Lockwood,  Light scattering in magnetic solids,  Wiley New York, 1986.
  * [43] Tianyu Liu, Xufeng Zhang, Hong X Tang, and Michael E Flatté,  “Optomagnonics in magnetic solids”,  Physical Review B, vol. 94, no. 6, pp. 060405, 2016.
  * [44] JA Haigh, A Nunnenkamp, and AJ Ramsay,  “Polarization dependent scattering in cavity optomagnonics”,  Physical Review Letters, vol. 127, no. 14, pp. 143601, 2021.
  * [45] R Hisatomi, A Noguchi, R Yamazaki, Y Nakata, A Gloppe, Y Nakamura, and K Usami,  “Helicity-changing brillouin light scattering by magnons in a ferromagnetic crystal”,  Physical Review Letters, vol. 123, no. 20, pp. 207401, 2019.
  * [46] Ryusuke Hisatomi, Alto Osada, Yutaka Tabuchi, Toyofumi Ishikawa, Atsushi Noguchi, Rekishu Yamazaki, Koji Usami, and Yasunobu Nakamura,  “Bidirectional conversion between microwave and light via ferromagnetic magnons”,  Physical Review B, vol. 93, no. 17, pp. 174427, 2016.
  * [47] Wei-Jiang Wu, Yi-Pu Wang, Jin-Ze Wu, Jie Li, and JQ You,  “Remote magnon entanglement between two massive ferrimagnetic spheres via cavity optomagnonics”,  Physical Review A, vol. 104, no. 2, pp. 023711, 2021.
  * [48] Victor Lazzarini, Joseph Timoney, and Thomas Lysaght,  “Asymmetric-spectra methods for adaptive fm synthesis”,  MURAL - Maynooth University Research Archive Library, 2008.
  * [49] Wei Li, Wen Ting Wang, Li Xian Wang, and Ning Hua Zhu,  “Optical vector network analyzer based on single-sideband modulation and segmental measurement”,  IEEE Photonics Journal, vol. 6, no. 2, pp. 1–8, 2014.
  * [50] S Shimotsu, S Oikawa, T Saitou, N Mitsugi, K Kubodera, T Kawanishi, and M Izutsu,  “Single side-band modulation performance of a linbo 3 integrated modulator consisting of four-phase modulator waveguides”,  IEEE Photonics Technology Letters, vol. 13, no. 4, pp. 364–366, 2001.
  * [51] Eyal Buks,  “Intermode coupling in a fiber loop laser at low temperatures”,  Journal of Lightwave Technology, vol. 42, pp. 2951, 2024.
  * [52] Mehmet C Onbasli, Lukáš Beran, Martin Zahradník, Miroslav Kučera, Roman Antoš, Jan Mistrík, Gerald F Dionne, Martin Veis, and Caroline A Ross,  “Optical and magneto-optical behavior of cerium yttrium iron garnet thin films at wavelengths of 200–1770 nm”,  Scientific reports, vol. 6, 2016.
  * [53] PC Fletcher and RO Bell,  “Ferrimagnetic resonance modes in spheres”,  Journal of Applied Physics, vol. 30, no. 5, pp. 687–698, 1959.
  * [54] S Sharma,  Cavity optomagnonics: Manipulating magnetism by light,  PhD thesis, Delft University of Technology, 2019.
  * [55] Tien Liang Jin,  “Design of a yig-tuned oscillator”,  New Jersey Institute of Technology, 1974.
  * [56] Ml Freiser,  “A survey of magnetooptic effects”,  IEEE Transactions on magnetics, vol. 4, no. 2, pp. 152–161, 1968\.
  * [57] Allan D Boardman and Ming Xie,  “Magneto-optics: a critical review”,  Introduction to Complex Mediums for Optics and Electromagnetics, vol. 123, pp. 197, 2003.
  * [58] Allan D Boardman and Larry Velasco,  “Gyroelectric cubic-quintic dissipative solitons”,  IEEE Journal of selected topics in quantum electronics, vol. 12, no. 3, pp. 388–397, 2006.
  * [59] DL Wood and JP Remeika,  “Effect of impurities on the optical properties of yttrium iron garnet”,  Journal of Applied Physics, vol. 38, no. 3, pp. 1038–1045, 1967\.
  * [60] Aleksandr Mikhailovich Prokhorov and Smolenskii,  “Optical phenomena in thin-film magnetic waveguides and their technical application”,  Soviet Physics Uspekhi, vol. 27, no. 5, pp. 339, 1984.
  * [61] Anil Prabhakar and Daniel D Stancil,  Spin waves: Theory and applications, vol. 5,  Springer, 2009.
  * [62] Eyal Buks,  Quantum mechanics - Lecture Notes,  http://buks.net.technion.ac.il/teaching/, 2024.
  * [63] Eyal Buks,  “Tunable multimode lasing in a fiber ring”,  Physical Review Applied, vol. 19, no. 5, pp. L051001, 2023.
  * [64] S Yamashita and K Hotate,  “Multiwavelength erbium-doped fibre laser using intracavity etalon and cooled by liquid nitrogen”,  Electronics Letters, vol. 32, no. 14, pp. 1298–1299, 1996.
  * [65] Haochong Liu, Wei He, Yantao Liu, Yunhui Dong, and Lianqing Zhu,  “Erbium-doped fiber laser based on femtosecond laser inscribed fbg through fiber coating for strain sensing in liquid nitrogen environment”,  Optical Fiber Technology, vol. 72, pp. 102988, 2022.
  * [66] J Lopez, H Kerbertt, M Plata, E Hernandez, and S Stepanov,  “Two-wave mixing in erbium-doped-fibers with spectral-hole burning at 77k”,  Journal of Optics, vol. 22, no. 8, pp. 085401, 2020.
  * [67] Julien Le Gouët, Jérémy Oudin, Philippe Perrault, Alaeddine Abbes, Alice Odier, and Alizée Dubois,  “On the effect of low temperatures on the maximum output power of a coherent erbium-doped fiber amplifier”,  Journal of Lightwave Technology, vol. 37, no. 14, pp. 3611–3619, 2019.
  * [68] Shi-Hai Wei, Bo Jing, Xue-Ying Zhang, Jin-Yu Liao, Hao Li, Li-Xing You, Zhen Wang, You Wang, Guang-Wei Deng, Hai-Zhi Song, et al.,  “Storage of 1650 modes of single photons at telecom wavelength”,  arXiv:2209.00802, 2022.
  * [69] Antonio Ortu, Jelena V Rakonjac, Adrian Holzäpfel, Alessandro Seri, Samuele Grandi, Margherita Mazzera, Hugues de Riedmatten, and Mikael Afzelius,  “Multimode capacity of atomic-frequency comb quantum memories”,  Quantum Science and Technology, vol. 7, no. 3, pp. 035024, 2022\.
  * [70] Duan-Cheng Liu, Pei-Yun Li, Tian-Xiang Zhu, Liang Zheng, Jian-Yin Huang, Zong-Quan Zhou, Chuan-Feng Li, and Guang-Can Guo,  “On-demand storage of photonic qubits at telecom wavelengths”,  arXiv:2201.03692, 2022.
  * [71] Lucile Veissier, Mohsen Falamarzi, Thomas Lutz, Erhan Saglamyurek, Charles W Thiel, Rufus L Cone, and Wolfgang Tittel,  “Optical decoherence and spectral diffusion in an erbium-doped silica glass fiber featuring long-lived spin sublevels”,  Physical Review B, vol. 94, no. 19, pp. 195138, 2016.
  * [72] Sara Shafiei, Erhan Saglamyurek, and Daniel Oblak,  “Hour-long decay-time of erbium spins in an optical fiber at milli-kelvin temperatures”,  in Quantum Information and Measurement. Optical Society of America, 2021, pp. F2A–4.
  * [73] A Liem, J Limpert, T Schreiber, M Reich, H Zellmer, A Tunnermann, A Carter, and K Tankala,  “High power linearly polarized fiber laser”,  in Conference on Lasers and Electro-Optics. Optica Publishing Group, 2004, p. CMS4.
  * [74] Javier Trujillo-Bueno, Fernando Moreno-Insertis, and Francisco Sánchez,  Astrophysical spectropolarimetry,  Cambridge University Press, 2002.
  * [75] Denis A Belyaev, Konstantin B Yushkov, Sergey P Anikin, Yuri S Dobrolenskiy, Aleksander Laskin, Sergey N Mantsevich, Vladimir Ya Molchanov, Sergey A Potanin, and Oleg I Korablev,  “Compact acousto-optic imaging spectro-polarimeter for mineralogical investigations in the near infrared”,  Optics Express, vol. 25, no. 21, pp. 25980–25991, 2017.


# Summed Parallel Infinite Impulse Response (SPIIR) Filters For Low-Latency
Gravitational Wave Detection

Shaun Hooper<EMAIL_ADDRESS>Australian International Gravitational
Research Centre, School of Physics, University of Western Australia, 35
Stirling Hwy, Crawley, WA 6009, Australia  ICRAR-Fairway M468, School of
Physics, The University of Western Australia, Crawley, WA 6009, Australia
Shin Kee Chung  Australian International Gravitational Research Centre, School
of Physics, University of Western Australia, 35 Stirling Hwy, Crawley, WA
6009, Australia  ICRAR-Fairway M468, School of Physics, The University of
Western Australia, Crawley, WA 6009, Australia  Jing Luan  Theoretical
Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125,
USA  David Blair  Australian International Gravitational Research Centre,
School of Physics, University of Western Australia, 35 Stirling Hwy, Crawley,
WA 6009, Australia  Yanbei Chen  Theoretical Astrophysics 350-17, California
Institute of Technology, Pasadena, CA 91125, USA  Linqing Wen
<EMAIL_ADDRESS>Australian International Gravitational Research Centre,
School of Physics, University of Western Australia, 35 Stirling Hwy, Crawley,
WA 6009, Australia  ICRAR-Fairway M468, School of Physics, The University of
Western Australia, Crawley, WA 6009, Australia

###### Abstract

With the upgrade of current gravitational wave detectors, the first detection
of gravitational wave signals is expected to occur in the next decade. Low-
latency gravitational wave triggers will be necessary to make fast follow-up
electromagnetic observations of events related to their source, e.g., prompt
optical emission associated with short gamma-ray bursts. In this paper we
present a new time-domain low-latency algorithm for identifying the presence
of gravitational waves produced by compact binary coalescence events in noisy
detector data. Our method calculates the signal to noise ratio from the
summation of a bank of parallel infinite impulse response (IIR) filters. We
show that our _summed parallel infinite impulse response_ (SPIIR) method can
retrieve the signal to noise ratio to greater than 99% of that produced from
the optimal matched filter. We emphasise the benefits of the SPIIR method for
advanced detectors, which will require larger template banks.

General Relativity, Gravitational Waves

###### pacs:

04.25.Nx, 04.30.Db, 04.80.Cc, 04.80.Nn, 95.55.Ym, 95.85.Sz

##  I Introduction

The interferometric gravitational wave (GW) detectors LIGO LIG , and Virgo VIR
have reached a sensitivity at which the detection of GWs is possible. The LIGO
detectors are currently undergoing a major upgrade to Advanced LIGO, for which
the sensitivity will be improved ten fold relative to Initial LIGO Smith and
for the LIGO Scientific Collaboration (2009). Hence Advanced LIGO will be able
to detect GW (GW) sources within a volume of space one thousand times larger
than that of initial LIGO, out to $\sim$200-300$\,\mathrm{Mpc}$ Adv (2007).

The emission of GWs produced by compact binary coalescence (CBC) can be
modelled with a high degree accuracy Abbott et al. (2009). When two compact
bodies, such as neutron stars or black holes are in orbit, Einstein’s
equations predict the generation of GWs. As the bodies spiral towards each
other a GW is created that increases in frequency over time until the bodies
merge, following what is known as the inspiral waveform. Ground based
detectors have frequency passbands that allow them to be sensitive to the
final stages of such events up to a total system masses of several hundred
$M_{\odot}$.

Neutron star binary mergers are widely thought to be the progenitors of short
hard gamma-ray bursts (short GRBs) Fox et al. (2005); Nakar (2007). The delay
between the final GW emission and the onset of the GRB is estimated to be as
short as 0.1 seconds or as long as tens to hundreds of seconds van Putten
(2009); Zhang and Mészáros (2004). The electromagnetic emission of the GRB
event is not well understood. Related to the initial GRB there is thought to
be a prompt emission in the X-ray and optical wavelengths followed by a
delayed afterglow of cascading wavelengths. Prompt optical emission may occur
tens to hundreds of seconds after the initial burst. The low-latency detection
of the GW associated with a neutron star merger could lead to the localisation
of a GRB source event on the sky, enabling fast moving telescopes to observe
the prompt optical emission. Data collected from a multitude of sources — GWs,
gamma-rays, X-rays and optical counterparts of the GRB — will lead to maximum
insight into these highly energetic events.

The standard strategy for searching for the existence of inspiral waveforms in
the detector data is based on matched filtering Abbott et al. (2009) (and
references therein). This method, based on Wiener optimal filtering, is a
correlation of an expected inspiral waveform template and the detector data,
weighted by the inverse noise spectral density of the detector Wainstein and
Zubakov (1962). In order to save computational costs, this correlation is
performed in the frequency domain, via a Fourier transform of a finite segment
of detector data. In previous LIGO searches, the detector data is split up
into “science blocks”, which are further divided into “data segments” chosen
to be at least twice the length of the longest waveform in the template bank
Abbott et al. (2008). Each proceeding data segment is chosen to overlap the
previous one by 50%. Each segment therefore must be matched filtered in a time
that is half the length of the segment for a real-time analysis, that is, the
filter output rate is equal to the data input rate. In this case, the matched
filter process has a minimum latency (from signal arrival to signal detection)
that is proportional to the longest template (see Luan (2011) for more
details). Advanced LIGO will have an increased bandwidth over Initial LIGO,
with the lower bound dropping from 40$\,\mathrm{Hz}$ to 10$\,\mathrm{Hz}$ Adv
(2007). GW signals from CBC events spend much more time at these lower
frequencies. Hence waveforms used for matched filtering in Advanced LIGO will
be much longer (1000s of seconds). This in turn means the segment length will
be increased, further increasing the latency. The latency of this method to
produce GW triggers is longer than the time to onset of prompt optical
emission after coalescence (10s to 100s of seconds). After this amount of
time, the early electromagnetic counterpart of a GRB event will be
significantly faded, and may be missed by telescopes altogether.

A low-latency GW detection method is required to trigger follow-up
electromagnetic observations of the prompt optical emission. So far two
frequency domain methods have been developed to solve this issue. The VIRGO
group has produced a low-latency pipeline based on _Multi-Band Template
Analysis_ (MBTA) Buskulic et al. (2010), and LIGO is also working on a new
method, _Low-Latency On-line Inspiral Data_ analysis (LLOID) method. In MBTA
the matched filtering technique is split over two frequency bands, and the
output is coherently added, reducing latency. A latency of less than 3 minutes
until the availability of a trigger using this method has been achieved
Buskulic et al. (2010). Low-latency in the LLOID method is achieved by first
down-sampling the incoming data into multiple streams and then applying
frequency domain finite impulse response (FIR) filters gst . The computational
cost of this pipeline is reduced by decreasing the number of templates via
singular value decomposition Cannon et al. (2010).

We introduce a new method to detect CBC signals in the time domain using
infinite impulse response (IIR) filters. Approximating an inspiral waveform by
a summation of time shifted exponentially increasing sinusoids enables us to
construct a bank of parallel single-pole IIR filters. Each IIR filter acts as
a narrow bandpass filter. When each appropriately delayed IIR filter is added
the coherent output approximates the matched filter output of the exact
waveforms. We call this the _summed parallel infinite impulse response_
(SPIIR) method. Figure 1 visually demonstrates the idea of using a bank of IIR
filters as narrow bandpass filters.

Figure 1: A schematic overview of the SPIIR method. The input is split into
different channels, time delayed by an amount $d$, then passed through a
narrow bandpass IIR filters, each with a different central frequency $f$.
Finally the output of each individual IIR filter is summed, giving the output
of the SPIIR method.

For a full explanation of the mathematical principles, see Luan (2011). In
this follow up paper, we numerically address the issues essential to the
practical use of this method for the upcoming advanced detectors. We calculate
the filter coefficients and demonstrate via numerical simulations how well our
method approximates the optimal matched filter as a function the number of
filters per bank using a range of parameters. We also show that the detection
rate of the SPIIR method is very similar to that of the matched filter method.
It has been shown theoretically that in order to get the same latency as the
SPIIR method, the frequency domain matched filter method would require greater
computational resources Luan (2011).

The structure of this paper is as follows: In section II we will go through
the formal introduction of the inspiral waveform and matched filtering, and
how to get from the continuous frequency domain matched filter to the time
domain discrete matched filter. This will lead to a demonstration on how it is
possible to approximate an inspiral signal by a sum of exponentially
increasing sinusoids. The methodology is explained in Section III and will
cover how we set up our simulation to test the efficiency of the SPIIR method
as opposed to the frequency domain matched filter. Section IV will analyse the
results of the simulation and Section V will discuss the implications of these
results for advanced detectors.

##  II Methodology

Gravitational wave interferometers output the strain induced by gravitational
waves incident on the detector, as well as inherent noise. In unitless strain,
the detector output will be,

$\displaystyle s(t)=\begin{cases}n(t)&\text{if signal is absent}\\\
n(t)+h(t)&\text{if signal is present}\end{cases}$ (1)

where $n(t)$ is the noise inherent in the detector. The sensitivity of the
instrument can be characterized by the (one-sided) strain power spectral
density $S_{n}(f)$,

$\displaystyle\left\langle\tilde{n}(f)\tilde{n}^{*}(f^{\prime})\right\rangle=\frac{1}{2}S_{n}(f)\delta(f-f^{\prime})$
(2)

where the tilde represents the forward Fourier transform,

$\displaystyle\tilde{q}(f)=\int_{-\infty}^{\infty}q(t)e^{-2\pi ift}dt.$ (3)

###  II.1 The Inspiral Waveform

The gravitational-wave strain incident at the interferometer is given by

$\displaystyle
h(t)=F_{+}(\theta,\phi,\psi)h_{+}(t)+F_{\times}(\theta,\phi,\psi)h_{\times}(t)$
(4)

where the detector response functions $F_{+}$ and $F_{\times}$ are functions
of $(\theta,\phi)$ \- the standard spherical polar coordinates measured with
respect to the Earth’s fixed frame, and $\psi$ is the polarisation angle. The
detector response function can be found in Anderson et al. (2001). The $+$ and
$\times$ polarisations of the waveform are,

$\displaystyle h_{+}(t)$
$\displaystyle=\left(\frac{1+\cos^{2}\iota}{2}\right)A(t)\cos\phi(t)$ (5)
$\displaystyle h_{\times}(t)$
$\displaystyle=\left(\cos\iota\right)A(t)\sin\phi(t)$ (6)

For non-spinning binaries with a chirp mass
$\mathcal{M}=((m_{1}m_{2})^{3}/(m_{1}+m_{2}))^{1/5}$ in the range of
$1-3M_{\odot}$ — we will hereafter assume — the waveforms can be modelled to
very high accuracy using the Restricted post-Newtonian (PN) expansion Abbott
et al. (2004); Blanchet et al. (1995, 1996) in the LIGO band (assumed to be
10-1500 $\,\mathrm{Hz}$ for advanced LIGO). For restricted waveforms, only the
leading order of the amplitude $A(t)$ is taken,

$\displaystyle A(t)$
$\displaystyle=\frac{G\mathcal{M}}{Dc^{2}}\left(\frac{t_{c}-t}{5G\mathcal{M}/c^{3}}\right)^{-1/4}$
(7)

and the post-Newtonian phase $\phi(t)$ is given by

$\displaystyle\phi(t)$
$\displaystyle=\phi_{c}-2\left(\frac{t_{c}-t}{G\mathcal{M}/c^{3}}\right)^{5/8}+\mbox{higher
order terms}$ (8)

In addition to the source masses $m_{1},m_{2}$, there are several unknown
parameters; the time of coalescence $t_{c}$, the phase at coalescence
$\phi_{0}$, distance from observer to source $D$, the inclination angle of the
binary’s orbital plane relative the line of sight $\iota$, and the
polarisation angle $\psi$. However by using the linear combination
trigonometric identity, one can re-express the strain (4) by splitting the
scaling factor due to distance, sky location and orientation to the mass
dependant time evolution of the waveform Fairhurst and Brady (2008),

$\displaystyle h(t)$ $\displaystyle=\frac{1\,\mathrm{Mpc}}{D_{\rm
eff}}\left[h_{c}(t)\cos\phi_{0}+h_{s}(t)\sin\phi_{0}\right]$ (9)

where the scalar factor $D_{\rm eff}$ is,

$\displaystyle D_{\rm
eff}=\frac{D}{\sqrt{F_{+}^{2}\left(1+\cos^{2}\iota\right)^{2}/4+F_{\times}^{2}\left(\cos\iota\right)^{2}}}$
(10)

which gives $\phi_{0}$, an unknown phase as,

$\displaystyle\phi_{0}=\phi_{c}+\arctan\frac{F_{\times}\left(2\cos\iota\right)}{F_{+}\left(1+\cos^{2}\iota\right)}$
(11)

We now define the terms $h_{c}$ and $h_{s}$ as the waveform at $\phi_{0}=0$
and $\phi_{\pi/2}$, scaled at 1$\,\mathrm{Mpc}$ as the so called “cosine” and
“sine” phases Brady and Fairhurst (2008),

$\displaystyle h_{c}(t)$ $\displaystyle=A_{1\,\mathrm{Mpc}}(t)\cos\phi(t)$
(12) $\displaystyle h_{s}(t)$
$\displaystyle=A_{1\,\mathrm{Mpc}}(t)\sin\phi(t)$ (13)

###  II.2 The Matched Filter

The matched filter $Q$ is a linear operator that maximises the ratio of
“signal” to “noise” present in the detector data $s$ Allen et al. (2005). It
is denoted by,

$\displaystyle z(t)$
$\displaystyle=2\int_{-\infty}^{\infty}\frac{\tilde{s}(f)\tilde{Q}^{*}(f)}{S_{n}(|f|)}e^{2\pi
ift}df=\left(s(t)\left|Q\right)\right.$ (14)

Where we have also defined the inner product $\left(a\left|b\right)\right.$.
The signal to noise ratio (SNR) is _generally_ defined as the ratio of
observed filter output to it’s expected root-mean square flucations or
standard deviation,

$\displaystyle\mathrm{SNR}$
$\displaystyle=\frac{z}{\sqrt{\left\langle(z-\left\langle
z\right\rangle)^{2}\right\rangle}}=\frac{z}{\sqrt{\left\langle
z\right\rangle^{2}}}=\frac{z}{\sqrt{\left(Q\left|Q\right)\right.}}$ (15)

Note that in the absence of a signal,
$\left\langle\mathrm{SNR}\right\rangle=0$ and
$\left\langle(\mathrm{SNR})^{2}\right\rangle=1$ independent of the
normalisation of the filter $Q$.

###  II.3 Two-Phase Filter

A convenient way to search for the unknown phase constant $\phi_{0}$ is to
filter both phases $h_{c}$ and $h_{s}$ separately and then combined to form a
complex signal. The two-phase filter is defined as,

$\displaystyle z(t)$
$\displaystyle=\left(s(t)\left|h_{c}\right)\right.+i\left(s(t)\left|h_{s}\right)\right.$
(16a)
$\displaystyle\begin{split}&=2\int_{-\infty}^{\infty}\frac{\tilde{s}(f)\tilde{h}_{c}^{*}(f)}{S_{n}(f)}e^{2\pi
ift}df\\\
&\qquad\qquad+i2\int_{-\infty}^{\infty}\frac{\tilde{s}(f)\tilde{h}_{s}^{*}(f)}{S_{n}(|f|)}e^{2\pi
ift}df\end{split}$ (16b)

The advantage of using the phases $h_{c,s}$ is that in the stationary phase
approximation Droz et al. (1999), $h_{c}$ and $h_{s}$ are exactly orthogonal
($\left(h_{c}\left|h_{s}\right)\right.=\left(h_{s}\left|h_{s}\right)\right.$,
$\left(h_{c}\left|h_{s}\right)\right.=0$). It then follows,
$\tilde{h}_{c}(f)=i\tilde{h}_{s}(f)$ for $f>0$. Generally, this is applied to
(16) to give the two-phase matched filter as,

$\displaystyle z(t)$
$\displaystyle=4\int_{0}^{\infty}\frac{\tilde{s}(f)\tilde{h}_{c}^{*}(f)}{S_{n}(|f|)}e^{2\pi
ift}df$ (17)

However in this paper, we prefer to maintain the form of the two-phase filter
in (16). In convention with the field, the amplitude signal to noise ratio of
the (quadrature) matched filter is defined as the absolute value of the two-
phase filter, divided by a normalisation constant that is equal to standard
deviation of the real and imaginary parts of the two-phase filter,

$\displaystyle\rho(t)=\frac{|z(t)|}{\sigma}$ (18)

where $\sigma^{2}$ is,

$\displaystyle\sigma^{2}=2\int_{-\infty}^{\infty}\frac{\left|\tilde{h}_{c}(f)\right|^{2}}{S_{n}(f)}df=\left(h_{c}\left|h_{c}\right)\right.$
(19)

Note that in the in the absence of a signal (just noise), the SNR $\rho$ (18)
is Rayleigh distributed with mean $\sqrt{\pi/2}$ and variance $1$, which is
identical to the Chi-distribution with two degrees of freedom (one for each of
the phases). This of course implies that the SNR squared, $\rho^{2}$ is Chi-
square distributed with two degrees of freedom. Hence the probability of
finding an SNR value greater than $\rho_{*}$ is Brady and Fairhurst (2008),

$\displaystyle P(\rho^{2}>\rho_{*}^{2})=e^{-\rho_{*}^{2}/2}.$ (20)

###  II.4 Digital Time Domain Filtering

The two-phase matched filter 16 is a cross correlation of phase $h_{c,s}(t)$
and the detector output $s(t)$, weighted by the inverse noise spectral density
$S_{n}(f)$. By defining the quantity $x$ as the _over_ -whitened strain data,

$\displaystyle x(t)=\int_{-\infty}^{\infty}\frac{\tilde{s}(f)}{S(f)}e^{2\pi
ift}df$ (21)

we can use the cross-correlation theorem to define the two-phase matched
filter in the time domain,

$\displaystyle z(t)$
$\displaystyle=2\int_{-\infty}^{t}x(t^{\prime})h_{c}(t^{\prime}-t)dt^{\prime}+i2\int_{-\infty}^{t}x(t^{\prime})h_{s}(t^{\prime}-t)dt^{\prime}$
(22)
$\displaystyle=2\int_{-\infty}^{t}x(t^{\prime})\hat{h}(t^{\prime}-t)dt^{\prime}$
(23)

where $\hat{h}=h_{c}(t)+ih_{s}(t)=A(t)e^{i\phi(t)}$.

The discrete form of the continuous time domain matched filter (23) is,

$z_{k}=2\sum_{j=-\infty}^{k}x_{j}\hat{h}_{j-k}\Delta t$ (24)

where $t=k\Delta t$. In practise, the inspiral waveform template $h_{i}$ is
bounded (because the detector is only sensitive over a bandwidth), and the
summation becomes finite, making this a _finite impulse response_ (FIR)
filter.

###  II.5 Infinite Impulse Response Filter

Now let us introduce an alternative digital filter, the _infinite impulse
response_ (IIR) filter. The difference equation of a general IIR filter is,

$\displaystyle y_{k}=\sum_{n=1}^{N}a_{n}y_{k-n}+\sum_{m=0}^{M}b_{m}x_{k-m}$
(25)

where $y_{k}$ is the filter output at time step $k$, ($t=k\Delta t$), $x_{k}$
is the filter input, and $a$’s and $b$’s are complex coefficients.

Examples of IIR filters in common usage are Chebyshev, Butterworth and
elliptic filters. IIR filters use much less computational resources than an
equivalent FIR filter. This is because they have “memory” — the previous
outputs are fed back into the filter. However _digital_ IIR filter design is a
more complex process than FIR design. Obtaining the coefficients is usually
done by first constructing an equivalent analog filter and applying well-known
methods, such as the bi-linear transform or impulse invariance. Multiple IIR
filters used together have different forms, such as direct form I & II,
cascade (series) and parallel. In a series configuration, the overall transfer
function is the multiplication of each IIR filter transfer function. In a
parallel bank of IIR filters, where the output is summed together, the overall
transfer function is the summation of the different transfer functions.

First, let’s analyse the simplest single-pole IIR filter. The difference
equation of this filter is

$y_{k}=a_{1}y_{k-1}+b_{0}x_{k}.$ (26)

Figure 2: A signal processing schematic showing the flow of data through a
digital single-pole IIR filter. The input, $x_{k}$ is multiplied by a complex
constant $b_{0}$, then added to the previous output that has been multiplied
by another complex constant $a_{1}$, resulting in the current output $y_{k}$.
It should be noted that this filter, in principle, should be have been run
forever.

A solution to this first-order linear inhomogeneous difference equation is

$\displaystyle y_{k}=\sum_{j=-\infty}^{k}x_{j}b_{0}a_{1}^{k-j}.$ (27)

By defining the complex coefficient $a_{1}$ in the form,

$\displaystyle a_{1}=e^{-(\gamma+i\omega)\Delta t}$ (28)

and comparing (24) and (27), it is easy to see that the output of the simple
filter (26) is the cross-correlation of $x_{k}$ and complex sinusoid $u_{n}$
with frequency $\omega$ and a magnitude that increases with an exponential
factor $\gamma$ for $n<0$:

$\displaystyle u_{n}=b_{0}e^{(\gamma+i\omega)n\Delta t}\Theta(-n)$ (29)

where $\Theta(-n)$ is the Heaviside function.

###  II.6 Approximation to an inspiral waveform

Since $\phi(t)$ is not linear in time, a complex sinusoid (29) cannot
approximate the $h_{c,s}$ phases of the inspiral waveform
$\hat{h}(t)=A(t)e^{i\phi(t)}$. However we can easily linearise the phases by a
first-order Taylor expansion about the time $t_{l}^{*}$:

$\displaystyle A(t)e^{i\phi(t)}\simeq
A(t_{l}^{*})e^{i\phi(t_{l}^{*})+i\dot{\phi}(t_{l}^{*})(t-t_{l}^{*})};$ (30)

since the amplitude $A(t)$ does not increase at the same rate as $\phi(t)$,
only a linear expansion of $\phi(t)$ is required. Multiplying by the window
function $e^{\gamma_{l}(t-t_{l})}\Theta(t_{l}-t)$ makes this approximation an
exponentially increasing constant frequency complex sinusoid with cutoff time
$t_{l}$:

$\displaystyle u_{l}(t)$
$\displaystyle=A(t_{l}^{*})e^{i(\phi(t_{l}^{*})+\dot{\phi}(t_{l}^{*})(t_{l}-t_{l}^{*}))}e^{(\gamma_{l}+i\dot{\phi}(t_{l}^{*})(t-t_{l})}\Theta(t_{l}-t).$
(31)

The expansion point $t_{l}^{*}$ is chosen to be near the cutoff time,
$t_{l}^{*}=t_{l}-\alpha T_{l}$, where $\alpha$ is a tunable parameter and the
interval $T_{l}$ is the duration in which the approximation is valid:

$\displaystyle|\frac{1}{2}\ddot{\phi}(t_{l})T_{l}^{2}|=\epsilon<1$ (32)

and $\epsilon$ is a tunable parameter chosen to be to small. Equation (31)
implies that the coefficient $b_{0}$ for the $l$th complex sinusoid is,

$\displaystyle b_{0,l}$
$\displaystyle=A(t_{l}^{*})e^{i(\phi(t_{l}^{*})+\dot{\phi}(t_{l}^{*})(t_{l}-t_{l}^{*}))}$
(33)

and the frequency $\omega_{l}=\dot{\phi}(t_{l}^{*})$.

In this paper, we chose the cutoff time $t_{l}$ of the first sinusoid to
correspond to the time at which the waveform has the highest frequency
detectable by the LIGO detector band. The next sinusoid is chosen by moving to
an earlier time, $t_{l+1}=t_{l}-T_{l}$. Since we want the $l$th sinusoid to be
mostly present on the interval $t_{l}-T_{l}<t<t_{l}$, we choose the damping
factor to be $\gamma_{l}=\beta/T_{l}$, where $\beta$ is a tunable parameter.
This procedure is repeated until the time $t_{l}$ corresponds to a time in the
waveform that has frequency below the LIGO detector band. Hence the number of
sinusoids is dependent on the value of $\epsilon$, the rate of frequency
change ($\ddot{\phi}(t)$), which is dependent on the masses of the system, and
the detector bandwidth. For more information on this procedure, see Luan
(2011).

We can now approximate the phases $\hat{h}(t)=A(t)e^{i\phi(t)}$ by an addition
of a series of damped sinusoids $u(t)$ with cutoff times $t_{l}$:

$\displaystyle A(t)e^{i\phi(t)}\simeq U(t)$ $\displaystyle=\sum_{l}u_{l}(t)$
$\displaystyle=\sum_{l}b_{0,l}e^{(\gamma_{l}+i\omega_{l})(t-t_{l})}\Theta(t_{l}-t).$
(34)

Figure 3 shows an illustration of how damped constant frequency sinusoids can
add to give an inspiral like waveform.

Figure 3: An illustrative diagram demonstrating the ability to linearly sum
exponentially increasing constant frequency sinusoids to approximate an
inspiral like waveform. The top three panels (a-c) show three example
sinusoids with different damping, frequency and cutoff time factors. Panel (d)
shows the linear addition of all the sinusoids (at different scales). Panel
(e) shows the exact inspiral-like waveform. Note that this figure is only for
illustrative purposes.

###  II.7 Summed Parallel IIR filtering

Each complex sinusoid $u_{l}(t)$ in equation (II.6) can be searched for in the
data $x$ using the single pole IIR filter (26). Here the cutoff time is
incorporated by running each filter on a delay, $d_{l}=t_{l}/\Delta t$. The
output of the $l$th filter at time $k$ is

$\displaystyle y_{k,l}=a_{1,l}y_{k-1,l}+b_{0,l}x_{k-d_{l}}.$ (35)

The linear summation of the output of _all_ filters is the cross-correlation
of the data $x$ and the approximate waveform $U(t)$ in (II.6):

$\displaystyle z_{k}\simeq 2\Delta t\sum_{l}y_{k,l}.$ (36)

Here $z$ is equivalent to the value computed by the discrete time domain two
phase filter (24) when using a template $\hat{h}(t)=U(t)$. From equation (18),
it follows that the absolute value of the summation (36) divided by
$\sigma_{U}$ is the SNR, which we term the output of the _Summed Parallel
Infinite Impulse Response_ (SPIIR). The normalisation factor $\sigma_{U}$ is
defined as

$\displaystyle\sigma_{U}^{2}=4\int_{0}^{\infty}\frac{\left|\tilde{U}_{\Re}(f)\right|^{2}}{S_{n}(f)}df.$
(37)

Where $\tilde{U}_{\Re}(f)$ is the Fourier transform of the real part of
$U(t)$, (which approximates $h_{c}(t)$). The similarity of the SPIIR output
and the matched filter output will depend on how well $U(t)$ approximates the
given template.

##  III Implementation for Performance Testing

###  III.1 IIR bank construction

To confirm the ability of the SPIIR method to recover a good SNR, it is first
required to show that the approximate inspiral waveform (II.6) is a good
“match” to the theoretical inspiral waveform (9). We define the _overlap_
$\Delta$ as the inner product of the approximate waveform $U$ and the template
$h$:

$\displaystyle\Delta=\left(\frac{h}{\sqrt{\left(h\left|h\right)\right.}}\left|\frac{U}{\sqrt{\left(U\left|U\right)\right.}}\right)\right.=\frac{\left(h\left|U\right)\right.}{\sqrt{\left(h\left|h\right)\right.\left(U\left|U\right)\right.}}$
(38)

We initially approximate a canonical 2PN 1.4-1.4 $M_{\odot}$ inspiral waveform
band limited to 10-1500$\,\mathrm{Hz}$ using the value of the tunable
parameters $\epsilon$, $\alpha$ and $\beta$ to be consistent with the high
overlap results of Luan (2011). With some minor variation of their values, we
aim to recover the highest overlap possible. Once a good choice of $\alpha$
and $\beta$ is found for the 2PN 1.4-1.4 $M_{\odot}$ template, we use the same
values for other templates, but vary the value $\epsilon$ (and consequently
the number of IIR filters in each bank) to see the effect on overlap.

###  III.2 Detector Data Simulation

To test the detection efficiency of the SPIIR method compared to the frequency
domain matched filter, we will filter two mock signals, one for which the
input data is just LIGO-like noise, and the other with the same noise plus an
inspiral waveform injection scaled to represent a source at a chosen effective
distance $D_{\rm eff}$.

For this test, we need to construct a finite segment of detector data to
filter. Because of the IIR filters should in principle should be run for an
infinite length of the input data, we need to run the IIR bank for a finite
“warm-up” period before the output is consistent with that of an IIR filter
that has been running for an infinite amount of time. In practise, we choose
to run each filter for 2 $e$-foldings of time before we accept the output as
being identical to one which has run for an infinite amount of time.
Additionally, since each IIR filter in the bank runs on a delay, the summed
output _of all the IIR filters_ will not be produced until after the longest
delay time ($d_{\rm max}$) has passed. The filter that has the longest delay
($d_{\rm max}$) is also the one that has the longest decay rate $\gamma_{\rm
max}$. In total, the input data must at least $d_{\rm max}+2\gamma_{\rm
max}^{-1}$ in length before any output is produced. Hence the length of the
input data is,

$N_{\rm input}=d_{\rm max}+2\gamma_{\rm max}^{-1}+N_{\rm analysis}$ (39)

where $N_{\rm analysis}$ is the length of analysis period, which we choose to
be 4 seconds. Hence the 4$\,\mathrm{s}$ SPIIR output will tell us whether
there is an injection that ended somewhere within those 4 seconds. At a sample
rate of 4096$\,\mathrm{Hz}$, the analysis period is $N_{\rm analysis}=16834$
data points long. In our simulation, we find $d_{\rm max}=4081683$ and
$2\gamma_{\rm max}^{-1}=149432$, resulting in $N_{\rm input}=4247499$.

####  III.2.1 Noise generation

The LIGO-like noise data is produced by creating a normally distributed white
noise time series of length $N_{\rm input}$, then colouring it by the
theoretical advanced LIGO noise spectrum $S_{n}(f)$ A. We then over-whiten
this time series using equation (21) to produce the waveform-free noise input
data $x$:

$\displaystyle x_{\rm noise}(t)=n^{\rm ow}(t).$ (40)

####  III.2.2 Waveform injection

We create our waveform injections by first producing an inspiral waveform
band-limited between 10 and 1500$\,\mathrm{Hz}$. The injection is padded with
zeros so that it has the length $N_{\rm input}$. The end of the waveform is
chosen so that it finishes somewhere after $d_{m}+2\gamma_{m}^{-1}$ data
points. The injection signal is then over whitened using equation (21). The
over-whitened injection can then be placed in the over-whitened noise signal,

$\displaystyle x_{\rm noise+injection}(t)=x_{\rm noise}(t)+h^{\rm ow}(t).$
(41)

####  III.2.3 Matched filter comparison

As a comparison, we will also perform a frequency domain correlation matched
filter. For this process, since the input data is already over-whitened, it
only needs to be cross-correlated with the waveform. Section II.2 outlines how
this is done. The cosine phase $h_{c}(t)$ gets pre-padded with enough zeros to
get to length $N_{\rm input}$. This ensures that $\tilde{h}_{c}(f)$ has the
same spectral resolution as $\tilde{s}(f)$. The matched filter (17) produces a
time series of $N_{\rm input}$ length. However the first $N_{\rm input}-N_{\rm
analysis}$ data points are erroneous wrap-around caused by the FFT. Only the
interval $[N_{\rm input}-N_{\rm analysis}+1,N_{\rm analysis}]$ is used to
determine if a waveform is present.

###  III.3 Detection Efficiency

To test the detection efficiency of the SPIIR method compared to the
traditional matched filter method we will construct several receiver operating
characteristic (ROC) curves for 2PN 1.4-1.4 $M_{\odot}$ waveforms injected for
different effective distances $D_{\rm eff}$. To create each ROC curve, we
first find the false alarm rate. The false alarm rate is found by realising an
$N_{\rm input}$ length LIGO-like noise time series, filtering this input data,
and analysing the output of the 4$\,\mathrm{s}$ analysis period (the SNR). We
will count this realisation as a false positive if at any point within the 4
seconds the SNR goes over a given SNR threshold. Several thresholds will be
chosen, giving the false positive as a function of threshold. After $>10^{6}$
noise realisations, the false alarm rate is simply the ratio of total number
of false positives to number of noise realisations. Likewise, to see if the
IIR filter doesn’t miss too many true positives, we inject a 2PN 1.4-1.4
$M_{\odot}$ waveform using the prescribed method in III.2.2 for a given
$D_{\rm eff}$ into LIGO-like noise. After filtering, if at any point within
the analysis period the SNR is above a given threshold, this realisation is
counted as a true positive. Again, after $>10^{6}$ noise realisations, we
calculate the detection rate as a ratio of the total number of true positives
to number of realisations. The plot of false alarm rate versus detection rate
gives the ROC curve.

##  IV Results

###  IV.1 Inspiral Waveform Overlap

Starting with the canonical 1.4-1.4 $M_{\odot}$ second order post-Newtonian
binary waveform band limited to be between 10 and 1500 $\,\mathrm{Hz}$ we
found, using the parameters $\epsilon=0.04$, $\alpha=0.99$, $\beta=0.25$ in
the procedure outlined in Section II.6, that can recover an overlap of 99%
using 687 IIR filters.

We find that increasing the value of $\epsilon$ will in general increase the
overlap, as the frequency space is more finely sampled. However there seems to
be a limit, as the damping factor $\gamma$ causes the adjacent IIR filters to
run into each other.

With this choice of $\alpha$ and $\beta$ we are able to recover a high overlap
for different mass pairs as well. Figure 4 shows the overlap as a function of
number of IIR filters for six different mass pairs.

Figure 4: The overlap between the exact inspiral waveform and the approximate
inspiral waveform as a function of number of damped sinusoids. In general the
greater the number of sinusoids per waveform, the greater the overlap. However
the choice of $\gamma_{l}$s greatly affects the overlap.

###  IV.2 Ability to Recover SNR

Figure 5 shows the SNR produced from both the matched filter technique and the
SPIIR method. The input time series is constructed following Section III.2.
The injection of a 2PN 1.4-1.4 $M_{\odot}$ waveform scaled for an effective
distance of 250$\,\mathrm{Mpc}$ is added to LIGO-like noise. The $x$-axis of
the plot is centred about end of the injection ($t=\tau_{c}$), which is
directly in the middle of the analysis period. Around this time, the SNR peaks
to 8.2, which is near the expected value of 7.9 for an injection at this
distance.

Figure 5: The SNR output of both the SPIIR method and a traditional matched
filter method. The plot is centred on $t-\tau_{c}$ where $\tau_{c}$ is the
time at which the injection ends. From the two curves, it is clear that the
SPIIR method can return a very similar SNR to that from the optimal filter.

This plot shows that the SPIIR method is capable of recovering a very similar
SNR to the matched filter at all times.

###  IV.3 Detection Efficiency

We analysed over $10^{6}$ independent noise realisations, for which the
waveform had been injected at $D_{\rm eff}$ of 250, 300, 350, 400
$\,\mathrm{Mpc}$. We performed both IIR filtering and traditional matched
filtering. Figure 6 shows that the SPIIR method recovers most of the same
events as the traditional matched filter method. At false alarm rates of
greater than $10^{-5}$, the SPIIR method recovers greater than 99% of the
injections recovered by the matched filter when searching for injections at an
effective distance of 250$\,\mathrm{Mpc}$ (SNR$\sim$7.9). Even in the worst
case, at a false alarm rate of $10^{-6}$, the SPIIR method catches 4.5% of
injections scaled at an extreme 400$\,\mathrm{Mpc}$ (SNR$\sim$5), whereas the
matched filter catches 5% of injections at this scale.

Figure 6: The receiver operating characteristics (ROC) of both the IIR filter
method and the traditional matched filter method. The x-axis shows the false
alarm rate, and the y-axis the detection rate. A one-to one relationship,
which is the worst case scenario, is shown by the boundary of the shaded area.
We show four different ROC curves, where the each curve represents the
detection rate as a function of false alarm rate for waveform injected at
effective distances of 250, 300, 350 and 400 $\,\mathrm{Mpc}$ (SNR $\sim$ 8,
6.6, 5.7 and 5 respectively).

##  V Summary and Discussion

The use of a bank of simple IIR filters for each template as opposed to the
matched filter method enables us get two extra processes for a minimal
additional cost. The first is that the individual IIR filter outputs can be
arranged into groups, such that their total summed output is roughly
independent and orthogonal to each other. This enables, with minimal extra
overhead, the calculation of a $\chi^{2}$ distributed statistic, giving a
secondary method of verification. We will demonstrate this in an upcoming
paper. The second natural advantage of using a parallel bank of single-pole
IIR filters is that they can easily be executed in parallel using multi-
threaded processors, such as graphics processing units (GPUs). Indeed, a side
study has shown that this is possible Liu (2011). This leads to the future
possibility that a single personal computer may be able to process the
detection of GWs.

A further way to reduce the computation of the IIR calculation is to split the
incoming data into differently down-sampled channels. The output of each IIR
filter in the bank is the correlation of a fixed frequency sinusoid and the
incoming data. For the sinusoids that have frequencies $<$124$\,\mathrm{Hz}$,
the incoming data need only be sampled at 256$\,\mathrm{Hz}$. The current
pipeline of LLOID uses a similar multi-channel down-sampling in their
detection pipeline. Their pipeline consists of the integration of the open-
source real-time multimedia handling software gstreamer and the LIGO Algorithm
Library (LAL) gst . This software library is an ideal platform to integrate
the SPIIR method. The total computation can also be further reduced by sharing
IIR filters (via interpolation) between different templates Luan (2011).

Although the design of the IIR filter so far only applies to chirping, post-
Newtonian approximation inspirals, we have performed preliminary tests using
more complicated combinations of single-pole IIR filters to replicate the
waveform of an inspiral with spin. If the amplitude/frequency beating of a
spinning inspiral waveform can be simulated by the linear addition of two
different non-spinning inspirals with different masses, then it can be
approximated by a linear addition of damped sinusoids. In this case, the SPIIR
method can produce the SNR for the beating waveform. There is also the
possibility of using higher order IIR filters, although designing the
coefficients can be very difficult.

##  VI Conclusion

We have demonstrated that the through the use of a parallel bank of single
pole IIR filters, it is possible to approximate the SNR derived from the
matched filter with greater than 99% overlap. The main advantage of our SPIIR
method is that it operates completely in the time domain, and in principle it
has zero latency (not taking into account whitening or computational time).
The SPIIR method recovers most of the injections the optimal matched filter
recovers.

We foresee that the use of IIR filters for time domain filtering of Advanced
LIGO will be ideal, as the waveforms will be much longer. The frequency domain
matched filter will take more time to calculate GW triggers, essentially
ruling out the possibility of triggering the detection the prompt optical
emission related to neutron star mergers (GRBs). We have shown that the use of
a parallel bank of IIR filters requires less computational cost, with minimal
detection rate loss, and most importantly can be calculated in the time-domain
with near zero latency.

##  VII Acknowledgements

We would like to thank Kipp Cannon, Drew Keppel and Chad Hanna for detailed
discussion on the design and implementation of low-latency detection
algorithms. This work was done in part during the LIGO Visiting Student
Researcher program, which was partially funded by the 2009 UWA Research
Collaboration Award. This research was supported by the Australian Research
Council. SH gratefully acknowledges the support of an Australian Postgraduate
Award.

## References

  * (1) URL http://www.ligo.org.
  * (2) URL http://www.virgo.infn.it/.
  * Smith and for the LIGO Scientific Collaboration (2009) J. R. Smith and for the LIGO Scientific Collaboration, Classical and Quantum Gravity 26, 114013 (2009), eprint 0902.0381.
  * Adv (2007) _Advanced LIGO Reference Design_ , LIGO Tech. Rep. M060056 (2007), URL http://www.ligo.caltech.edu/docs/M/M060056-08/M060056-08.pdf.
  * Abbott et al. (2009) B. P. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. S. Amin, S. B. Anderson, W. G. Anderson, M. A. Arain, et al. (LIGO Scientific Collaboration), Phys. Rev. D 79, 122001 (2009).
  * Fox et al. (2005) D. B. Fox, D. A. Frail, P. A. Price, S. R. Kulkarni, E. Berger, T. Piran, A. M. Soderberg, S. B. Cenko, P. B. Cameron, A. Gal-Yam, et al., Nature (London) 437, 845 (2005), eprint arXiv:astro-ph/0510110.
  * Nakar (2007) E. Nakar, Phys. Rep. 442, 166 (2007), eprint arXiv:astro-ph/0701748.
  * van Putten (2009) M. H. P. M. van Putten, ArXiv e-prints (2009), eprint 0905.3367.
  * Zhang and Mészáros (2004) B. Zhang and P. Mészáros, International Journal of Modern Physics A 19, 2385 (2004), eprint arXiv:astro-ph/0311321.
  * Wainstein and Zubakov (1962) L. A. Wainstein and V. D. Zubakov, _Extraction of Signals from Noise_ (Prentice-Hall, 1962).
  * Abbott et al. (2008) B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson, W. G. Anderson, M. Arain, et al. (The LIGO Scientific Collaboration, http://www.ligo.org), Phys. Rev. D 77, 062002 (2008).
  * Luan (2011) J. Luan, Phys. Rev. D (2011), (to be submitted).
  * Buskulic et al. (2010) D. Buskulic, Virgo Collaboration, and LIGO Scientific Collaboration, Classical and Quantum Gravity 27, 194013 (2010).
  * (14) URL https://www.lsc-group.phys.uwm.edu/daswg/projects/gstlal.html.
  * Cannon et al. (2010) K. Cannon, A. Chapman, C. Hanna, D. Keppel, A. C. Searle, and A. J. Weinstein, Phys. Rev. D 82, 044025 (2010).
  * Anderson et al. (2001) W. G. Anderson, P. R. Brady, J. D. Creighton, and É. É. Flanagan, Phys. Rev. D 63, 042003 (2001), eprint arXiv:gr-qc/0008066.
  * Abbott et al. (2004) B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen, R. Amin, S. B. Anderson, W. G. Anderson, M. Araya, H. Armandula, et al., Phys. Rev. D 69, 122001 (2004), eprint arXiv:gr-qc/0308069.
  * Blanchet et al. (1995) L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Phys. Rev. Lett. 74, 3515 (1995), eprint gr-qc/9501027.
  * Blanchet et al. (1996) L. Blanchet, T. Damour, and B. R. Iyer, Phys. Rev. D 54, 1860 (1996).
  * Fairhurst and Brady (2008) S. Fairhurst and P. Brady, Class. Quant. Grav. 25, 105002 (2008), eprint 0707.2410.
  * Brady and Fairhurst (2008) P. R. Brady and S. Fairhurst, Classical and Quantum Gravity 25, 105002 (2008), eprint 0707.2410.
  * Allen et al. (2005) B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. E. Creighton (2005), eprint gr-qc/0509116.
  * Droz et al. (1999) S. Droz, D. J. Knapp, E. Poisson, and B. J. Owen, Phys. Rev. D 59, 124016 (1999), eprint arXiv:gr-qc/9901076.
  * Liu (2011) Y. Liu, J. Comput. Phys. (2011), (in preparation).

##  Appendix A Noise Spectral Density

We use an algebraic expression for the noise spectral density of Advanced LIGO
detectors defined by,

$\displaystyle\begin{split}S_{h}(f)=S_{0}\left\\{\left(\frac{f}{f_{0}}\right)^{-4.14}-5\left(\frac{f_{0}}{f}\right)^{2}+\right.\\\
\left.111\left(\frac{1-\frac{f}{f_{0}}^{2}+0.5\frac{f}{f_{0}}^{4}}{1.+0.5\frac{f}{f_{0}}^{2}}\right)\right\\};\end{split}$
(42)

where, $f_{0}=215$Hz and $S_{0}=10^{49}$.