OBSERVATIONAL SELECTION EFFECTS WITH GROUND-BASED GRAVITATIONAL WAVE DETECTORS

The detection of binary black holes with the Laser Interferometer Gravitational-wave Observatory (LIGO; Collaboration et al. 2015; Abbott et al. 2016a, 2016e) has ushered in the age of gravitational wave astronomy. A particularly promising avenue for exploring new physics is multi-messenger astronomy, involving the joint detection of GW sources, electromagnetic (EM) signals, or astrophysical particles (Kulkarni et al. 1998; Gehrels et al. 2005). During the last years of initial LIGO and Virgo, as well as in advanced LIGO's first observing run (O1; The LIGO Scientific Collaboration et al. 2016a), a large consortium of EM observers followed up GW candidates (Abadie et al. 2012a; Aasi et al. 2014; Abbott et al. 2016f), and there has been substantial effort to plan and optimize EM follow-up. Previous work compared the advantages and disadvantages of different telescopes (Kasliwal & Nissanke 2014). Regardless of a facility's hardware, observatories at different locations will have systematically different opportunities to follow-up GW events due to properties of the GW detector network.

Even if astrophysical sources are distributed isotropically on the sky, GW detections with the two LIGO detectors will not be. Because of the detectors' locations (Hanford, WA and Livingston, LA) and duty cycles, detectable GW sources preferentially come from certain locations on the celestial sphere, and the preferred regions vary with the seasons. Follow-up of EM counterparts with emission timescales of less than a few weeks will be especially affected by these biases. In particular, the Earth's rotation limits ground-based EM follow-up facilities. Different sites can access different parts of the GW localization maps at different times, and therefore, have different expectations for the fraction of counterparts they can detect. These GW selection effects have other important ramifications for follow-up efforts, such as the amount of time before an average source will be accessible and the average air mass expected.

GW detectors do not operate continuously (Aasi et al. 2015; Abbott et al. 2016b). Previous studies considered their duty cycle when estimating detection rates (Singer et al. 2014; Abbott et al. 2016c; The LIGO Scientific Collaboration et al. 2016b). However, their operation is not uniformly distributed in time, and instead shows a strong preference for acquiring data during their local night. We describe this behavior quantitatively and assess its influence on GW detections as well as EM follow-up.

We first explain the sources of biases and their impact on GW detections in Section 2. We show the effects of the bias on ground-based EM follow-up facilities in Section 3, and conclude in Section 4.

Ground-based GW interferometers do not have isotropic sensitivity. Over time, the antenna pattern produces systematic preferences for the locations of detectable sources. These can be split into a dependence on the source's declination, which will not improve until additional detectors are added to the network, and a dependence on the source's right ascension, which can be mitigated with existing facilities.

2.1. Dependence on Declination

Because of projection effects, GW detectors are most sensitive to signals coming from above or below the plane defined by their arms (Thorne 1987). The two LIGO detectors, for example, are most sensitive, and equally sensitive, to sources directly above North America and above the Indian Ocean (Figure 1(a)). As the Earth rotates, the antenna pattern is swept across the celestial sphere, creating preferred bands in the mid-latitudes of both hemispheres (Figure 1(b)). This is determined primarily by the detectors' relative sensitivity and the geometry of the network. For the two LIGO detectors, we expect more detections at mid-declinations compared to low- or high-declinations.

Zoom In Zoom Out Reset image size

2.2. Dependence on Right Ascension

The antenna pattern also introduces a dependence on longitude, which translates to a dependence on the right ascension. If the detectors operated uniformly in time, then this right ascension dependence would average away as the Earth rotates. However, because the duty cycle is not uniform in time, this engenders a net bias in favor of certain right ascensions.

Ground-based GW detectors' data acquisition shows a clear diurnal cycle. In previous runs (S6—Abadie et al. 2012b; Aasi et al. 2013—and earlier) this was due primarily to anthropogenic noise; the detectors observed lower ambient noise at night because humans were less active. During O1, we observed a similar diurnal cycle. Although anthropogenic noise was mitigated by improved seismic isolation, commissioning activities still preferentially occurred during the day at the detector sites. Furthermore, the detectors achieved an overall coincident duty cycle of 50% (Abbott et al. 2016). This generated a non-trivial preference for acquiring data during the site's night (Figure 1d). We model this preference with a probability distribution for the time of day during which both interferometers are likely to be operating:

$$\begin{eqnarray}&&p(t| \mathrm{operating})=\displaystyle \frac{1}{\mathrm{day}}\left(1+A\sin \left(\displaystyle \frac{2\pi t}{\mathrm{day}}-\displaystyle \frac{\pi }{8}\right)\right).\end{eqnarray} \tag{ 1 }$$

The amplitude (A) reflects the extent of the day/night bias and is typically near 0.4. This means that the detectors are 2.3 times more likely to record data at their local midnight than at noon. We note that this bias may be decreased by reducing the amount of commissioning activity during the day or by increasing the overall duty cycle of the instruments (e.g., reducing the duration of downtime due to events uncorrelated with the diurnal cycle). As the duty cycle increases, the relative importance of the diurnal behavior will decrease (see the Appendix for more details). However, we always expect some small diurnal cycle to be present.

The diurnal cycle preserves the dependence on right ascension as the Earth rotates (Figure 1(c)), which persists over timescales of days to weeks. Figure 2 shows the dependence averaged over a month for several months throughout the year. However, over the course of a year, the Earth's orbit will average away this dependence. Nonetheless, as currently scheduled, GW detectors do not operate year-round, instead only recording data over a few consecutive months during observing runs. If observing runs are scheduled during the same season repeatedly, the right ascension dependence introduced by the diurnal cycle can persist for years.

Zoom In Zoom Out Reset image size

Furthermore, the galactic plane intersects the preferred directions from May through September (Figure 2). Depending on the desired target, this could be an advantage or a hindrance. Detected compact binaries are expected to be extragalactic (Abbott et al. 2016a, 2016e) and the galactic plane will serve as a foreground for any EM follow-up, significantly complicating the removal of transient contaminants. However, detectable core-collapse supernova are expected to be primarily galactic (Gossan et al. 2016), and therefore GW detectors may have their best chance of observing such events during the north's summer. Regardless of the source, a full understanding of the distribution of detectable signals across both declination and right ascension will be crucial when considering any isotropy or homogeneity measurements using GW observations alone.

Selection effects associated with ground-based GW detectors impact EM follow-up facilities. In particular, we focus on three possible effects: the localization probability that a telescope can survey (observable probability), the source's closest angle of approach to an observatory's zenith while the observatory can observe (mimimum zenith distance), and the time until a GW source becomes observable (delay time). We focus on EM follow-up timescales of up to a few days or weeks, and therefore neglect seasonal modulations of the sky over the duration of the EM observations for each individual GW event. This is an appropriate timescale for short gamma-ray burst afterglows (Berger 2014) and kilonovae (Barnes & Kasen 2013), two promising EM counterparts of compact binary coalescences involving at least one neutron star. Radio transients may persist over longer timescales and therefore these effects may be less relevant. For these EM follow-up timescales, we consider low-latency GW alerts. In O1, it took a couple of days to issue alerts, but we expect this to be reduced to a few minutes in O2 and beyond (Abbott et al. 2016f). Throughout this paper, we assume the diurnal cycle modeled in Equation (1) as well as $18^\circ$ of astronomical twilight (Patat et al. 2006). We also assume that the observatories can observe within $90^\circ$ of their zenith.

In the limit of a large number of detections, the combined posterior distributions trace out the network antenna pattern. In what follows, we use the antenna pattern to approximate limiting distribution of many events (see the Appendix for more details). To quantify departures from this limit when only a few events are available, we simulate collections of events based on binary black hole localizations with two detectors (Essick et al. 2015; Essick 2015). We find good agreement.

3.1. Observable Probability

GW localizations are driven by triangulation and, for networks of two detectors, the localization is characterized by large rings that can span hundreds of square degrees regardless of source morphology. Furthermore, these error regions typically have support at antipodal points on the sky, making it difficult for a single EM observatory to access the entire skymap. This is compounded by solar occlusion, which renders certain parts of the sky inaccessible. We ignore lunar occlusion because, unlike solar occlusion, it is not thought to systematically correlate with when the detectors operate. Because we focus on timescales of a few days, we assume that the Earth will revolve at least once. Therefore, the observable region depends only on the observatory's latitude.

Figure 3 shows the observable probability as a function of latitude for a year-long average and for the solstices. We define

$$\begin{eqnarray}{p}_{\mathrm{obs}}({\mathrm{lat}}_{\mathrm{site}}) & = & \displaystyle \int {dt}\,\left(\phantom{{{l}_{{i}_{i}}}_{{i}_{i}}\,}p(t| \mathrm{operating})\right.\\ & & \left.\times \displaystyle \int d{\rm{\Omega }}\,{p}_{\mathrm{GW}}({\rm{\Omega }},t){{\rm{\Theta }}}_{\mathrm{obs}}({\rm{\Omega }},t,{\mathrm{lat}}_{\mathrm{site}})\right),\end{eqnarray} \tag{ 2 }$$

where ${p}_{\mathrm{GW}}({\rm{\Omega }},t)$ is the probability density of the GW source coming from Ω (at time t), and ${{\rm{\Theta }}}_{\mathrm{obs}}$ is the region accessible from a particular latitude at t, respectively. We see sharp declines near latitudes of $\pm 50^\circ$ corresponding to the Arctic/Antarctic circles and astronomical twilight, but otherwise p_obs follows a smooth distribution favoring equatorial observatories. This is because equatorial observatories can systematically access sources in both hemispheres, whereas other observatories may be confined to only the probability within their own. While equatorial facilities are favored overall, mid-latitudes have larger p_obs than they would for isotropically distributed detections. This will persist regardless of the day/night cycle and is driven solely by the GW detectors' latitudes. We note that the optimal observing months of the northern and southern hemispheres are out of phase (Figure 3). For example, if LIGO operates from September to February, the northern hemisphere will have a better chance of observing counterparts than the southern hemisphere.

Zoom In Zoom Out Reset image size

While p_obs does not depend on the longitude, the diurnal cycle can still introduce a systematic bias. This is because the Sun will be systematically out of phase with the northern maximum of the antenna pattern (Figure 2). We find that this typically produces an increase in p_obs of a few percent for mid-latitudes in the north, even though their analogs in the south have larger p_obs because of the shorter solar exposure over an entire year.

We note that p_obs reflects the amount of localization probability that is observable in the limit of many detections. Outside of this limit, we consider an analogous quantity (${\hat{p}}_{\mathrm{obs}}$) defined for a finite number of detections (N_d) with the correspondence ${p}_{\mathrm{obs}}={\mathrm{lim}}_{{N}_{d}\to \infty }\,{\hat{p}}_{\mathrm{obs}}$ (see the Appendix for more details). Statistical fluctuations in ${\hat{p}}_{\mathrm{obs}}$, caused by variations in which events occur, can be large, particularly with $\lesssim 10$ events. We expect the uncertainty in this estimate to scale inversely with $\sqrt{{N}_{d}}$, and Figure 3 reports ${\mathrm{lim}}_{{N}_{d}\to \infty }\,{({N}_{d}/10)}^{1/2}{\sigma }_{{\hat{p}}_{\mathrm{obs}}}$ as error bars along with the mean. Typical values are between 6% and 10%. We also note that the intrinsic distribution for single events (N_d = 1) may not be Gaussian but, in the limit of many detections, the distribution of the mean will be.

3.2. Minimum Zenith Distance

While equatorial observatories may be able to access the largest integrated probability, and thereby have the largest probability of being able to image a source, this does not necessarily imply that they will have the best conditions for observing. An important consideration is the closest approach of each field to an observatory's zenith (minimum zenith distance). Sources at large angles from an observatory's zenith can be difficult to observe because of high air mass and mechanical limitations. In Figure 4, we present the mean minimum zenith distance, as a function of latitude, weighted by the probability that the source actually comes from each location:

$$\begin{eqnarray}{D}_{\mathrm{zen}}({\mathrm{lat}}_{\mathrm{site}}) & = & \displaystyle \frac{1}{{p}_{\mathrm{obs}}}\displaystyle \int {dt}\,(p(t| \mathrm{operating})\\ & & \times \displaystyle \int d{\rm{\Omega }}\,{p}_{\mathrm{GW}}({\rm{\Omega }},t){{\rm{\Theta }}}_{\mathrm{obs}}({\rm{\Omega }},t,{\mathrm{lat}}_{\mathrm{site}})\\ & & \times {D}_{\mathrm{zen}}({\rm{\Omega }},t,{\mathrm{lat}}_{\mathrm{site}}))\end{eqnarray} \tag{ 3 }$$

${{\rm{\Theta }}}_{\mathrm{obs}}$ accounts for solar occlusion and ${D}_{\mathrm{zen}}({\rm{\Omega }},t,{\mathrm{lat}}_{\mathrm{site}})$ incorporates when the EM facility will actually be able to observe. We find that observatories at extreme latitudes ($\pm 90^\circ$) will almost always have large D_zen, with a gradual transition to lower values at lower latitudes. There is a $\sim 10 \%$ difference in D_zen between observatories at mid-latitudes, with northern sites preferred. This is because the diurnal cycle makes the Sun preferentially overlap the southern antenna pattern and forces southern facilities to observe closer to sunrise and sunset than their northern counterparts. This behavior is particularly evident at the northern winter solstice.

Zoom In Zoom Out Reset image size

Figure 4 also shows the fluctuations in the analgous stastic defined for a finite number of detections (${\hat{D}}_{\mathrm{zen}}$). We typically find ${\mathrm{lim}}_{{N}_{d}\to \infty }{({N}_{d}/10)}^{1/2}{\sigma }_{{\hat{D}}_{\mathrm{zen}}}\sim 4$–10$^\circ$ (see the Appendix for more details). Furthemore, D_zen corresponds to the mean of many events. For a single event, the mode of ${\hat{D}}_{\mathrm{zen}}({N}_{d}=1)$'s distribution falls near $20^\circ$ at mid-latitudes in the north and $60^\circ$ in the south, whereas the mode is near $50^\circ$ at both poles.

3.3. Delay Time

For counterparts with timescales of days, EM observatories' longitudes can play an important role. This is because observatories will have to wait to begin observing until the source rises at their location. This can be exacerbated by the position of the Sun, which will systematically correlate with the diurnal cycle.

We expect that the amount of time an observatory must wait before commencing observations (delay time) will depend on both the observatory's longitude and latitude and define

$$\begin{eqnarray}{D}_{\mathrm{del}}({{\rm{\Omega }}}_{\mathrm{site}}) & = & \displaystyle \frac{1}{{p}_{\mathrm{obs}}}\displaystyle \int {dt}\,(p(t| \mathrm{operating})\\ & & \times \displaystyle \int d{\rm{\Omega }}\,{p}_{\mathrm{GW}}({\rm{\Omega }},t){{\rm{\Theta }}}_{\mathrm{obs}}({\rm{\Omega }},t,{\mathrm{lat}}_{\mathrm{site}})\\ & & \times {D}_{\mathrm{del}}({\rm{\Omega }},t,{{\rm{\Omega }}}_{\mathrm{site}}))\end{eqnarray} \tag{ 4 }$$

where ${D}_{\mathrm{del}}({{\rm{\Omega }}}_{\mathrm{site}},t,{\rm{\Omega }})$ accounts for sunrise and sunset along with the source's relative position to the observatory. We restrict ourselves to only the parts of the skymap that are actually accessible from each site. Figure 5 shows the dependence on the observatory's longitude and latitude. Although the shape changes, we see a reasonably uniform distribution of D_del throughout the globe (there are slightly longer delay times in the north compared to their analogs in the south away from the blob over Africa) with the notable exception of the Southern Atlantic, Africa, and the Indian Ocean, particularly in the year-long average (Figure 5(e)). Because the majority of detections will occur during North America's night, these locations are likely to already be in daylight and will therefore have to wait for sunset. Furthermore, observatories in the south will have to wait until their zenith is very close to the source's right ascension before they observe, which explains why the blob is wider in the south. From simulations outside the limit of many detections, we also note that the distribution of the analogous statistic (${\hat{D}}_{\mathrm{del}}({N}_{d}=1)$) is very non-Gaussian and very skewed to the right. However, the mean is still Gaussian as ${N}_{d}\to \infty$ with an associated variance ${\mathrm{lim}}_{{N}_{d}\to \infty }{({N}_{d}/10)}^{1/2}{\sigma }_{{\hat{D}}_{\mathrm{del}}}\,\sim$ 80–240 minutes, with most values near 100 minutes, depending on the observatory's location (see the Appendix for more details).

Zoom In Zoom Out Reset image size

We note that if the diurnal cycle were to be strongly reduced or eliminated, the blob above Africa would become less prominent and an analogous blob would appear in the Northern Pacific. The latter corresponds to observatories that are in daylight when a source is detected above the southern antenna pattern. The overall effect of removing the diurnal cycle, however, is toward a more uniform distribution of delay times across the entire globe.

In this paper, we enumerate three effects that imprint selection effects on the distribution of GW detections, which will be relevant for all types of GW sources: (1) the detector's locations introduce a preference for mid-declinations, (2) a diurnal cycle modulates when the detectors operate and produces a sky sensitivity with a dependence on the right ascension over short timescales, and (3) if detections are made only during relatively short observing runs, the right ascension bias can be imprinted over longer timescales. While the effects are important when modeling the expected distribution of GW detections on the sky, they can also have significant implications for EM follow-up.

In general, ground-based EM observatories located at latitudes comparable to those of the GW detectors are preferred. While equatorial sites can access the largest fraction of the sky, mid-latitude sites in the north will be able to observe more events closer to their zenith without sacrificing much coverage. Furthermore, there is a preference for EM observatories located near the same longitude as the GW detectors, or slightly west thereof, assuming the EM facilities can begin observing quickly and that counterparts decay with timescales of hours. This leads us to the conclusion that EM facilities located near the GW detectors are favored by how the GW detectors' actually operate.

In addition to informing follow-up efforts for GW candidates, our results may help guide EM activities in advance of GW observations. For example, because we know which parts of the celestial sphere are most likely to host detectable GW events, we can build up relevant templates for image subtraction before GW observations even begin, focusing particularly on the highest probability regions. Galaxy catalog constructions could be focused similarly. In addition, surveys may focus their observations near the peak of the instantaneous antenna pattern, thereby anticipating where detectable GW sources are most likely to occur and facilitating target of opportunity follow-ups. This will also increase the probability of serendipitous detection of prompt EM counterparts. We note that EM observatories located near GW detectors will naturally survey a maximum of the antenna pattern because their zenith lies near that maximum automatically.

The behavior of GW observatories can become more complicated as the global network of detectors expands. A larger network increases the uniformity of sensitivity to GW signals across the sky, and therefore reduces the selection effects in both declination and right ascension. Furthermore, detectors located around the globe will likely experience diurnal cycles that are out of phase, further reducing any preference for certain right ascensions. However, because all planned detectors lie within a relatively confined band of latitudes (Aso et al. 2013; Fairhurst 2014; Acernese et al. 2015), the bias on declination may persist at some level. In addition, the LIGO detectors will likely have the best sensitivity in the network for at least the next few years (Abbott et al. 2016c). This means that they will generally provide the dominant contribution to detections, and the observed distribution of sources will follow their antenna pattern. Thus, our analysis serves as a reasonable prediction for both the distributions of GW detections as well as their impact on EM follow-up efforts for the next few years.

The authors want to thank Lisa Barsotti, Matt Evans, and Rich Middleman for discussions about locking behavior, commissioning, and anthropogenic noise sources, James Annis for discussions about EM observations, and Marica Branchesi for very helpful comments. They also acknowledge the support of the National Science Foundation and the LIGO Laboratory. LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-0757058. The authors would like to acknowledge the LIGO Data Grid clusters. Specifically, we thank the Albert Einstein Institute in Hannover, supported by the Max-PlanckGesellschaft, for use of the Atlas high-performance computing cluster.

D.E.H. and H.Y.C. were supported by NSF CAREER grant PHY-1151836. They also acknowledge support from the Kavli Institute for Cosmological Physics at the University of Chicago through NSF grant PHY-1125897 as well as an endowment from the Kavli Foundation. This is LIGO document P1600219.

Throughout the main paper, we use the antenna pattern to approximate limiting distributions of many events. We also study the distributions by drawing simulated localization maps. We repeatedly draw N_d maps, distributed through time according to $p(t| \mathrm{operating})$, to obtain the distributions of our statistics for N_d detections. In this approach, the observable probability, the minimum zenith distance, and the delay time are calculated as

$$\begin{eqnarray}&&{\hat{p}}_{\mathrm{obs}}({\mathrm{lat}}_{\mathrm{site}})=\displaystyle \frac{1}{{N}_{d}}\displaystyle \sum _{i=1}^{{N}_{d}}\int d{\rm{\Omega }}\,{p}_{\mathrm{sky},i}{{\rm{\Theta }}}_{\mathrm{obs}}({\rm{\Omega }},{t}_{i},{\mathrm{lat}}_{\mathrm{site}}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\hat{D}}_{\mathrm{zen}}({\mathrm{lat}}_{\mathrm{site}})=\displaystyle \frac{1}{{N}_{d}\,{\hat{p}}_{\mathrm{obs}}}\displaystyle \sum _{i=1}^{{N}_{d}}\int d{\rm{\Omega }}\,{p}_{\mathrm{sky},i}{{\rm{\Theta }}}_{\mathrm{obs},i}{D}_{\mathrm{zen}}({\rm{\Omega }},{t}_{i},{\mathrm{lat}}_{\mathrm{site}}),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\hat{D}}_{\mathrm{del}}({{\rm{\Omega }}}_{\mathrm{site}})=\displaystyle \frac{1}{{N}_{d}\,{\hat{p}}_{\mathrm{obs}}}\displaystyle \sum _{i=1}^{{N}_{d}}\int d{\rm{\Omega }}\,{p}_{\mathrm{sky},i}{{\rm{\Theta }}}_{\mathrm{obs},i}{D}_{\mathrm{del}}({\rm{\Omega }},{t}_{i},{{\rm{\Omega }}}_{\mathrm{site}}),\end{eqnarray} \tag{ 7 }$$

where ${p}_{\mathrm{sky},i}$ is the localization map probability for the ${i}_{\mathrm{th}}$ detection and ${\hat{p}}_{\mathrm{obs}}$ in Equations (6) and (7) is computed using the same set of skymaps as the explicit sum. We note that using the antenna pattern yields the mean of these statistics when simulating many detections. Namely, if we calculate χ using the antenna pattern and $\hat{\chi }$ using sets of simulated maps, we expect

$$\begin{eqnarray}\chi & = & \mathop{\mathrm{lim}}\limits_{{N}_{d}\to \infty }\hat{\chi }\\ & = & \mathop{\mathrm{lim}}\limits_{{N}_{d}\to \infty }\displaystyle \frac{1}{{N}_{d}}\displaystyle \sum _{i=1}^{{N}_{d}}\displaystyle \int d{\rm{\Omega }}\,{p}_{\mathrm{sky},i}\chi ({t}_{i}).\end{eqnarray} \tag{ 8 }$$

We note that $\chi ({t}_{i})$ only depends on the time the event occurs through the position of the Sun; the dependence will be the same for all detections that occur at the same time. If we break the sum into small segments of time, we can write

$$\begin{eqnarray}&&\chi =\mathop{\mathrm{lim}}\limits_{{N}_{d}\to \infty }\displaystyle \frac{1}{{N}_{d}}\displaystyle \sum _{j=1}^{{N}_{t}}\,\displaystyle \sum _{\mathop{{t}_{i}\in [{t}_{j},{t}_{j}+{\rm{\Delta }}t)}\limits^{i=1}}^{{N}_{d,j}}\int d{\rm{\Omega }}\,{p}_{\mathrm{sky},i}\chi ({t}_{i}),\end{eqnarray} \tag{ 9 }$$

where ${\sum }_{j}{N}_{d,j}={N}_{d}$. Now, when ${N}_{d}\to \infty$, we can make the segments as small as we like while maintaining a large number of detections in each bin. We then obtain

$$\begin{eqnarray}\chi & = & \mathop{\mathrm{lim}}\limits_{\mathop{{N}_{t}{\rm{\Delta }}t=T}\limits^{{N}_{t}\to \infty }}\displaystyle \sum _{j=1}^{{N}_{t}}\,\mathop{\mathrm{lim}}\limits_{{N}_{d,j}\to \infty }\displaystyle \frac{1}{{N}_{d}}\displaystyle \sum _{\mathop{{t}_{i}\in [{t}_{j},{t}_{j}+{\rm{\Delta }}t)}\limits^{i=1}}^{{N}_{d,j}}\displaystyle \int d{\rm{\Omega }}\,{p}_{\mathrm{sky},{\rm{i}}}\chi ({t}_{i})\\ & = & \mathop{\mathrm{lim}}\limits_{\mathop{{N}_{t}{\rm{\Delta }}t=T}\limits^{{N}_{t}\to \infty }}\displaystyle \sum _{j=1}^{{N}_{t}}\,\mathop{\mathrm{lim}}\limits_{{N}_{d,j}\to \infty }\left(\displaystyle \frac{{N}_{d,j}}{{N}_{d}}\right)\displaystyle \int d{\rm{\Omega }}\left(\displaystyle \frac{1}{{N}_{d,j}}\displaystyle \sum _{\mathop{{t}_{i}\in [{t}_{j},{t}_{j}+{\rm{\Delta }}t)}\limits^{i=1}}^{{N}_{d,j}}{p}_{\mathrm{sky},{\rm{i}}}\right)\chi ({t}_{i})\\ & = & \mathop{\mathrm{lim}}\limits_{\mathop{{N}_{t}{\rm{\Delta }}t=T}\limits^{{N}_{t}\to \infty }}\displaystyle \sum _{j=1}^{{N}_{t}}(p(t| \mathrm{operating}){\rm{\Delta }}t)\displaystyle \int d{\rm{\Omega }}({p}_{\mathrm{GW}}({\rm{\Omega }},{t}_{j}))\chi ({t}_{j})\\ & = & \displaystyle \int {dt}\,p(t| \mathrm{operating})\displaystyle \int d{\rm{\Omega }}\,{p}_{\mathrm{GW}}({\rm{\Omega }},t)\chi (t)\end{eqnarray} \tag{ 10 }$$

where T is the length of the observing season. We assumed ${t}_{i}\to {t}_{j}$ for all events within each bin and used the fact that many localization posteriors stacked on top of one another will average to the network antenna pattern in equatorial coordinates (${p}_{\mathrm{GW}}({t}_{j})$), which depends on the bin's time. We have also used the fact that the fraction of events occuring in each bin is equal to the probability that the detectors are operating throughout that bin ($p(t| \mathrm{operating}){\rm{\Delta }}t$). By approximating these integrals, we obtain the limits of the means much more efficiently than through direct simulation. This proceedure, or an equivalent, is used in Figures 3–5.

We also note that we normalize by the total observable probability in Equations (3), (4), (6), and (7). This is because we restrict ourselves to only the fraction of the probability that is actually observable. The particular form of our normalization (dividing by ${\hat{p}}_{\mathrm{obs}}$ for a set of events rather than each event separately) guarantees that we sample the antenna pattern in the limit.

Section 3 presents the limit of our statistics when many detections are present. In this limit, the distributions of these statistics will be Gaussian, but it is also informative to examine the distributions when there are only a few events. This is particularly important when considering how an observatory could be impacted for any individual event, rather than for a collection of events. We also consider these distributions in the limit of extremely well localized sources: ${p}_{\mathrm{sky},i}\to \delta ({\rm{\Omega }}-{{\rm{\Omega }}}_{i})$. This point-source limit describes the distributions obtained when only observing the location of the true source.

These distributions were computed using bootstrapped simulations of BBH detections based off of Essick et al. (2015). The library developed to perform these simulations is publicly available (https://github.com/reedessick/selectionEffects) and readers are encouraged to use it to determine distributions for their favorite observatory.

B.1.  ${\hat{p}}_{\mathrm{obs}}$

Of all the statistics we consider, ${\hat{p}}_{\mathrm{obs}}$ is the most Gaussian for small N_d for year-long averages. However, because it is bounded from above and below, the distribution does deviate at times. Figure 6 shows the distributions for a few latitudes and a few values of N_d.

Zoom In Zoom Out Reset image size

We note that, in the limit ${N}_{d}\to \infty$, p_obs is the fraction of true counterparts that an observatory can observe. When considering extremely well localized sources, the observatory will either be able to observe the true source or not, and each trial will have nearly the same probability of success, modulo variations caused by the Sun's declination. Therefore, the fraction of true sources each observatory can observe will be nearly binomially distributed. Figure 7 plots the point-source limit of these distributions (${\hat{p}}_{\mathrm{obs}}^{(\mathrm{pt}\,\mathrm{src})}$), again with several values of N_d.

Zoom In Zoom Out Reset image size

B.2.  ${\hat{D}}_{\mathrm{zen}}$

${\hat{D}}_{\mathrm{zen}}({N}_{d}=1)$ is also fairly Gaussian for some latitudes, but there can often be non-trivial deviations therefrom. In particular, mid-latitudes may show interesting skew right distributional shapes. Figure 8 shows these distributions for a few latitudes as a function of N_d. We note that these distributions are much narrower than the point-source limit in which all events are well localized (${\hat{D}}_{\mathrm{zen}}^{(\mathrm{pt}\,\mathrm{src})}$), shown in Figure 9. There is more "shot noise" for ${\hat{D}}_{\mathrm{zen}}^{(\mathrm{pt}\,\mathrm{src})}$, which broadens the distributions, but the general distributional shapes are similar for both ${\hat{D}}_{\mathrm{zen}}$ and ${\hat{D}}_{\mathrm{zen}}^{(\mathrm{pt}\,\mathrm{src})}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

B.3.  ${\hat{D}}_{\mathrm{del}}$

We note in Section 3.3 that the distribution of ${\hat{D}}_{\mathrm{del}}$ is very skewed right when N_d is small. Figure 10 demonstrates this. Typically, there is an extremely large lobe near zero, corresponding to events that are immediately observable, and a long tail comprised of events that require waiting. When we consider the point-source limit of well localized events, this behavior is enhanced. Figure 11 demonstrates that. When we integrate over typical skymaps, instead of point sources, we find that the peak is smeared out to longer delay times. This is because different parts of the skymap may become observable at different times, and that fuzz tends to smooth the distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While Equation (1) is fairly self evident from the inspection of Figure 1, we can derive its form from more basic assumptions about how human activity may cause downtime. We posit two states of a network of detectors: up (u) and down (d), with science-quality data available only in the up state. Furthermore, we posit two possible causes for detectors being in thedown state: random causes (r) that are uncorrelated with time and cyclic causes (c) which are correlated with time, usually through a diurnal cycle. We note that the random and cyclic models are not necessarily mutually exclusive and the detector could be down for mutliple reasons at the same time.

Therefore, we have

$$\begin{eqnarray*}p(d| t) & = & p(d| c,t)p(c| t)\\ & & +p(d| r,t)p(r| t)-p(d| c\cap r,t)p(c\cap r| t),\end{eqnarray*}$$

which implies

$$\begin{eqnarray*}p(d) & = & \displaystyle \int {dtp}(t)p(d| t)\\ & = & \displaystyle \int {dt}\,p(t)[p(d| c,t)p(c| t)+p(d| r,t)p(r| t)\\ & & -p(d| c\cap r,t)p(c\cap r| t)]\end{eqnarray*}$$

Furthermore, because the up and down states are mutually exclusive and span the space of possible detector states at any single time, we have

$$\begin{eqnarray*}&&p(u| t)=1-p(d| t)\Rightarrow p(u)=1-p(d)\end{eqnarray*}$$

and Bayes theorem yields

$$\begin{eqnarray*}&&p(t| u)=\displaystyle \frac{p(u| t)p(t)}{p(u)}=\displaystyle \frac{(1-p(d| t))p(t)}{1-p(d)}.\end{eqnarray*}$$

We measure p(u) and $p(t| u)$ in a straightforward manner from the data used in Figure 1. Typically, we assume some periodicity in $p(t| u)$ to generate a histogram with enough samples to be statistically meaningful, but we expect to be able to identify the periodic elements of $p(t| u)$ through Fourier analysis as well.

However, we are also interested in slightly different probabilities. Namely, we would like to know

$$\begin{eqnarray*}&&p(c| d,t)=\displaystyle \frac{p(d| c,t)p(c| t)}{p(d| t)},\end{eqnarray*}$$

$$\begin{eqnarray*}&&p(r| d,t)=\displaystyle \frac{p(d| r,t)p(r| t)}{p(d| t)}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&p(c\cap r| d,t)=\displaystyle \frac{p(d| c\cap r,t)p(c\cap r| t)}{p(d| t)}\end{eqnarray*}$$

which express the probabilities that the detector is down due to a particular cause at a specific time. If we only care about long-term averages (over timescales much longer than the cyclic model's periodicty), then we can marginalize away the time dependence:

$$\begin{eqnarray*}&&p(c| d)=\int {dtp}(t)\displaystyle \frac{p(d| c,t)p(c| t)}{p(d| t)}\end{eqnarray*}$$

$$\begin{eqnarray*}&&p(r| d)=\int {dtp}(t)\displaystyle \frac{p(d| r,t)p(r| t)}{p(d| t)}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&p(c\cap r| d)=\int {dtp}(t)\displaystyle \frac{p(d| c\cap r,t)p(c\cap r| t)}{p(d| t)}.\end{eqnarray*}$$

These equations should hold regardless of the specific form of the cyclic and random models.

By definition, we assume that the random causes are uncorrelated with time so that $p(d| r,t)=p(d| r)$. Furthermore, we assume some periodic function for the cyclic model so that $p(d| c,t)=p(d| c,t+\tau )$ for some τ. Specifically, we expand the periodic function in terms of the oscillating (AC) and constant (DC) components

$$\begin{eqnarray*}&&p(d| c,t)={p}_{\mathrm{DC}}(d| c)+{p}_{\mathrm{AC}}(d| c,t)\end{eqnarray*}$$

such that

$$\begin{eqnarray*}&&{p}_{\mathrm{DC}}={\int }_{0}^{\tau }{dt}\,p(d| c,t).\end{eqnarray*}$$

Clearly, we require ${p}_{\mathrm{DC}}(d| c)\geqslant -{p}_{\mathrm{AC}}(d| c,t)\,\forall t$. We also assume the priors for the causes do not depend on time ($p(r| t)=p(r)$ and $p(c| t)=p(c)$) and that the priors are equal for the two causes ($p(r)=p(c)=p$).

Note that by measuring p(d), we only extract the combination of ${p}_{\mathrm{DC}}(d| c)+p(d| r)$ and cannot separate these terms further. However, by measuring $p(t| u)$ as well, we are able to determine ${p}_{\mathrm{AC}}(d| c,t)$ from which we can determine ${p}_{\mathrm{DC}}(d| c)$ by requiring that ${\min }_{t}\{p(d| c,t)\}=0$. Any other DC component to the cyclic model is indistinguishable from the random model and therefore we lump it together with the random model.

If we allow the cause models to overlap but require them to be independent, we obtain $p(d| c\cap r,t)=p(d| c,t)p(d| r,t)$. This implies

$$\begin{eqnarray*}&&p(r)+p(c)-p(r)p(c)=1\Rightarrow p(r)=p(c)=p=1\end{eqnarray*}$$

and the two observables are

$$\begin{eqnarray*}&&p(d)=p(d| r)p+(1-p(d| r)p){p}_{\mathrm{DC}}(d| c)p\end{eqnarray*}$$

$$\begin{eqnarray*}&&p(t| u)=p(t)\displaystyle \frac{1-p(d| r)p-(1-p(d| r)p){p}_{\mathrm{DC}}(d| c)p-(1-p(d| r)p){p}_{\mathrm{AC}}(d| c,t)p}{1-p(d| r)p-(1-p(d| r)p){p}_{\mathrm{DC}}(d| c)p}=p(t)\left(1-\displaystyle \frac{(1-p(d| r)p)p}{1-p(d)}{p}_{\mathrm{AC}}(d| c,t)\right).\end{eqnarray*}$$

This yields

$$\begin{eqnarray*}&&p(c| d,t)\\ &&\quad =\displaystyle \frac{{p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c)}{({p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c))(1-p(d| r)p)+p(d| r)}\end{eqnarray*}$$

$$\begin{eqnarray*}&&p(r| d,t)\\ &&\quad =\displaystyle \frac{p(d| r)}{({p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c))(1-p(d| r)p)+p(d| r)}\end{eqnarray*}$$

$$\begin{eqnarray*}&&p(c\cap r| d,t)\\ &&\quad =\displaystyle \frac{({p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c))p(d| r)p}{({p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c))(1-p(d| r)p)+p(d| r)}\end{eqnarray*}$$

and their marginalized counterparts

$$\begin{eqnarray*}&&p(c| d)=\displaystyle \int {dt}\,p(t| d)\\ &&\quad \times \displaystyle \frac{{p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c)}{({p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c))(1-p(d| r)p)+p(d| r)}\\ &&\quad =\displaystyle \frac{{p}_{\mathrm{DC}}(d| c)p}{p(d)}\end{eqnarray*}$$

$$\begin{eqnarray*}&&p(r| d)=\displaystyle \int {dt}\,p(t| d)\\ &&\quad \times \displaystyle \frac{p(d| r)}{({p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c))(1-p(d| r)p)+p(d| r)}\\ &&\quad =\displaystyle \frac{p(d| r)p}{p(d)}\end{eqnarray*}$$

$$\begin{eqnarray*}&&p(c\cap r| d)=\displaystyle \int {dt}\,p(t| d)\\ &&\quad \times \displaystyle \frac{({p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c))p(d| r)p}{({p}_{\mathrm{AC}}(d| c,t)+{p}_{\mathrm{DC}}(d| c))(1-p(d| r)p)+p(d| r)}\\ &&\quad =\displaystyle \frac{{p}_{\mathrm{DC}}(d| c)p\cdot p(d| r)p}{p(d)}.\end{eqnarray*}$$

A reasonable ansatz is $p(d| c,t)=B\left(1+\sin \left(\tfrac{2\pi t}{\tau }-\phi \right)\right)$, in which case we obtain

$$\begin{eqnarray*}&&p(d)=p(d| r)p+{Bp}(1-p(d| r)p)\end{eqnarray*}$$

$$\begin{eqnarray*}p(t| u) & = & p(t)\left(1+\displaystyle \frac{(1-p(d| r)p){pB}}{1-p(d)}\sin \psi \right)\\ & = & p(t)(1+A\sin \psi ),\end{eqnarray*}$$

which is exactly the form of Equation (1). We also see why a low duty cycle (small $p(u)=1-p(d)$) can amplify the amplitude of they day/night bias.

If we are interested in just the probability associated with the cyclic model regardless of the random model, we obtain

$$\begin{eqnarray*}&&p(c| d)=\displaystyle \frac{1-p(d)}{p(d)}\left(\displaystyle \frac{A}{1-p(d)+(1-p(d))A}\right)\approx 0.49,\end{eqnarray*}$$

which implies that nearly half the time cyclic causes were at least partly responsible for bringing down the detector network during O1. If we restrict ourselves to times when the cyclic causes were the sole cause of the downtime, we obtain

$$\begin{eqnarray*}&&p(c| d)-p(c\cap r| d)=\displaystyle \frac{(1-p(d))A}{p(d)}\approx 0.29,\end{eqnarray*}$$

which suggests that we could reduce the downtime by 30% if we completely removed cyclic causes of downtime.