Table Of ContentDraft version January 23, 2015
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
NANOGRAV CONSTRAINTS ON GRAVITATIONAL WAVE BURSTS WITH MEMORY
Z. Arzoumanian1, A. Brazier2, S. Burke-Spolaor3,20, S. J. Chamberlin4, S. Chatterjee2, B. Christy5,
J. M. Cordes2, N. J. Cornish6, P. B. Demorest3, X. Deng7, T. Dolch2, J. A. Ellis8,21, R. D. Ferdman9,
E. Fonseca10, N. Garver-Daniels11, F. Jenet12, G. Jones13, V. M. Kaspi9, M. Koop7, M. T. Lam2, T. J. W. Lazio8,
L. Levin11, A. N. Lommen5, D. R. Lorimer11, J. Luo12, R. S. Lynch14, D. R. Madison2,†, M. A. McLaughlin11,
S. T. McWilliams11, D. J. Nice15, N. Palliyaguru11, T. T. Pennucci16, S. M. Ransom17, X. Siemens4, I. H. Stairs10,
D. R. Stinebring18, K. Stovall19, J. Swiggum11, M. Vallisneri8, R. van Haasteren8,21, Y. Wang12, & W. W. Zhu10.
(NANOGrav Collaboration(cid:63))
Draft version January 23, 2015
5
1
ABSTRACT
0
2 Among efforts to detect gravitational radiation, pulsar timing arrays are uniquely poised to de-
tect “memory” signatures, permanent perturbations in spacetime from highly energetic astrophysical
n
events such as mergers of supermassive black hole binaries. The North American Nanohertz Obser-
a
vatory for Gravitational Waves (NANOGrav) observes dozens of the most stable millisecond pulsars
J
using the Arecibo and Green Bank radio telescopes in an effort to study, among other things, grav-
1
itational wave memory. We herein present the results of a search for gravitational wave bursts with
2
memory (BWMs) using the first five years of NANOGrav observations. We develop original methods
for dramatically speeding up searches for BWM signals. In the directions of the sky where our sen-
]
A sitivity to BWMs is best, we would detect mergers of binaries with reduced masses of 109 M(cid:12) out
to distances of 30 Mpc; such massive mergers in the Virgo cluster would be marginally detectable.
G
We find no evidence for BWMs. However, with our non-detection, we set upper limits on the rate at
. which BWMs of various amplitudes could have occurred during the time spanned by our data–e.g.,
h
BWMs with amplitudes greater than 10−13 must occur at a rate less than 1.5 yr−1.
p
- Subject headings: pulsars, gravitational waves
o
r
t
s 1. INTRODUCTION
a
Due to the intrinsically nonlinear nature of Ein-
[ (cid:63)Authororderalphabeticalbysurname
stein’s equations, all systems that radiate gravitational
†Contactauthor: D.R.Madison,[email protected]
1 1Center for Research and Exploration in Space Science and waves(GWs)areanticipatedtoproduce“memory,”non-
v Technology and X-Ray Astrophysics Laboratory, NASA God- oscillatory components of the gravitational waveforms
3 dardSpaceFlightCenter,Code662,Greenbelt,MD20771,USA (Smarr 1977; Bontz & Price 1979; Braginskii & Thorne
4 2Department of Astronomy, Cornell University, Ithaca, NY 1987; Christodoulou 1991; Blanchet & Damour 1992).
14853,USA
3 3National Radio Astronomy Observatory, 1003 Lopezville Supermassive black hole binaries (SMBHBs), during the
5 Rd.,Socorro,NM87801,USA final few orbits preceding their merger, are expected to
0 4Center for Gravitation, Cosmology and Astrophysics, De- generate GW bursts with memory (BWMs) with suffi-
. partment of Physics, University of Wisconsin-Milwaukee, P.O.
1 Box413,Milwaukee,WI53201,USA ciently large amplitudes to make them potentially de-
0 5DepartmentofPhysicsandAstronomy,Franklin&Marshall tectable with pulsar timing arrays (PTAs; Favata 2009;
5 College,P.O.Box3003,Lancaster,PA17604,USA Seto 2009; Pshirkov et al. 2010; van Haasteren & Levin
1 6DepartmentofPhysics,MontanaStateUniversity,Bozeman, 2010;Cordes&Jenet2012;Madisonetal.2014). Cutler
MT59717,USA
v: 7Department of Astronomy and Astrophysics, Pennsylvania etal.(2014)recentlysingledoutmemoryasakeydetec-
i StateUniversity,UniversityPark,PA16802,USA tiontargetforPTAsasaprobeofexoticandunexpected
X 8Jet Propulsion Laboratory, California Institute of Technol- GW sources like phase transitions in the early Universe.
ogy,4800OakGroveDrive,Pasadena,CA91106,USA To facilitate GW detection, several international con-
r 9Department of Physics, McGill University, 3600 University
a sortiaarecurrentlyusingsensitiveradiotelescopespaired
St.,Montreal,QCH3A2T8,Canada
10Department of Physics and Astronomy, University of with pulsar-optimized hardware to realize the PTA con-
BritishColumbia,6224AgriculturalRoad,Vancouver,BCV6T cept (Hellings & Downs 1983; Foster & Backer 1990).
1Z1,Canada
The European Pulsar Timing Array (EPTA; Kramer &
11DepartmentofPhysicsandAstronomy,WestVirginiaUni-
versity,P.O.Box6315,Morgantown,WV26506,USA Champion2013),theNorthAmericanNanohertzObser-
12Center for Gravitational Wave Astronomy, University of vatory for Gravitational Waves (NANOGrav; McLaugh-
TexasatBrownsville,Brownsville,TX78520,USA lin 2013), and the Parkes Pulsar Timing Array (PPTA;
13Department of Physics, Columbia University, New York,
Hobbs 2013) are pushing precision pulsar timing to its
NY10027,USA
14NationalRadioAstronomyObservatory,P.O.Box2,Green limits and developing new data analysis techniques to
Bank,WV24944,USA usherintheeraofPTAGWastronomy. Bypoolingtheir
15Department of Physics, Lafayette College, Easton, PA data and expertise, these consortia have formed the In-
18042,USA
16UniversityofVirginia,DepartmentofAstronomy,P.O.Box ternational Pulsar Timing Array (IPTA; Manchester &
400325,Charlottesville,VA22904,USA IPTA2013)whichispoisedtobecomethemostsensitive
17National Radio Astronomy Observatory, 520 Edgemont of all PTAs.
Road,Charlottesville,VA22903,USA
2
In recent years, the various PTAs have begun to place by changes in dispersion measure (DM) to be accounted
astrophysically meaningful upper limits on continuous for in the timing model fit (Lam et al. 2014).
GWs from individually resolvable SMBHBs (Arzouma- Three of the pulsars comprising the five-year
nian et al. 2014; Zhu et al. 2014) and a stochastic back- NANOGrav data set, J1853+1308, J1910+1256, and
ground of GWs (van Haasteren et al. 2011; Demorest B1953+29, do not have sufficient dual-frequency cover-
et al. 2013; Shannon et al. 2013). Wang et al. (2015) age to correct for timing errors from variations in DM
have searched for BWMs in the first approximately six overtime; weexcludethesepulsarsfromouranalysisfor
years of PPTA data; they detected nothing, but deter- thisreason. WealsoexcludethedatasetforJ1600−3053
minedwith95%confidencethatBWMswithamplitudes fromouranalysisbecauseitiscomparativelyveryshort,
greater than 10−13 occur at a rate less than 0.8 yr−1. spanning just two years. Searching for BWMs in such
Inthispaper, wesearchforBWMsinthefirstapprox- shortdatasetsisfeasibleinprinciple,butweavoiditfor
imately five years of NANOGrav data using new tech- two reasons. First, the minimum detectable BWM am-
niquesthatleadtoconsiderablecomputationalspeedups plitude in a particular data set scales approximately as
over the search methods described by Madison et al. T−3/2 whereT isthespanofthedataset(vanHaasteren
(2014) and those used by Wang et al. (2015). In Sec- &Levin2010),sowedonotanticipatethatthisdataset
tion 2, we describe the data used in our BWM search. will greatly improve our sensitivity to BWMs. Second,
In Section 3, we describe the BWM signal model we use in such short data sets, many timing model parameters
forouranalysis. InSection4, wediscussmodelsofnoise are highly covariant with each other (e.g., spin and as-
and how our sensitivity to BWMs is influenced by them. trometricparameters;Madisonetal.2013)andwithany
In Sections 5 and 6, we describe our search techniques, BWM signal present.
differentiatingbetweensearchesforso-calledpulsar-term Finally we exclude the data for J1643−1224 from our
and Earth-term events. In Sections 7 and 8, we present analysis because we believe it contains chromatic tim-
the results of our pulsar-term and Earth-term searches, ing biases from unaccounted-for phenomenology in the
respectively. In Section 9, we place upper bounds on interstellar medium (ISM; see Figure 2). The DM of
BWM rates and amplitudes. In Section 10 we summa- J1643−1224isapproximately62pccm−3,nearlyafactor
rize our key results and offer concluding remarks. We oftwogreaterthananyotherpulsarinoursample. With
provideasummaryofourimportantnotationalelements high DM pulsars, chromatic timing errors that deviate
in a table of key symbols in an appendix. from the ν−2 scaling expected from cold plasma disper-
sion alone become more significant (Cordes & Shannon
2. PULSARTIMINGDATASET 2010). Furthermore,thispulsarisdirectlybehindacom-
In this section, we discuss several aspects of the data plex region of H-II associated with ζ-Ophiuchi, a mas-
that are relevant to our analysis. For a more thorough sive, runaway O-type star spinning very near breakup
descriptionofthedata,seeDemorest(2007)andDemor- (Gaustad et al. 2001; Villamariz & Herrero 2005). This
est et al. (2013). intervening H-II region could conceivably contribute to
Thefive-yeardatasetconsistsofpulsetimes-of-arrival non-trivial and currently unaccounted-for chromatic ef-
(TOAs)collectedatapproximatelymonthlyintervalsfor fects on the timing behavior of J1643−1224.
each of 17 millisecond pulsars (MSPs). All observations
weredonewiththeAreciboradiotelescopeandtheGreen 3. SIGNALMODEL
Bank Telescope (GBT). Of these 17 MSPs, those visi-
For a given pulsar, the timing perturbation from a
ble to Arecibo were observed with Arecibo; all others
BWM of amplitude h is well-modeled as
were observed with the GBT. One pulsar, J1713+0747, B
wasobservedwithbothtelescopes. Allobservationswere
∆t(t)=h B(θ,φ)[(t−t )Θ(t−t )−
doneusingoneoftwoidenticalbackendsystems: theAs- B 0 0
(t−t )Θ(t−t )]. (1)
tronomical Signal Processor (ASP) and the Green Bank 1 1
Astronomical Signal Processor (GASP). These backends
The function B(θ,φ) = (1/2)(1−cosθ)cos(2φ) ranges
performedreal-timecoherentdedispersionoverbandsup
between −1 and 1 and is common to all pulsar timing
to 64 MHz wide and recorded the results averaged over
efforts to detect point-like sources of GWs (Estabrook
channels of width 4 MHz each.
& Wahlquist 1975; Hellings & Downs 1983; Lee et al.
For each observing epoch, several TOAs are reported
2011). The angle between the direction the burst prop-
from various frequency channels of the 64 MHz band.
agates and the line of sight from Earth to the pulsar is
At Arecibo, observations were typically conducted at
θ; φ is the angle between the principal polarization vec-
two widely-separated frequencies (usually 430 MHz and
tor of the wave and the projection of the line of sight
1400 MHz) within one day; at the GBT, observations at
from the Earth to the pulsar onto the plane normal to
two frequencies (usually 820 MHz and 1400 MHz) were
the wave propagation direction. The BWM encounters
conducted within several days of each other. Approx-
the Earth at a time t and is observed from Earth to
imately contemporaneous observations at multiple fre- 0
encounter the pulsar at a time t = t +(l/c)(1+cosθ)
quenciesallowepoch-to-epochtimingfluctuationscaused 1 0
wherel isthedistancefromtheEarthtothepulsar(van
Haasteren & Levin 2010; Cordes & Jenet 2012; Madison
18 Department of Physics and Astronomy, Oberlin College, et al. 2014). The function Θ is the Heaviside step func-
Oberlin,OH44074,USA
19 Department of Physics and Astronomy, University of New tion. The amplitude of a BWM coming from a SMBHB
Mexico,Albuquerque,NM87131,USA mergerofreducedmassµ≡M1M2/(M1+M2)(M1 and
20 JanskyFellow M are the masses of the black holes in the binary) with
2
21 EinsteinFellow a typical inclination angle of I = π/3 at a luminosity
3
distance D from Earth is (Madison et al. 2014) that the radiometer noise adds to a fixed pulse shape
L
undertheassumptionsofmatchedfilteringagainstavery
(cid:18) (cid:19)(cid:18) (cid:19)
µ 1 Mpc
h ≈1.5×10−12 . (2) high S/N template pulse profile.
B 109 M D The second random contribution to t , (cid:15) , is from
(cid:12) L ν,τ Jν,τ
phase jitter in single pulses. Jitter is highly correlated
The distances to the pulsars in the NANOGrav array
between frequency channels, but is known to decorre-
are on the order of kiloparsecs. Our timing baseline is
late over widely-separated frequencies; the probability
approximately five years. Unless θ differs from π by less
distribution function of the phase of single-pulse cen-
than ∼ 3 degrees for a particular pulsar, we expect that
troids evolves with frequencies similar to pulse profiles
t andt willnotbothfallwithinourobservingwindow.
0 1 (Shannon et al. 2014). Jitter noise can be modeled as
Since we consider only 12 pulsars in our analysis, less
than 1% of the sky is within 3 degrees of one of our pul- (cid:104)(cid:15) (cid:15) (cid:105)=σ2(N)δ ρ ≈σ2(N)δ , (6)
sars. AssumingBWMsoccurisotropically,thereisaless Jν,τ Jν(cid:48),τ(cid:48) J ττ(cid:48) Jν,ν(cid:48) J ττ(cid:48)
than 1% chance of both t and t occurring within our whereσ ∝N−1/2. Thequantityρ isapproximately
five-year observing span if0 a BW1M occurs at all. And unity unJless ν and ν(cid:48) are very wideJlyν,νs(cid:48)eparated.
if a BWM occurred with a small enough angular sepa- Usually, σ (cid:46) σ . Jitter noise begins to dominate
J S/N
ration from a pulsar in our array that we could see the
radiometer noise only when the S/N of single pulses be-
timingperturbationbothturnonatt andturnoffatt ,
0 1 gins to exceed unity. For our data, collected with ASP
it would only be observed to turn off in that one pulsar.
and GASP, jitter noise should be significantly subdomi-
So, in each of our pulsar timing data sets, we need only
nanttoradiometernoiseforallpulsars. Bothradiometer
to look for evidence of timing perturbations of the form
andjitternoisearesignificantlylargerthancontributions
∆t(t)=hp(t−tB)Θ(t−tB), (3) from (cid:15)DISSν,τ, a random error associated with diffrac-
tive interstellar scintillation (DISS; Cordes 1990; Cordes
wherehp =±hBB(θ,φ),whatwecalltheprojectedburst & Shannon 2010). In a detailed study of J1713+0747,
amplitude,andtB iseithert0ort1,whatwecalltheburst Shannon & Cordes (2012) found that for single pulses,
epoch. Burstsarrivingatanindividualpulsaronlyinflu- σ ≈27±1µs,whiletimingerrorsfromDISSwereonlya
J
encethetimingbehaviorofthatpulsar;werefertothese fewnanoseconds. Arecent24-hourstudyofJ1713+0747
as pulsar-term bursts (see Section 5). Bursts arriving (Dolch et al. 2014) confirmed the jitter measurement of
at the Earth will simultaneously begin to influence the Shannon & Cordes (2012). Errors from DISS can be
timing behavior of all pulsars in our sample; we refer to substantially more important in the timing of high DM
these as Earth-term bursts (see Section 6). pulsars like B1937+21. We expect DISS to be a sub-
dominant component of our timing error budget for the
4. NOISEMODEL MSPs we include in our sample, so we ignore it in our
ConsiderTOAsmeasuredusingpulseprofilesobtained analysis. In summary, the noise covariance matrix for
by synchronously averaging N pulses for each of several a pulsar in our array, with these anticipated sources of
radio-frequency channels. The pulse profiles are func- noise, can be approximated as
tions of pulse phase, ϕ, channel frequency, ν, and ob-
(cid:104) (cid:105)
serving epoch, τ. A TOA from a particular frequency C =(cid:104)(cid:15) (cid:15) (cid:105)≈δ δ σ2 (N)+σ2(N) ,
νν(cid:48),ττ(cid:48) ν,τ ν(cid:48),τ(cid:48) ττ(cid:48) νν(cid:48) S/N J
channel and observing epoch can be written as
(7)
t =t +t +t +(cid:15) +(cid:15) +(cid:15) ,
ν,τ ∞,τ DMν,τ Cν,τ S/Nν,τ Jν,τ DISSν,τ where (cid:15)ν,τ is the sum of all noise influencing the TOA
(4) from observing epoch τ and frequency channel ν.
where t is the TOA at infinite frequency, t is the
∞,τ DM 4.1. Empirical Noise Models
dispersive delay from propagation through ionized inter-
stellar plasma, t is an additional, non-dispersive chro- In a previous analysis of the first five years of
C
matic perturbation, such as intrinsic profile evolution NANOGrav data, Arzoumanian et al. (2014) searched
with frequency, pulse broadening from multipath prop- for continuous GWs using a different noise model:
agation, and interstellar refraction (Cordes & Shannon (cid:104) (cid:105)
2010). If unaccounted for, tDM and tC can lead to sys- Cνν(cid:48),ττ(cid:48)=δττ(cid:48) δνν(cid:48)(Q2+E2σS2/N(N))+J2 +
tematic errors in TOA estimates. With multi-frequency
C (A,γ). (8)
observations at each observing epoch, effects from t red
DM
are mitigated in the NANOGrav data set. Potential ef- ThetermsQandE arecommonlyreferredtoasEQUAD
fectsfromt arecombattedbyfittingforconstantinter- (a source of Gaussian white noise with time units added
C
channel offsets unique to each pulsar. Unlike t and in quadrature to radiometer noise) and EFAC (a dimen-
DM
t , the ‘(cid:15)’ terms in Equation 4 are random errors from a sionless constant multiplier on the anticipated amount
C
variety of noise sources. of radiometer noise), respectively. The term J mimics
The term (cid:15) is uncorrelated between frequency the correlations induced between TOAs from different
S/N
channels and betwν,eτen observing epochs, i.e., frequency channels of the same observing epoch by jit-
ter, but qualitatively differs from jitter in that its value
(cid:104)(cid:15)S/Nν,τ(cid:15)S/Nν(cid:48),τ(cid:48)(cid:105)=σS2/N(N)δνν(cid:48)δττ(cid:48), (5) doesnotchangetoreflectthenumberofpulses,N,aver-
aged together to yield a TOA. Most observations in the
where σ ∝ (S/N)−1 ∝ N−1/2. The quantity σ five-year NANOGrav data set we considered had similar
S/N S/N
is the TOA uncertainty from radiometer noise assuming integration times, so N does not fluctuate substantially
4
between epochs and J is thus very much like a jitter-
induced correlation. Different values of Q, E, and J
were used for each widely separated observing frequency
of each pulsar.
The Arzoumanian et al. (2014) noise model also in-
cluded red noise with a power-law power spectrum
P(f) = A[f/(1 yr−1)]γ. This component of the noise
model is not to be thought of as part of the TOA er-
ror budget; TOAs measured with extreme accuracy may
still differ from the expectations of a timing model be-
cause of spin noise intrinsic to a pulsar, which is well
modeled as a stochastic process with a power-law spec-
trum (Shannon & Cordes 2010). In the Arzoumanian
etal.(2014)analysis,undertheassumptionthatnoGW
signal was present in the data set, the noise model pa-
rameters Ξ=[A,γ,E,Q,J] were determined by finding
the maximum of the likelihood function
exp(cid:2)−1(R−Mδp)TC−1(R−Mδp)(cid:3)
L(R|δp,Ξ)= 2 ,
(cid:112)
(2π)NTOAdetC
Fig. 1.— Uncertainty on the projected amplitude of a BWM
(9) at various trial burst epochs with the NANOGrav five-year data
set for PSR J1713+0747 using three different noise models. The
where R are the timing residuals from an initial timing lowest curve assumes that only radiometer noise is present in the
model,δparedeviationsofthetimingmodelparameters data and is identical to the results of Madison et al. (2014). The
from the initial timing model, M is the timing model middle curve assumes that only radiometer noise and jitter noise
are present in the data (as in Equation 7); the scale of the jitter
design matrix, and N is the total number of TOAs
TOA contribution is consistent with Shannon & Cordes (2012). The
in the data set. most conservative curve is based on the fixed noise model used
WeappliedthenoisemodelassessmentdonebyArzou- in Arzoumanian et al. (2014) (as in Equation 8) sans a red noise
manian et al. (2014) to two simulations of the five-year component.
NANOGravdatasetforPSRJ1909−3744,onethatcon- physically motivated noise model of Equation 7 and the
tainedonlysimulatedradiometernoiseandonethatcon- empiricalnoisemodelofEquation8isnotcurrentlywell
tained radiometer noise and a bright BWM occurring at understood. Detailed studies like Dolch et al. (2014) are
themidpointofthedataset. Theresultingbest-fitnoise being carried out to better understand the noise budget
modelswereidenticalexceptforintherednoiseparame- of pulsar timing experiments and any systematic effects
tersAandγ. TheBWMsignalisheavilycovariantwith thatinfluencethetimingprocedure,butbridgingthegap
the red noise component of the noise model, resulting in between the top and middle curves of Figure 1 is an on-
a large value for A and a poorly constrained value of γ. going area of research. The Arzoumanian et al. (2014)
An independent analysis of the five-year NANOGrav noisemodelisthemostconservativeofthethreewecon-
data set by Perrodin et al. (2013) found that except sider so we use it for the remainder of our analysis.
for J1643−1224 and J1910+1256, pulsars we already ex-
cludedfromouranalysis,allpulsarswereconsistentwith 4.3. Epoch Averaging
white noise alone. Because the demonstrated covariance
InFigure2,wedepicttheepoch-averagedtimingresid-
between red noise and BWM signatures would compli-
ualsfromthe12pulsarsinoursample(andtheexcluded
cate detection of a BWM (Cordes & Jenet 2012) and
pulsarJ1643−1224). Residualsareobtainedindividually
because of the lack of strong evidence for red noise, we
formultiplefrequencychannelsforeachintegrationlast-
have opted to consider the Arzoumanian et al. (2014)
ing15to45minutes. Theresidualsfromthatintegration
fixed noise models without the red noise component.
are combined to form a single epoch-averaged residual.
The noise model is central to this procedure. We define
4.2. Comparing Noise Models
an operator,
In Figure 1, we show how the uncertainty on the pro-
jected BWM amplitude, σh,p, varies in the NANOGrav A=(cid:0)UTC−1U(cid:1)−1UTC−1, (10)
five-year data set for PSR J1713+0747 over several trial
burst epochs under three different noise models. The which maps the raw residuals, R, to epoch-averaged
quantity σh,p is a direct measure of how sensitive a par- residuals, RE = AR. The matrix U is the “exploder”
ticular data set is to BWMs and it is directly influenced matrixdiscussedinArzoumanianetal.(2014)thatmaps
bythenoisemodel;itisdiscussedby,e.g.,vanHaasteren epochs to the full set of TOAs. The uncertainties on the
& Levin (2010) and Madison et al. (2014). epoch-averaged residuals are the square roots of the di-
ThebottomcurveofFigure1isbasedonanoisemodel agonal entries of the matrix
that only includes radiometer noise. The middle curve
C =(cid:104)R RT(cid:105)=(cid:0)UTC−1U(cid:1)−1. (11)
uses the noise model described by Equation 7 with the E E E
jittermeasurementfromShannon&Cordes(2012). The
top curve is based on the noise model used by Arzouma- 5. SEARCHESFORBWMSINTHEPULSARTERM
nian et al. (2014) and described in Equation 8 (without Informationinthepulsartermsofthe12pulsarsinour
the red noise component). The discrepancy between the samplecomesfromcausallydistinctregionsofspacetime
5
cult to rule out pulsar-specific phenomena. Nonetheless,
we carry out a search for BWMs in each of our individ-
ual pulsars as non-detections in many pulsars allow us
to place constraints on otherwise inaccessible regions of
BWM parameter space.
Madison et al. (2014) describe techniques for search-
ing for BWMs in individual pulsar timing data sets and
forassessingtheminimalprojectedamplitudedetectable
at a particular epoch. The timing perturbation from a
BWM is deterministic and can be included as part of
the timing model. The projected amplitude enters the
timing model as a linear parameter and can be fit for
using least-squares methods (Gregory 2010), yielding an
estimate for the projected amplitude, hˆ , and its uncer-
p
tainty,σ . Wedealwiththenonlinearparametert by
h,p B
searching over a grid of trial burst epochs. For all of our
timing model fits and calculation of timing residuals, we
use the software package TEMPO2 (Edwards et al. 2006).
The modification to the timing solution caused by in-
cludingaburstofprojectedamplitudeh attimet can
p B
beassessedbycomputingthelikelihoodratioforamodel
with and without a burst:
(cid:18) (cid:19)
Γ(h ,t )=exp −1(cid:2)χ2(h ,t )−χ2(0,t )(cid:3) .(12)
p B 2 p B B
In this expression,
χ2(h ,t )=RT(h ,t )C−1R(h ,t ), (13)
p B p B p B
whereR(h ,t )arethetimingresidualswhenaburstof
p B
projected amplitude h at time t is included as part
p B
of the timing model. The least-squares estimator for
the projected burst amplitude at a trial burst epoch t ,
B
hˆ (t ), maximizes Γ at that trial burst epoch; call this
p B
valueΓˆ(t ). Forafixedtrialburstepoch,ifthedataare
B
consistent with the noise model, Dˆ ≡ 2lnΓˆ, the reduc-
tion in the χ2 of the residuals caused by introducing a
projectedburstamplitudetothetimingmodelfitwillbe
a random variable following a χ2 distribution with one
degree of freedom:
Fig. 2.—Epoch-averagedtimingresidualsforthe12pulsarsused
inouranalysis(andtheexcludedpulsarJ1643−1224). Thevarious f (Dˆ)=(2πDˆ)−1/2exp(−Dˆ/2). (14)
colors indicate different observing frequencies: 327 MHz is green, 1
430 MHz is blue, 820 MHz is red, 1400 MHz is black, and 2300
Since the burst amplitude is a linear parameter in the
MHz is cyan. We show the residuals for J1643−1224 to illustrate
theapparentchromaticproblemswiththedata: theresidualsfrom timing model, we expect the χ2 value of the residuals to
820 MHz and 1400 MHz are consistently anti-correlated and are respond quadratically to the burst amplitude, i.e.,
oftenextremeoutliers.
(cid:16) (cid:17) (cid:104) (cid:105)2
χ2(0,t )=χ2 hˆ (t ),t + hˆ (t )/σ (t ) ,
(Cordes & Jenet 2012). A BWM is 12 times more likely B p B B p B h,p B
to encounter a single pulsar in our array than it is to
(15)
encountertheEarth. However,itisnottruethatBWMs
encounteringindividualpulsarsinthearrayare12times or Dˆ(t ) = [hˆ (t )/σ (t )]2. The σ values are the
B p B h,p B h,p
more likely to be confidently detected. Different pulsars
least-squares1-σ amplitudeuncertaintiesonhˆ ,e.g.,the
aretimedwithvaryingdegreesofprecisionandwithdif- p
quantities plotted in Figure 1.
fering observing cadences making them unequally sensi-
tive in searches for BWMs. A particular BWM may be IfwecomputeDˆ fortwotrialepochsthatareveryclose
polarized in such a way or from such a part of the sky together, the results will be correlated. Because of these
that the projection factor, B(θ,φ), makes its influence correlations,whencomputingDˆ alongadenselysampled
onthetimingbehaviorofaparticularpulsarvanishingly grid of N trial burst epochs, the number of effectively
t
small. Furthermore, non-BWM phenomena such as in- independent trial burst epochs tested is N < N . The
I t
trinsic pulsar spin noise (Shannon & Cordes 2010) and probability distribution function for the maximum value
microglitches (Cognard & Backer 2004) can be confused of Dˆ along the grid, Dˆ , is
max
asBWMs. Withoutthesignalappearingconcurrentlyin
multiple pulsars as in an Earth-term BWM, it is diffi- f (Dˆ )=N f (Dˆ )F (Dˆ )(NI−1), (16)
NI max I 1 max 1 max
6
epochs to test. With N fixed at five, inverting Equa-
I
tion 19 allows us to compute Dˆ for any desired al-
thresh
lowable false alarm probability; false alarm probabilities
of approximately five and one percent are expected for
Dˆ equal to 6.60 and 9.54, respectively.
thresh
6. SEARCHESFORBWMSINTHEEARTHTERM
SearchesforBWMsintheEarthtermhaveseveralad-
vantages over searches in pulsar terms. The timing per-
turbationfromaBWMwillturnonsimultaneouslyforall
pulsarsinthePTAifitoccursintheEarthterm,provid-
ingapowerfulmeansbywhichpulsar-specificphenomena
(e.g., glitches) can be ruled out. Variations in the pro-
jected BWM amplitude from pulsar to pulsar provide
informationaboutthelocationoftheGWsourceandits
polarization. As such, for a trial BWM source location,
polarization,andepoch,theresidualsofallpulsarsinthe
PTA can be combined into a coherent global fit for the
true burst amplitude, h . Such global fits, described in
Fig. 3.— Cumulative distributions of Dˆmax values from 1000 Madison et al. (2014), hBave better amplitude sensitivity
simulationsofnoise-liketimingresidualsforthreeofthepulsarsin
thanwhatisattainablewithanyonetimingdatasetand
oursample. Overlaidaretheanticipatedcumulativedistributionsif
NI,theeffectivenumberoftrialburstepochstested,isone,five,or canbecarriedoutwithasimilarleast-squaresapparatus
twenty. Forallpulsarsinoursample,theanalyticallyanticipated as is used in pulsar term searches.
cumulativedistributionbestfitstheresultsofoursimulationswith Global least-squares fitting techniques were recently
NI between4and5. used by the PPTA to search for BWMs (Wang et al.
where F is the cumulative distribution of f ,
1 1 2015). The PPTA data set those authors analyzed con-
(cid:18)(cid:113) (cid:19) tained many fewer TOAs than the NANOGrav data set
F1(Dˆ)=erf Dˆ/2 . (17) we are considering here owing to NANOGrav’s practice
of reporting many TOAs from a single observing epoch
The cumulative distribution for Dˆ associated with butfromdifferentfrequencychannels. Withthesemulti-
max
frequencyTOAs,NANOGravincludesaspartofitstim-
f is then simply
NI ing models numerous chromatic parameters to account
(cid:18)(cid:113) (cid:19) forprofileevolutionacrossourwide-bandwidthreceivers
F (Dˆ )=erfNI Dˆ /2 . (18)
NI max max and time-variable DM, so, we fit for many more timing
model parameters than the PPTA.
Theanticipatedfalsealarmprobabilityfornoisealoneto Suppose N is the number of timing model parame-
exceed a threshold value of Dˆ , Dˆ , is just ters being fit for each of M pulsars in a PTA (in reality,
max thresh
N varies from pulsar to pulsar). Then carrying out a
FNI(Dˆthresh)=1−FNI(Dˆthresh). (19) global fit for the amplitude of a BWM requires the in-
version of an (NM +1)×(NM +1) matrix, a proce-
For a fixed allowable false alarm probability, Dˆ
thresh durethatrequirescomputationtime∝(NM+1)3. This
grows logarithmically with N if N (cid:29)1.
I I globalfithastobedoneoverN N N trialsinthefour-
To estimate N , for each pulsar, we generated 1000 Ω ψ t
I dimensionalphasespaceofskyposition,burstepoch,and
simulatedsetsofTOAsthatmatchedtherealdatasetin
burstpolarization(N ,N ,andN are,respectively,the
numberofTOAs, observingschedule, andtimingmodel, Ω ψ t
number of sky positions, polarization angles, and burst
but yielded timing residuals consistent with our noise
epochs tested). Methods requiring these global fits are
models. We then fit for hˆ (t ) along an equispaced grid
p B thuscomputationallyonerousforNANOGravbecauseof
of twenty trial burst epochs between MJD 53500 and
the large number of timing model parameters for which
54900 and computed Dˆ(tB) to get Dˆmax. we fit. More importantly, the total number of suitably
In Figure 3, we show the cumulative distribution of stableMSPsbeingtimedbyNANOGravisgrowingwith
the1000Dˆ valuesfromoursimulationsforthreepul- time as projects like the Arecibo PALFA survey (e.g.,
max
sars: J0030+0451andJ2145−0750,observedbyArecibo Cordes et al. 2006; Swiggum et al. 2014), the Arecibo
and the GBT, respectively, with rms timing residuals of All-sky327-MHzDriftPulsarSurvey(e.g.,Denevaetal.
∼ 100 ns, and J1713+0747, observed by both Arecibo 2013), and the GBT Northern Celestial Cap Pulsar Sur-
and the GBT with an rms timing residual of ∼ 30 ns. vey (e.g., Stovall et al. 2014) continue to find more pul-
We also plot the theoretically anticipated curves for N sars. A search procedure based on global fits will not
I
equal to one, five, and twenty. In fitting the anticipated scale well to future PTA endeavors. We use a different
cumulative distributions of Dˆ with N as a free pa- technique that avoids having to carry out any global fits
max I
rametertotheresultsofoursimulations,wefindthatfor and is computationally much more efficient.
allpulsars,N isbetweenfourandfive. AsN increases,
I I
large values of Dˆ occur more frequently in the pres- 6.1. Accelerated Earth-term Searches
max
ence of pure noise. We take the conservative approach Over a five-dimensional grid of trial burst sky posi-
and assume that there are five independent trial burst tions, Ω , polarizations, ψ , epochs, t , and ampli-
ij k B,l
7
tudes, h , we search for the maximum in the global likelihood ratio surfaces will be faster than global fit-
B,m
likelihood ratio surface, ting schemes. In practice, we have set N = 300, test-
p
ing150trialprojectedamplitudeslogarithmicallyspaced
ΓG(Ωij,ψk,tB,l,hB,m) (20) between 5×10−17 and 10−12 and their negatives. Uti-
M lizing the precomputed two-dimensional likelihood ratio
(cid:89)
= Γ [B (Ω ,ψ )h ,t ], surfaces is thus faster than using global fitting by a fac-
K K ij k B,m B,l
tor of over thirteen thousand and is the only reason why
K=1
searching such a densely sampled grid is feasible.
where BK is the projection factor, B(θ,φ), for the Kth Once the global, five-dimensional likelihood ratio sur-
pulsar. TheΓK quantitiesareanalogoustothelikelihood face is computed, for fixed Ωij, ψk, and tB,l, we isolate
ratiosinEquation12. Theanglesθ andφdependonthe thevalueofh thatmaximizesΓ ;wecallith˜ . This
location of the Kth pulsar and the trial burst location, B,m G B
four-dimensional surface, which we call Γ˜ , is approxi-
Ω . The trial burst polarization angle, ψ , influences φ. G
Thijetotalnumberoftrialsinthisfive-dimkensionalsearch matelyequaltoΓˆG,whatwewouldgetifwehadcarried
is N N N N , where N is the number of trial burst out the computationally costly global fits we have taken
Ω ψ t h h
amplitudes tested. great care to avoid, differing from Γˆ only because of
G
Ostensibly, the procedure appears to require that for the finite resolution of our trial burst amplitude grid.
each grid point in the five-dimensional space and for We can then, as in the pulsar-term case, consider the
each pulsar in our sample, we incorporate the signa- false alarm statistics of the quantity D˜ = 2lnΓ˜ . For
G G
tureofaBWMofprojectedamplitudeBK(Ωij,ψk)hB,m any fixed grid point in our four-dimensional parameter
occurring at time tB,l into the timing model of the space, if the data are consistent with our noise models,
KThtihs wpuolusladr,tarkeefitcotmheputtimatiinogn tmimodee∝l, NanΩdNcψoNmtpNuhtMe ΓNK3., dD˜oGm.agWainithfonllooiwses-lχik2esdtaattais,tiicfswweitchonosnideerdetghreeeproofbafrbeiel--
a speed up over the global fitting scheme so long as ity distribution of D˜ , the maximum value of D˜
N < (NM +1)3/MN3. The average N for our sam- G,max G
h over the whole grid, also follows Equation 16, but with
ple of M =12 pulsars is approximately 90, meaning N
h a larger number of effective independent trials sampled
needs to be less than approximately 140 for this tech-
than in the pulsar-term case owing to a search not just
nique to be faster.
over time, but also over polarization angle and sky posi-
However, thereisagreaterspeed-uptobehad. Varia-
tion. We call the effective number of independent trials
tionsinΩ ,ψ ,andh onlyaltertheprojectedBWM
ij k B,m in the global search N . We discuss N in more detail
amplitude along different pulsar lines of sight. We can G G
in Section 8.
precomputethetwo-dimensional(projectedburstampli-
tude and burst epoch) likelihood ratio surface for each 6.2. Assessing BWM Amplitude Uncertainty
pulsar in our sample, a procedure requiring computa- A similar relationship exists between D˜ , h˜ , and σ
tion time ∝N N MN3, where N is the number of trial G B h
p t p as exists between Dˆ, hˆ and σ as discussed in Section
projected burst amplitudes considered. We then con- p h,p
struct the full five-dimensional likelihood ratio surface 5, i.e., Dˆ = (hˆ /σ )2. Since Γ is just the product
p h,p G
from Equation 20 by figuring out what the projected of all the pulsar-specific likelihood ratio surfaces (with
burstamplitudeintheKthpulsarwouldbeforapartic- appropriate projection factors, as in Equation 20),
ular choice of Ω , ψ , and h , looking to the pre-
ij k B,m M M
computed two-dimensional likelihood ratio surface for D˜ = (cid:88) D (B h˜ )=h˜2 (cid:88)(B /σ )2. (21)
the Kth pulsar, and interpolating between the nearest G K K B B K h,p,K
projectedburstamplitudestestedandtheprojectedam- K=1 K=1
plitude of interest from the five-dimensional search; this So, we have
step requires no additional timing model fits (the com-
D˜ =(h˜ /σ )2, (22)
putationally costly step), is very fast, and is essentially G B h
an exercise in efficiently searching lookup tables. where
Using precomputed two-dimensional likelihood ratio
surfacesforaglobalEarth-termBWMsearchaswehave (cid:34)(cid:88)M (cid:35)−1/2
just described will lead to a speed up over the global σh = (BK/σh,p,K)2 . (23)
fitting scheme if Np < NΩNψ(NM + 1)3/MN3. For K=1
our search, we have tested Nt = 40 trial burst epochs Equation 23 is equivalent to Equation 32 from van
evenly spaced between MJDs 53541 and 54995, Nψ =17 Haasteren & Levin (2010). While those authors ana-
trial burst polarization angles evenly spaced between 0 lyzed idealized pulsar timing data sets with equispaced
and π, and NΩ = 1598 trial burst sky positions isotrop- observing epochs, uniform TOA uncertainties, and very
ically distributed on the sky. As a brief but important simple timing models, their result holds more generally
aside,wehavechosentosearchwithintheparticularwin- and applies here.
dow of dates just mentioned because for each pulsar in
our sample, there is at least one collection of multifre- 6.3. Cross-checks with Bayesian Methods
quency observations before and after this window, with As an independent check, we also carry out the
the limiting pulsar on the early side being J0613−0200 Bayesian method of van Haasteren & Levin (2010), im-
and J2145−0750 on the late side. With the grid pa- plemented in the software package Piccard21, developed
rameters we have chosen, so long as N is less than ap-
p
proximately4×106,usingprecomputedtwo-dimensional 21 https://github.com/vhaasteren/piccard
8
TABLE 1
Summary of Results from Pulsar Term BWM Search
PSR hˆmpax/10−13 tmB,aix Dˆmax log10(F5)
(MJD)
J0030+0451 −0.87±0.77 54100 1.2 −0.11
J0613−0200 0.68±0.50 54584 1.8 −0.20
J1012+5307 1.47±0.79 54062 3.4 −0.54
J1455−3330 −2.57±3.35 54211 0.5 −0.02
J1640+2224 −2.65±1.45 53801 3.3 −0.52
J1713+0747 0.08±0.07 54808 1.4 −0.13
J1744−1134 −2.06±1.19 53541 2.9 −0.45
B1855+09 1.00±0.73 53578 1.8 −0.21
J1909−3744 −1.29±0.74 54994 3.0 −0.45
J1918−0642 −8.29±3.39 54994 5.9 −1.15
J2145−0750 26.60±11.51 54994 5.3 −0.99
J2317+1439 −2.40±1.20 53541 3.9 −0.67
Note. — Description of most significant BWM-like signal de-
tectedinthepulsartermofeachof12NANOGravdatasets. From
lefttoright,thecolumnsarethenameofthepulsar,themostsig-
nificant projected BWM amplitude calculated from least-squares
fitting, the trial burst epoch at which the most significant ampli-
tudewasfound,thecorrespondingDˆmax value,andthelogarithm
of the anticipated false alarm probability for that value of Dˆmax
assuming NI, the number of effectively independent trial burst
epochstested,isfive.
independently of the Madison et al. (2014) method. In
brief,thisBayesianmethodtakesthelikelihoodofEqua-
tion 9 as a function of the noise parameters Ξ, the tim-
ing model parameters δp, and the BWM parameters
(Ω ,ψ ,t ,h ). Similar to what we do for our fre-
ij k B,l B,m
quentist search, we keep the noise parameters fixed to
the values obtained in Arzoumanian et al. (2014). Our Fig. 4.— Results of a search for a BWM in the NANOGrav
prior distributions are flat in all stated BWM parame- five-year data set for PSR J1713+0747, the single most sensitive
ters, where we note that the prior in Ω is taken flat datasetinoursampleforsuchsearches. Thebottompanelshows
ij
over the sphere. We have been able to confirm that the theepoch-averagedtimingresidualsforJ1713+0747astheyarein
resultswedescribeinthispaperareidenticalforthefre- Figure2. Thesecondpanelfromthebottomshowshˆp±σh,pforthe
40trialburstepochswetested. Thethirdpanelfromthebottom
quentist and Bayesian methods. showsthedetectionstatisticDˆ varyingoverthespanoftrialburst
epochs tested. The remaining panels show a “heat map” of the
7. PULSARTERMRESULTS two-dimensional likelihood ratio surface, as in Equation 12, that
weuseforourEarth-termanalysis.
For all twelve of the pulsar data sets we analyze, our
search for pulsar-term BWMs occurring between MJDs
53541 and 54994 yields results that are entirely consis- datasetareslightoutliers,allowingforaninstantaneous
tent with our noise models. In Table 1, we summarize change in the spin period of the pulsar near the begin-
thekeyresultsinoursearchforpulsar-termBWMs. For ning or end of the data set can bring them more in line
each pulsar in our sample, we list the pulsar name, the withzeroandleadtoamodestreductionoftheχ2 value
most significant value of hˆ (according to the Dˆ num- of the residuals.
p Theindividualdatasetinoursamplethatmosttightly
ber) along with its 1-σ uncertainty, the trial burst epoch
constrainsBWMsisforJ1713+0747. Thisisinlinewith
at which this most significant amplitude was found, the
theexpectationsofMadisonetal.(2014). InFigure4,we
Dˆ value for that pulsar, and the logarithm of the
max show the detailed BWM search results for J1713+0747.
false alarm probability anticipated from noise alone for
Thebottompanelissimplytheepoch-averagedresiduals
that value of Dˆmax assuming NI = 5. The most sig- asinFigure2;weplotthemagainheretoexplicitlyshow
nificant event we find is in the data for J1918−0642; it howtheintervaloftrialburstepochstestedoverlapswith
is consistent with approximately seven percent of noise the TOA coverage. The second panel from the bottom
realizations. shows hˆ ±σ at the 40 trial burst epochs we tested.
The epochs at which the most significant projected p h,p
Thebest-fitprojectedamplitudeiscompletelyconsistent
BWMamplitudesoccuraredistributednearlyuniformly
with zero at each trial burst epoch considered. Further
throughout the range of trial epochs we tested with re-
echoingthis,inthethirdpanelfromthebottom,weshow
peat values occurring only at the very first and last
the value of Dˆ as it varies with trial burst epoch. Never
epochs tested. Clustering at the edges of the window
of tested epochs is not surprising. If there are few tim- does Dˆ approach the values necessary to be inconsistent
ing residuals outside of the window of tested trial burst with noise at the five- or one-percent level. In the re-
epochs, when testing the first or last trial epoch, the maining panels of Figure 4, we show a “heat map” of
pre-orpost-bursttimingmodelisconstrainedbyasmall the two-dimensional likelihood ratio surface, as in Equa-
numberofdatapoints. Iftheresidualsattheedgeofthe tion 12, that we use in our Earth-term search.
9
cation of J1713+0747.
The single largest value of D˜ we find is D˜ =
G G,max
6.03, which corresponds to h˜ =2.46σ . We mentioned
G h
at the end of Section 6.1 that the effective number of
trials tested in an Earth-term BWM search, N , will
G
be larger than N = 5 because of the search over not
I
just many trial burst epochs, but also over many trial
burst source locations and polarizations, implying that
for a fixed allowable false-alarm probability we have to
raisethethresholdweimposeonD˜ abovethethreshold
G
used on the detection statistic in a pulsar-term search.
However, D˜ is small enough that even if N were
G,max G
only five, this result would be entirely consistent with
more than 95% of realizations of our noise.
IntheBWManalysiscarriedoutbyWangetal.(2015),
when faced with marginally high values of their test
Fig. 5.—AtimeandpolarizationsliceofourEarth-termsearch.
statistic, they conducted an extensive suite of simula-
The plotted quantity, D˜1/2, is equivalent to the best-fit BWM
G tions to assess N , comparable to what we have done
amplitude, h˜B, in units of the 1-σ uncertainty on the amplitude, to assess N (andGdepicted in Figure 3), and were able
σ . Atthelocationindicatedbytheblackcircle,σ =6.07×10−15, I
thhesmallestvalueofσ inourentiresearch. Inohurwholesearch, to justifiably raise the detection threshold on their test
h
themaximumvalueofD˜1/2wefindisapproximately2.46,entirely statisticandruleoutadetection. Wefindnocomparably
consistent with our noisGe models at the 95% level even if NG is large values of D˜G that exceed our detection threshold
no greater than 5 as we used in our pulsar-term searches. The even if we underestimate N as five.
G
diamonds indicate the positions of the 12 pulsars in our analysis.
WedostillwantanestimateofN asitwillallowusto
ThelargestdiamondrepresentsJ1713+0747. G
applyupperlimitstothepopulationofBWMsgivenour
non-detection. Cornish & van Haasteren (2014) recently
demonstratedthattheresponseofaPTAtoGWsofany
waveform can be decomposed into a linear combination
of a finite number of modes, or sky maps. The number
ofmodesrequiredisequaltotwicethenumberofpulsars
in the array (the factor of two accounts for the two pos-
sible polarization modes of GWs). We use their result
to estimate that the number of statistically independent
samples in our Earth-term search over source-position
andpolarizationspaceis24,twicethenumberofpulsars
used in our analysis. Given five independent samples in
time, we thus adopt N =120.
G
Adopting N = 120 is a conservative estimate as our
G
sensitivity to BWMs is so strongly dominated by a sin-
glepulsar,sotheeffectivenumberofpulsarsinourarray
is fewer than 12. As mentioned following Equation 19,
for a fixed allowable false alarm probability, if N (cid:29) 1,
Fig. 6.—ThemaximumluminositydistanceofaSMBHBmerger D˜ onlydivergeslogarithmicallywithN ,sGoover-
causing a BWM detectable with 95% confidence given our sensi- G,thresh G
tivityaveraged overburst epoch andpolarizationangle. We have estimates of NG are not exceedingly deleterious for the
assumed the binary had a typical inclination angle of π/3 and a purpose of setting upper limits. With N =120, setting
G
reducedmassof109 M(cid:12);theplotteddistancesscalelinearlywith D˜ =12.4 assures a false alarm probability of less
this fiducial reduced mass. The diamonds indicate the positions G,thresh
of the 12 pulsars in our analysis. The largest diamond represents thanfivepercent. Thisisequivalenttorequiringthath˜
B
J1713+0747. Thegreentrianglerepresentsthepositionofthecen- be greater than 3.52σ in order for it to be inconsistent
h
ter of the Virgo Cluster. Just 16.5 Mpc from Earth, the Virgo
with 95% of realizations of noise. Again, in our Earth-
Cluster is near the edge of the volume in which we can detect
BWMswiththeseproperties. termsearch, we findnosignalthatmeetsor exceedsthis
level of significance.
8. EARTHTERMRESULTS In Figure 6, we have averaged σ over trial burst po-
h
In Figure 5, we show a two-dimensional slice of our larization angles and epochs and shown the maximum
Earth-term search results. For fixed trial burst polar- luminosity distance at which a SMBHB merger with an
ization angle and epoch, we depict the quantity D˜1/2 as inclinationangleofπ/3andµ=109M(cid:12)(consistentwith
G Equation 2) would be detectable with our data set with
it varies with trial burst source location. The quantity
D˜1/2 is the best-fit BWM amplitude, h˜ , in units of the 95% confidence, or where h˜B = 3.52(cid:104)σh(cid:105)ψ,t. Our sen-
G B sitivity is worst near the position of PSR J0613−0200.
uncertaintyontheamplitude,σ (seeEquation22). The
h This has little to do with J0613−0200, but is instead
slice we show contains the smallest value of σ found in
h because these sky positions are antipodal to our great-
our entire search, i.e., the point in our parameter space
estconcentrationofpulsars,especiallyourmostsensitive
where we are most sensitive to a BWM; its location on
pulsar, J1713+0747. The green triangle in Figure 6 in-
the sky is indicated by the black circle very near the lo-
10
dicates the position of the Virgo Cluster. Just 16.5 Mpc
from Earth, the Virgo Cluster falls very near the edge of
the volume over which we are sensitive to BWMs from
µ=109 M , I =π/3 binary mergers.
(cid:12)
9. BWMRATE-AMPLITUDECONSTRAINTS
To synthesize the results from both our Earth- and
pulsar-term analyses, we derive constraints on Λ(>h ),
B
the annual rate of BWMs from any part of the sky with
any polarization having amplitudes greater than h , as-
B
suming they occur as a Poisson process. This quantity
is readily relatable to astrophysical models of processes
producing BWMs, i.e., Equation 15 of Cordes & Jenet
(2012),Equation21ofMadisonetal.(2014),orEquation
17 of Wang et al. (2015).
Toward this end, define the quantities,
∆t∆Ω∆ψ (cid:88)Nt (cid:88)NΩ (cid:88)Nψ
τ (h ,n)= Θ(h −nσ ),
E B 4π2 B h,ijk
i j k Fig. 7.— 95% confidence Earth-term upper bound on the rate,
(24) Λ(>hB),ofBWMsoccurringwithamplitudesatoraboveampli-
tudeshB. ThetwocurvescomefromourEarth-termandpulsar-
and term analyses. The pulsar-term probes lower-rate events at high
amplitudesbecauseindividualpulsartermscontaincausallyinde-
τ (h ,n)= ∆t (cid:88)M (cid:88)Nt (cid:90) Θ(cid:18)h − nσh,p,K(cid:19)dψdΩ, ppuenlsdaern-tteirnmforcmuravteiosna.tTlohwenBeWarMtotaamlcpolnitvuedregsendceemoofntshteraEtaersthth-aatnda
P B 4π2 B B(θ,φ) singlepulsar,J1713+0747,isdominatingoursensitivity.
K i
Sincetherearefewereffectivelyindependenttrialsina
(25)
pulsar-termsearch,weshouldusealowervalueofnfora
where ∆t, ∆Ω, and ∆ψ describe the grid spacing in the 95% confidence upper limit. But, by treating the Earth-
Earth-term search. The quantities τ and τ are the and pulsar-terms comparably, we are able to highlight
E P
totaltimethatthePTAhadn-σsensitivitytoaBWMof an important fact about our data set. The near-total
amplitudeatleasth intheEarth-termandpulsar-term, convergence of our Earth- and pulsar-term constraints
B
respectively,weightedbythefractionofthetotalsource- for low values of h indicates that our upper limits are
B
location and polarization angle space over which that almost entirely dominated by a single pulsar data set:
sensitivity was achieved. Our definitions for τ and τ J1713+0747. This was anticipated by Madison et al.
E P
differ from nearly identical definitions in Madison et al. (2014) and holds true even with the more sophisticated
(2014) by an overall factor of two. In that work, it was noise modeling we have employed in this analysis.
mistakenly assumed that only BWM polarization angles
between0andπ/2mustbeconsidered. Weherecorrectly 10. CONCLUSION
carryoutasearchthattestspolarizationanglesbetween In this paper, we have conducted a search for BWMs
0 and π. in the first five years of NANOGrav data. We did not
With the definitions for τE and τP in place, we can detect any BWMs. Based on our current understanding
derive constraints on Λ(>hB) from our Earth-term and ofSMBHBs,themostconventionalanticipatedsourceof
pulsar-term analyses: bright BWMs, it is unsurprising that we did not. Wang
et al. (2015) predict that BWMs from SMBHB mergers
ln(1−Q)
Λ(>hB)<− τ , (26) exceeding amplitudes of 10−14 occur at a rate of just a
i few every 105 yr. However, Cordes & Jenet (2012) con-
where Q is the probability that at least one BWM oc- cludethatbecauseoflargeuncertaintiesinthingssuchas
curring at rate Λ(>h ) encounters the PTA during the theinspiralrateofSMBHBs,theactualrateofBWMsis
B
time τ and i is a placeholder for either E or P. essentially unconstrained. Also, we stress that memory
i
In Figure 7, we show our constraints on Λ(> h ) is a very general feature of bright GW events and large-
B
from our Earth- and pulsar-term analyses. We have set amplitude BWMs may be produced by wholly unex-
Q=0.95. We have set n=3.52 for both our Earth- and pected phenomena occurring at unconstrained rates; on-
pulsar-termconstraints. Thisnumbercomesfromrequir- goingsearchesforBWMsarecrucialforutilizingtheraw
ing a false alarm probability of less than five-percent in discovery potential of PTA observations (Cutler et al.
an Earth-term search with N = 120. Our estimate of 2014). The methods developed in this paper are read-
G
N is likely significantly larger than it needs to be. A ily generalizable to future, more sensitive data sets and
G
detailedsuiteofsimulationscouldgiveusamorerealistic the accelerated search techniques we have developed for
assessment of N which would lead to more constrain- Earth-term analyses will greatly expedite future BWM
G
ing upper limits, but for a fixed allowable false alarm searches.
probability, the amplitude a BWM must exceed to be Our Figure 7 and Figure 10 from Wang et al. (2015)
inconsistent with the noise scales as the square root of show constraints on the rate of BWMs from initial
the natural logarithm of N , so the improvement to the NANOGrav and PPTA data releases, respectively. The
G
upper limit from this analysis would be very slight. NANOGrav constraints and the PPTA constraints are
