Optimal Optical Search Strategy for Finding Transient in Large Sky Error Region Under Realistic Constraints

TL;DR

This paper develops optimized search strategies for rapidly identifying optical counterparts of gravitational wave sources within large sky error regions, improving efficiency by combining tiling, galaxy targeting, and scheduling constraints.

Contribution

It introduces two novel algorithms that incorporate airmass and slew optimization into telescope scheduling, enhancing detection probability and coverage over previous methods.

Findings

Slew-optimization reduces total slew angle by hundreds of degrees.

Airmass-weighted algorithm increases probability acquisition by up to 20%.

The methods improve coverage and efficiency in simulated sky maps.

Abstract

In order to identify the rapidly-fading, optical transient counterparts of gravitational wave (GW) sources, an efficient follow-up strategy is required. Since most ground-based optical observatories aimed at following-up GW sources have a telescope with a small field-of-view (FOV) as compared to the GW sky error region, we focus on a search strategy that involves dividing the GW patch into tiles of the same area as the telescope FOV to strategically image the entire patch. We present an improvement over the optimal telescope-scheduling algorithm outlined in Rana et al. (2016), by combining the tiling and galaxy-targeted search strategies, and factoring the effects of the source airmass and telescope slew, along with setting constraints, into the scheduling algorithm in order to increase the chances of identifying the GW counterpart. We propose two separate algorithms: the airmass-weighted algorithm, a specific solution to the Hungarian algorithm that maximizes probability acquired, while minimizing the image airmass, and the slew-optimization algorithm that minimizes the overall slew angle covered between images for the given probability acquired by the optimal telescope-scheduling algorithm in Rana et al. (2016). Using the observatory site of the GROWTH-India telescope as an example, we generate 100s of skymaps to test the performance of our algorithms. Our results indicate that slew-optimization can reduce the cumulative slew angle in the observing schedule by 100s of degrees, saving several of minutes of observing time without the loss of tiles and probability. Further, we demonstrate that as compared to the greedy algorithm, the airmass-weighted algorithm can acquire up to 20 % more probability and 30 sq. deg. more in areal coverage for skymaps of all sizes and configurations.

Tables

Table 1. Table 1 Properties and parameters of the optical observatories used in simulations to test the slew optimization algorithm.

Telescope Name	Location	FOV size	latitude	longitude	Altitude
		(deg2)			(m)
GROWTH-India	Ladakh, IN	0.49	32.78	78.96	4500
Himalayan Chandra Telescope	Ladakh IN	0.25	32.78	78.96	4500
Swope	Atacama Region, CL	3.68	-29.02	-70.67	2380
PTF	Mt. Palomar, CA	7.80	33.35	-116.86	1712

Table 2. Table 2 Telescope performance with varying slew rate. The table above displays the difference between the performance of the optimal and slew optimization algorithms with slew rates (v slew ) of 2, 5, and 10 deg./s in terms of cumulative probability, slew angle in degrees, and slew time (t slew ) in minutes. The slew times here have already been factored into the tile observation times. The sample observatory is at the location of Swope, with a one square degree field of view. With increasing slew rate, the overall slew angle increases, as does the probability, while the slew time decreases. As expected, slew optimization is more critical for telescopes with slower slew rates.

$v_{slew}$ (deg./s)	algorithm	prob.	${\sigma }_{slew}(deg.)$	$t_{slew}(min.)$
2	optimal	0.683	2459	23.42
	slewopt	0.714	2095	17.47
5	optimal	0.728	3099	10.33
	slewopt	0.730	2501	8.33
10	optimal	0.737	3283	5.47
	slewopt	0.739	2616	4.31

Table 3. Table 3 Telescope performance with varying FOV. The table above displays the difference between the performance of the optimal and slew optimization algorithms with telescope FOVs of 1.0, 2.0 and 3.0 square degrees in terms of cumulative probability, slew angle in degrees, and slew time in minutes. The slew times here have already been factored into the tile observation times. The sample observatory is at the location of Swope, with a slew rate of 2 deg/s. With increasing FOV, the overall slew angle increases, as does the probability, and slew time. This is because larger FOV instruments can cover a given patch faster than a small FOV instrument, and therefore can afford to spend more of its total observing time in slewing. We find that slew optimization is more important for smaller FOV telescopes.

FOV (deg2)	algorithm	prob.	${\sigma }_{slew}$ (deg.)	$t_{slew}(min.)$
1.0	optimal	0.706	2810	23.4
	slewopt	0.714	2095	17.5
2.0	optimal	0.888	4891	32.5
	slewopt	0.890	4399	36.7
3.0	optimal	0.945	3899	40.8
	slewopt	0.949	3373	28.1

Equations

$${\sigma }_{slew}=co{s}^{-1}(\sin al{t}_1\ast \sin al{t}_2+\cos al{t}_1\ast \cos al{t}_2\ast \cos \Delta az)$$

$$\begin{array}{r}{t}_{gained}=(20/v_{slew})\ast N_{tiles}-t_{slew,exact}+\\ ((T>20{)}_{optimal}-(T>20{)}_{slew})\end{array}$$

$$am(t)=\{\begin{array}{ll}\frac{1}{cos(90-alt(t))},& \text{if }alt(t)\ge alt(horizon)\\ \infty ,& \text{if }alt(t)<alt(horizon)\end{array}$$

$$p_{amw}=Norm\left[\frac{p_{gw}}{am}\right]$$

$$t_{exp}=\frac{d_{gal}^2}{d_0^2}{t}_{exp0}$$

$$t_{exp}=am(t)t_{exp0}$$

Full text

1. Introduction

Within the next few years of gravitational-wave (GW) astrophysics, current estimates predict the identification of tens of compact binary coalescence (CBC) sources with electromagnetic (EM) counterparts detectable by ground-based optical telescopes (LIGO Scientific Collaboration, 2018). Specifically, amongst CBC sources, the most promising candidates to contain a visible electromagnetic (EM) counterpart are binary neutron star (BNS) and neutron star-black hole (NSBH) systems. The recent detection of the binary neutron star system GW170817 using the Hanford-Livingston-Virgo detector network (LIGO Scientific Collaboration et al., 2017a) demonstrated that a wealth of information can be gained from combining the electromagnetic and gravitational wave emission from the same source. For the 2017 detection, GW-EM information was used to identify the source host galaxy (LIGO Scientific Collaboration et al., 2017a), probe properties of the progenitor and remnant (LIGO Scientific Collaboration, 2017b, c), estimate the Hubble constant (LIGO Scientific Collaboration et al., 2017d), constrain possible models for the merger and emission (LIGO Scientific Collaboration, 2017e), and study the r-process nucleosynthesis processes resulting from the collision of the two neutron stars (Smartt et al., 2017; Coughlin et al., 2018b). The 2017 BNS detection also marked the first confirmed detection of a kilonova, the bright, infrared or UV emission hypothesized to result from the r-process nucleosynthesis occurring during a BNS or NSBH merger (LIGO Scientific Collaboration & Virgo Collaboration et al., 2017). This study focuses on the problem of identifying optical counterparts to GW or GRBs localized to large sky-error regions.

Due to the proximity of the binary neutron star system to Earth, the gravitational-wave source was localized to a sky area and volume of 28 sq. deg. and 380 Mpc3 respectively, making it the most well-localized amongst all past gravitational wave detections (LIGO Scientific Collaboration et al., 2017a). The fact that the source was localized to the nulls of Virgo’s antenna pattern significantly reduced the localization area from an ordinary two-detector localization (LIGO Scientific Collaboration et al., 2017a). IM2H, the observing team that first detected the optical counterpart to GW170817 about 10 hours after the gravitational wave detection, employed a strategy of targeting known galaxies within the source’s localization volume (Coulter, D. A., 2017; LIGO Scientific Collaboration & Virgo Collaboration et al., 2017). Though we anticipate that some future BNS detections will be as nearby and as well-localized as GW170817, not all the GW detections of BNS systems will have a bright optical counterpart. In fact, only a small fraction of sources will be optimally oriented, maximizing the GW amplitude and favoring the detection of a coincident on-axis afterglow (Metzger & Berger (2012); Petrillo et al. (2013); Barnes & Kasen (2013); Kasen et al. (2015); Metzger (2017)). Thus, we do not anticipate the search for optical counterparts to BNS sources to be as straightforward as it was during the GW170817 detection.

The approach for detecting counterparts to GW events differs based on the wavelength of the emission, due to the fact that each type of emission is observable for a different duration of time. Optical counterparts to BNS and NSBH sources, which could either be kilonovae (Metzger, 2017), or optical afterglows (Ascenzi et al., 2018; Ghosh & Bose, 2013), fade rapidly, and could last anywhere from hours to days after the gravitational wave source is detected (Metzger & Berger, 2012). Using optical telescopes to observe the source as early as possible after the GW detection will maximize the information gained from the source spectrum. During the upcoming third GW observing run, the two advanced LIGO (LIGO Scientific Collaboration, 2014) interferometers and the advanced Virgo (Virgo Scientific Collaboration, 2014) interferometer will be “online” and taking data for one year. Based on the expected number and sensitivity of gravitational wave detectors detecting BNS and NSBH sources in this Advanced Detector Era (ADE), GW localization sky error regions are expected to span a few tens to a few hundreds of square degrees (Nissanke et al., 2013; Singer et al., 2014), while most wide-field optical observatories following up GW sources tend to have fields-of-view of less than a degree to a few square degrees, requiring an extensive search within the sky localization region in order to locate the optical transient (Singer et al., 2014; Cornish & Littenberg, 2015; Essick et al., 2015; Klimenko et al., 2016; Rana et al., 2016). Modern optical surveys with large FOV telescopes such as ZTF (Bellm, 2014), ATLAS (Shanks et al., 2015), Pan-STARRS (Chambers et al., 2016), and LSST (Ivezic et al., 2008) are considered ideal in this scenario; however, implementing optimized search strategies for smaller-FOV optical telescopes will increase the odds of identifying a coincident optical counterpart to a given BNS/NSBH trigger by strengthening the overall telescope follow-up network. In the remainder of the paper, our discussion of “optical telescopes” will primarily concern telescopes for which a scheduling strategy is most relevant (i.e. FOV $<$ 10 deg2), though the same methods could be applied to optimize the scheduling of larger-FOV telescopes.

Optical telescopes usually employ either the galaxy-targeted or the tiling strategies to search for GW counterparts. The former strategy requires identification of all of the galaxies within the sky error region (and therefore, a near-complete galaxy catalogue in the vicinity of the source), and determines the probability of the galaxy containing the GW trigger, based on its mass or luminosity. Then, ranking the galaxies in order of their properties, the algorithm will provide the telescope with the coordinates of each of the galaxies within the region and the order in which the telescope should point at them (Gehrels, 2016).

The tiling strategy has four main steps: placing tiles, allocating time, scheduling, and evaluating efficiency (Coughlin et al., 2018b). Amongst different tiling methods described in Coughlin et al. (2018b), the ranked tiling method (Ghosh et al., 2016) that we employ in this paper involves dividing the GW sky error region into “tiles” that are the size and shape of the field-of-view of the observing telescope. After the placement of tiles, the next step is to perform time allocation. Salafia et al. (2017) and Chan et al. (2017) propose time allocation based on counterpart lightcurve models, while Coughlin & Stubbs (2016) derive scaling relations for time allocation based on GW likelihood and galactic extinction. The two algorithms we present in this work use two different methods of allocating exposure time to tiles. Once exposure time has been allocated to each tile, one can schedule the tiles for observation. Scheduling involves running an optimization algorithm on the tiles to maximize probability and patch area coverage, and generating a list of the coordinates of the observable tiles and the order in which the telescope should observe them (Rana et al., 2016). Finally, one can evaluate the process by running simulations and determining the amount of GW probability acquired by tiling the patch. Our optimization algorithms combine the tiling and the galaxy-targeted strategies.

Past work on optimizing telescope scheduling (Rana et al., 2016), uses the tiling method for the purpose of optical follow-up of GW transients, and describes three main algorithms for ranking tiles within the GW sky error region in order to generate an observing sequence. They are outlined in brief, below:

• •

greedy algorithm - ranks tiles in descending order according to probability; doesn’t account for setting tiles

• •

setting algorithm - selects the highest probability tiles in each setting window, starting from the earliest setting tiles

• •

optimal algorithm - reorders the observing sequence generated by the setting algorithm to observe the highest probability tiles first, without losing setting tiles

Rana et al. (2016) demonstrates that while the optimal and setting array algorithms have equivalent performance in terms of probability and area coverage within the patch, the optimal array provides an advantage over setting array as it schedules higher probability tiles to be observed first. This work builds on the work of Rana et al. (2016) by considering the effects of slew and airmass on the existing algorithms, and optimizing over both parameters.

Slew is the process of rotating a telescope to observe different regions of the sky. The slew angle between two tile observations is equivalent to the solid angle between two points on a spherical surface, simply given by the spherical law of cosines as:

[TABLE]

Here, alt1 and alt2 and $\Delta\rm az=|az_{2}-az_{1}|$, correspond to the two altitude coordinates and the difference between the azimuthal coordinates of each of the points on the sky, assuming a celestial sphere Earth-based coordinate system. The previous optimization algorithms neither take into account the effect of telescope slewing within the observing sequence nor factor in the time spent slewing between individual tile exposures. Thus, in this paper, we demonstrate the effect of slewing in the existing algorithms, and propose an alternate algorithm that minimizes the slew between tile observations without loss of tiles or probability.

We also calculate and compare the cumulative airmass amongst different algorithms for various patch-observatory configurations. Airmass is the path length for light from a celestial source to pass through the atmosphere; near the horizon, where light is attenuated by scattering and absorption, the airmass is at its maximum, while at the observatory’s zenith, it is unity. However, many of the tiles scheduled to be observed with high airmass will require a long exposure time in order to resolve the source, while tiles closer to the zenith could be observed within a brief exposure time. Thus we implement a modification to the the optimal algorithm that will minimize the airmass at which tiles are observed and adjust each tile exposure time based on tile airmass.

In addition to the algorithms presented in Rana et al. (2016), we discuss the following algorithms in this paper:

• •

modified optimal algorithm - this algorithm is a modification of the optimal algorithm mentioned earlier that factors in the slew time into the observing schedule

• •

slew-optimization algorithm - we propose a new algorithm to optimize over slew and setting constraints to acquire the maximum probability in the patch

• •

airmass-weighted algorithm - we propose another new algorithm to optimize over the airmass and setting constraints to acquire the maximum airmass-weighted probability in the patch

Like in the previous algorithms, all of the algorithms discussed in this paper account for ground-based visibility constraints and only schedule tiles that are above the horizon. We run simulations at the locations of the optical observatories GROWTH-India (GROWTH Collaboration, 2019a, b) and ZTF to systematically compare the performance of the slew optimization and the airmass-weighted algorithm with that of the modified optimal and greedy algorithms.

2. Description of Problem

While previous methods operated under the assumption that the slew time was accounted for in the tile exposure, our calculations demonstrate that the time the telescope spends slewing can cut significantly into the time available for observing tiles in the greedy, setting, or optimal sequence, preventing the telescope from covering all of the tiles provided by the optimal array. Depending on the shape of the patch, the telescope FOV, and the slew rate, the telescope can spend anywhere from 10s of minutes to over an hour in slewing to tile patches of a few hundreds of square degrees. In general, this problem of slew times exceeding a tile exposure resulting in loss of tiles is relevant for telescopes with camera readout times that are less than the average slew time, and telescope-patch configurations that require several telescope pointings in order to cover the 95 percent credible region. For the rapidly fading transients that optical telescopes seek to follow up, the loss of even a few minutes could make the difference between detection and non-detection of the EM counterpart of a GW source (Rana et al., 2016).

In order to determine how tiles lost to slewing affect the overall performance of our existing algorithms, we run the optimal algorithm using the FOV, location, and other properties of the GROWTH-India telescope (see Table 1) to tile a 94 square degree patch. First we run the algorithm without accounting for the slew time, and the second time, the algorithm calculates the overall time spent slewing, and removes tiles appropriately from the observing sequence. We refer to the optimal algorithm with slew accounted for as the modified optimal algorithm. Neglecting to account for slew time within the telescope scheduling could negatively affect the overall algorithm performance. If the total time spent slewing exceeds a single or several tile exposures, one or more tiles could be scheduled after they have already set. When we account for slewing, we omit these already-set tiles; as a result, accounting for slew demonstrates that the total probability acquired at the end of the observation is, in reality, lower than predicted by the original code. As demonstrated by Figure 1, in certain cases, the original optimal algorithm grossly overestimates the telescope’s ability to tile the patch when it does not account for time lost due to slewing. The comparison plots of the optimal algorithm tile coverage and cumulative probability are shown in Figures 1 and 2.

One of the limiting factors in being able to identify potential transient counterparts to gravitational wave events is the ability to resolve faint and distant sources. Depending on the size of the telescope field-of-view, within a given tile, there could be multiple galaxies, each with varying distance to Earth. In order to conduct a less-biased search, we choose to adopt a strategy such that the faintest known galaxies within the tile are resolvable. Furthermore, an optical survey is likely to miss even bright sources within tiles that are about to set, as the effects of atmospheric extinction attenuate the light from the source, making it more difficult to observe. For these reasons, when generally conducting optical observations, each target imaged is allotted a different exposure time. Therefore, the tile exposure time, though previously uniformly assigned to all tiles, must change depending on the most distant galaxies within the tile and the tile airmass.

3. Methods

3.1. Slew Optimization

Past work on slew optimization involved the radio follow-up of blazars (Max-Moerbeck, 2014), as the several hours large radio interferometers can spend slewing cuts significantly into the overall source observing time. In the case of radio follow-up, because the main constraints imposed are the rising and setting times of the patch, optimizing the slew involved applying a slight modification to the well-known “traveling salesman” algorithm to minimize the slew path between the tiles. However, as additional factors such as the telescope FOV, sunrise, and sunset are pertinent for optical follow-up, we cannot simply apply the traveling salesman algorithm on the entire patch to minimize the slew. Instead, we characterize the patches based on the patch visibility, determine whether slew optimization will be possible, and then apply the algorithm on them, with the setting and rising constraints imposed. In general, slew optimization will only be possible for cases in which all tiles in a given patch rise above the telescope’s horizon. As with the optimal array algorithm (Rana et al., 2016), our starting point is the setting array algorithm. We sketch the process of slew optimization for optical follow-up of GW transients as follows:

Convolve a galaxy catalogue with the GW probability distribution on the patch. We modify the tile probability based on the total mass of galaxies contained within each tile, eliminating tiles from the patch that do not contain any galaxies. This step is optional for the slew-optimized algorithm. 2. 2.

Run the setting array algorithm (Rana et al., 2016) on the patch. The setting array algorithm provides the set of maximum probability tiles within the available observation time. 3. 3.

Calculate the total slew angle and airmass of the tile array. We assign the variable texpo, a single tile exposure time, such that it will include the slew time of any tiles under 20 degrees when initially selecting tiles for the optimal sequence. The algorithm operates under the assumption that all slews less than 20 degrees will take the same amount of time as a 20 degree slew, as it is difficult for telescopes to maintain exact slew times for small slew angles ($\sigma_{\rm slew}\leq 20$), due to the acceleration and deceleration of the telescope. Therefore we allot a time gap such that $\rm t_{\rm gap}=\rm t_{\rm expo}+\rm T_{20}$, where T20 is the approximate time taken for slews less than 20 degrees. Here, $\rm T_{20}=(20\deg)/(v_{\rm slew})$, where $\rm v_{\rm slew}$ is the telescope slew rate. 4. 4.

Calculate amount of gained time and use it to acquire additional probability. For cases in which slew optimization is possible, after running the traveling salesman algorithm on the patch, we apply the time gained from optimization to search for additional unobserved tiles. We define the gained time as follows:

[TABLE]

where the first term is the difference between the accounted-for slew time and the exact amount of time spent slewing, and the second term records difference between time spent performing slews greater than 20 degrees, before and after slew-optimization. 5. 5.

Re-shift the tiles by their exact slew times. Because the slew optimization algorithm makes use of the gained time to observe additional visible tiles when possible, we space the tile observation times by their exact slew times, accounting for the time spent slewing between successive observations within the schedule itself. Both our modified optimal and slew-optimization algorithms space tile observations by the slew times between tiles such that the overall slewing time in the final schedule does not exceed a single tile exposure time. 6. 6.

Remove already-set tiles. The loss of total coverage and cumulative probability after accounting for slew arises solely from the fact that time-shifting the tile observations by their slew time results in a number of tiles being scheduled after they have already set. Thus we only include tiles in the schedule where the observation time is earlier than the tile set time. For tiles scheduled before they have risen, we shift their observation times by one tile exposure. This process is performed iteratively throughout the optimization. Because we have included the slewing time within the observation time, some tiles that could previously be observed using the original optimal algorithm will be lost due to setting.

3.2. Airmass and Exposure Time

Our objective is to maximize the total probability coverage based on the setting time and the airmass of the tiles on the GW localization. The tiles on the localization move as the Earth rotates, and the airmass of the tiles change as the altitude of the tiles change with time. The airmass of a given tile changes at each observation time. Therefore, we propose an algorithm to maximize the total probability coverage over the setting and airmass of the GW localization. We use the Hungarian optimization algorithm (Munkres et al., 1957) to get the optimal solution for this problem. In our algorithm, we define a new probability called the airmass-weighted probability. The airmass-weighted probability of a tile is the ratio between the galaxy-gw probability and the airmass of that tile. Although the galaxy-gw convolved (or GW) probability is fixed at all times for each tile in the localization, the airmass-weighted probability of each tile varies with time. We determine the airmass-weighted probability for each tile at all different times, using the time-dependent airmass, and use our algorithm to select the schedule of tiles that will maximize the airmass-weighted probability. We demonstrate the algorithm point wise below:

Convolve the GW probability with a galaxy catalogue. We modify the GW probability based on the total mass of galaxies contained in the given localization. We make a grid of tiles based on the convolved probability, eliminating those tiles that do not contain any galaxies. This is an optional step in our algorithm. 2. 2.

Account for airmass. The airmass is given by

[TABLE]

where $am(t)$ is the airmass as a function of time and $alt(t)$ is the altitude in degrees as a function of time and $alt(horizon)$ is the altitude of the observing horizon of the telescope. The airmass-weighted probability of a tile is inversely proportional to its airmass:

[TABLE]

where $p_{gw}$ is the GW (or galaxy-gw) probability, $p_{amw}$ is the airmass-weighted probability, and $Norm$ is the normalization constant chosen such that the total airmass-weighted probability is unity. Our strategy here is to devote less telescope time to tiles at higher airmass, so that low-airmass high-probability tiles are prioritized for observation. At any given time, if the altitude of a tile becomes less than the observing horizon, we make the airmass value infinity so that the airmass-weighted probability will become zero. Using Eqns. 2 and 3, we account for the setting and rising of all the tiles in the localization.

To illustrate how the tile airmass affects pamw, we can consider an example where two tiles with the same galaxy-gw probability are at different altitudes. If one tile is at an altitude of 30 degrees from the horizon and the other is located at the zenith, the first tile will have double the airmass of the tile at the zenith, so the airmass-weighted probability of that tile will be half of the airmass-weighted probability of the tile at the zenith.

We calculate $p_{amw}$ for every tile for all time steps from the observation start time to the observation end time. We make a table, where the rows represent the time steps (one time step is equivalent to one exposure) and the columns represent the tiles in the localization. One element in the $ith$ row and $jth$ column is the airmass-weighted probability of finding the source at the $jth$ tile at the $ith$ time step. 3. 3.

Use the Hungarian algorithm to maximize the airmass-weighted probability. If there exists more than one optimal solution, we choose the solution where higher probability tiles with lower airmass are scheduled earlier. 4. 4.

Allocate exposure time. For a galaxy targeted hunt, we adjust each tile exposure time based on the most distant galaxy within that tile, as the source’s flux decreases as its distance squared:

[TABLE]

where $t_{exp}$ is the adjusted exposure time of the tile, $t_{exp0}$ is the starting exposure time assigned to all tiles, d0 is the maximum distance within the 3D GW localization within the tile, and dgal is the distance to the farthest galaxy in the tile. If the telescope FOV is large enough, all tiles will contain both nearby and distant galaxies such that the required exposure time to resolve distant galaxies in each tile will average out.

Then, we apply a second adjustment to the tile exposure time based on its airmass:

[TABLE]

where am(t) is the tile airmass, and texp0 could either be the starting exposure time or the galaxy-adjusted exposure time for the tile, depending on whether the galaxy convolution step was performed. In allocating our exposure time, we do not account for the fading counterpart lightcurve.

The starting exposure time is a user-defined input into our code. One should set the starting exposure time by first selecting one tile in the localization, and then determining the amount of exposure time needed in order to observe a source at some desired magnitude out to the maximum distance within the tile (based on the 3D GW localization). Our code assigns that input exposure time as the initial exposure time for all tiles, and then applies the adjustments described above. In our simulation, we use the highest probability tile within the localization to set our starting exposure time.

Combining the airmass-weighted method with the slew optimization method is very complicated, as the airmass-weighted method breaks the order determined by the slew optimization algorithm that will minimize the slew between tiles. For this reason, we do not attempt to combine the slew and airmass optimizations. Instead, we show our results for the airmass-weighted algorithm separately from the slew optimization method, for which we only adjust the exposure times of the tiles appropriately in the finalized sequence. In this section, the two new methods we have presented are to be considered as alternative methods by which to schedule observations.

The question of which algorithm to use in a given situation now becomes relevant. The most direct way to determine which scheduling method is most suitable for a given telescope-patch configuration is to run both algorithms on the patch, and choose the algorithm that covers the most GW probability. On a standard 2.4 GHz processor, the slew-optimization algorithm takes $\sim$ 3-5 minutes to complete and the airmass-weighted method takes less than a minute to complete, so running both algorithms before the start of the observation does not present any significant overheads in time. The cumulative GW probability acquired is directly comparable between both algorithms, as the GW probability relies only on the GW skymap to inform about the true location of the source.

For the galaxy targeted search, an additional complication is the incompleteness of galaxy catalogs such as CLU and Glade (Cook et al., 2017; Dalya et al., 2018) out to distances larger than $\sim$ 200 Mpc. However since most optical telescopes are only sensitive out to about 200 Mpc, this should not significantly bias the results.

4. Results

In this section, we use both case studies and simulations to test the performance of our two algorithms. For the purpose of this study, the telescopes we select (GROWTH-India, HCT, Swope, and PTF) are wide-field optical telescopes that have small to medium FOV sizes and are apt for rapid follow-up of transient events. Table 1 shows the locations and parameters of each of the telescopes used.

We revisit the example referred to earlier in Sec. 2 to see how the slew-optimization algorithm performance differs from that of the modified optimal algorithm for the GROWTH-India telescope. We can compare the cumulative probability in each case using Figure 3. Without taking slewing time into account the telescope is capable of covering all of the tiles in the patch, covering a cumulative probability of 93.7% (see Fig. 1). On the other hand, the modified optimal algorithm can only cover a probability of 82.1%. The slew optimization algorithm is designed such that it cannot perform worse than the modified optimal algorithm, in terms of probability and coverage; for the cases in which the slew cannot be further optimized, the algorithm will return the same tile schedule as the optimal algorithm. In this case, the slew-optimization algorithm acquires about the same amount of probability (82.1%) as the modified optimal algorithm, though both algorithms far surpass the amount of probability covered by the greedy algorithm. The key advantage of the slew-optimized algorithm is demonstrated in Figure 4 by the cumulative slew angle being reduced from above 600 degrees to about 140 degrees as a result of slew-optimization. The time saved in slewing aids in more rapid observation of the patch, which is important especially when there are multiple targets to observe in one night.

By comparing cumulative slew angles between the slew-optimized array and the modified optimal array, we can see that the curve shapes differ quite drastically. The optimal method schedules tiles to be observed with small slew jumps between observations, resulting in a relatively smooth increase in overall slew angle with tile number. On the other hand, the slew-optimized method is mostly linear with two jumps around the 20th and 130th tile in the schedule. In effect, the slew-optimization method attempts to schedule as many tiles consecutively as possible before jumping to a different region of the patch. Figure 4 displays the behavior of the slew optimization algorithm as expected.

The telescope slew rate also influences the amount of improvement the slew optimization algorithm has over the modified optimal algorithm. Table 2 indicates that slew optimization is vital for slowly slewing telescopes that are more likely only to cover a small portion of the patch within the time constraint due to time spent slewing. Faster slewing telescopes tend to have a larger cumulative slew angle at the end of the slew-optimized schedule because they can afford to perform larger slews between tiles than slowly slewing telescopes, within the time constraint. Table 2 shows that even for telescopes with fast slew rates of 10 deg./s, slew optimization can reduce the overall time the telescope spends slewing.

In a similar fashion we check how the telescope FOV influences the performance of slew-optimization compared to modified optimal. We vary the FOV of the telescope, assuming a slew rate of 2 deg/s, and fixing the location and other parameters to that of Swope’s. Table 3 demonstrates the relative improvement of the slew optimization algorithm over the optimal algorithm; we note that there is more improvement for smaller FOV sizes. Therefore, the problem of slew optimization is especially relevant for smaller FOV telescopes with slow slew rates attempting to tile large GW sky error regions.

Using a different set of patch and telescope parameters, we now perform a comparison between the airmass-weighted and greedy algorithms, observing a marked improvement when using the airmass-weighted algorithm. Figure  5 shows the cumulative probability coverage by the greedy and airmass-weighted methods for a localization by the PTF telescope. Around tile 35 the optimal algorithm deviates from greedy, appearing to acquire probability more gradually at first, but eventually surpassing the greedy algorithm in tiles and probability covered. The greedy method covers 32.8% and airmass-weighted method covers 36.5% probability. We do not account for airmass in the greedy method. When comparing Figure 5 to the cumulative probability acquired by the original optimal algorithm, we observe that the curve corresponding to the optimal algorithm has the same shape as the airmass-weighted algorithm, but deviates from the greedy curve earlier in tile number. This indicates that the airmass-weighted algorithm prioritizes observing higher probability tiles earlier on in the observing schedule.

In order to test the performance of our airmass-weighted algorithm, we generate a distribution of GW sky localizations for binary neutron stars, assuming that the sources are distributed uniformly in volume. As our study pertains to the ADE, in which we expect our detectors to have improved in sensitivity over the current state of GW detectors, the sky patches we use are distributed between areas of 15-225 deg2. Though our algorithm can be run with various different input parameters for the FOV, coordinates, altitudes, and other observatory-specific parameters, we choose the location and parameters of the GROWTH-India telescope, changing the original field of view to a 1$\arcdeg\times 2\arcdeg$ square field of view to determine in which cases the airmass-weighted algorithm shows an improvement over the greedy algorithm. We keep the observing horizon at an altitude of $30\arcdeg$ (airmass $\sim 2$), and the observations are performed when the sun is below $-12\arcdeg$ altitude. We assume that the luminosity of the optical counterpart does not change with time, which fixes an exposure time of 300 seconds to see a source at a distance of 100 Mpc. The exposure time changes for different localizations based on their average distance. We ran the airmass-weighted and greedy methods on 900 of these 3D-localizations for the GROWTH-India optical telescope. 616 out of 900 localizations were visible from the GROWTH-India observatory; the localizations at the far south were not visible.

Figure 6 presents a scatter plot of the covered total probability by the airmass-weighted and greedy methods. The x–axis and y–axis display the probability covered by the greedy algorithm $(P_{gr})$ and probability covered by airmass-weighted algorithm $(P_{amw})$ respectively. The dashed line represents equal coverage by both of the methods. Most of the points fall above the dashed line, indicating that in most cases, $P_{amw}>P_{gr}$ regardless of the amount of probability acquired by the greedy method. In Figure 8, we compare the total areal coverage by the airmass-weighted method with the greedy method. Here, the x–axis is area covered by greedy algorithm, $(A_{gr})$, and the y–axis is the difference between the area covered by airmass-weighted algorithm and the greedy algorithm, $(A_{amw}-A_{gr})$. Each additional tile acquired by the airmass-weighted method over greedy adds an area of 2 square degrees because of the 2 sq. deg. field-of-view we use for these simulations. The general trend we observe from Figures 6 and 8 is that the improvement in areal or probability coverage between the airmass-weighted and greedy algorithms increases with increasing Pgr or Agr. This is because when patches are visible for a longer duration of time, both algorithms can cover a larger area of the patch, and the airmass-weighted algorithm can determine a more optimal solution to maximize the airmass-weighted probability.

To demonstrate the relative improvement in the total probability coverage between the two methods, we show Figure 7, in which the x–axis represents the probability covered by the greedy algorithm, $P_{gr}$, while the y–axis represents the difference between the probabilities covered by the airmass-weighted method and the greedy method, expressed as a percentage of the greedy probability coverage $((P_{amw}-P_{gr})/P_{gr})$. When the probability covered by the greedy method is less than $\sim 10\%$ of the total GW probability in a sky patch, then we see that the airmass-weighted method can cover more than double the probability acquired by the greedy method in the same patch. However, even when the greedy method covers a large percentage ($\sim$80%) of the GW localization, the airmass-weighted method is capable of covering up to 10 % more probability than greedy, which is consistent with the earlier trend we observed. In short, our results from simulations run to test the airmass-weighted method are a strong indication of the robustness of the algorithm.

5. Discussion

5.1. Summary

The optimal algorithm presented in Rana et al. (2016) provides an optimized telescope scheduling method for observing large sky-error regions that accounts for several observational constraints including the setting and rising of tiles, sun, moon and telescope time constraints. In this paper, we demonstrate the need to improve this algorithm by optimizing over slew and airmass, as the time spent slewing and additional exposure time required to resolve high airmass tiles detract from the overall time available for imaging the patch. We present two algorithms - airmass-weighted, and slew-optimized - and describe criteria for choosing which algorithm to use to schedule observations. One important modification we make to the original optimal algorithm is adding in an optional step in both algorithms to convolve the GW patch with a galaxy catalogue, combining the tiling and galaxy-targeted search strategies. We demonstrate that both algorithms presented in this work are more realistic that the original optimal algorithm, and improve over the greedy algorithm not only in terms of cumulative airmass and slew, but also in overall probability acquired.

5.2. Caveats

Many of the caveats of the original optimal algorithm remain true for our newly proposed algorithms, as do the procedures for overcoming these caveats. We briefly summarize these caveats below (see Rana et al. (2016) Sec. 4.2 for a more detailed discussion).

We assume that the source light curve stays flat throughout the duration of the observation. GW170817 is the only BNS merger event observed so far, and its optical counterpart was first seen $\sim 11$ hours after the merger. Therefore, we do not yet have any observational information about the evolution of the light curve during the 11 hours immediately following the merger. Theoretically, the optical and infrared light curves of the counterpart of the BNS merger might vary depending on the neutron star masses and equation of state. Similarly, no NSBH mergers have been observed conclusively either. Consequently, we have chosen to be agnostic about this aspect and opted for a flat light curve here. Future studies may extend this work based on any new information that arises from anticipated merger observations involving neutron stars in the coming years.

We do not account for cloudiness in our scheduling since it is outside the scope of this work. This could be a potential weakness in most telescope scheduling algorithms – the cloudiness constraint should be addressed in a future work.

As mentioned in the previous work, synoptic surveys will often image based on a pre-defined grid for comparison of previously imaged fields; again, here offsetting our grid tiles (placed based on the maximum GW probability) to match the pre-defined grid should not significantly affect the performance of either algorithm.

The true tile airmass, though parameterized as purely a function of the tile altitude, depends on atmospheric visibility as well as the color of filter used to conduct the observation. We do not account for “seeing” or filters into our airmass calculation as it would make our algorithm very telescope-specific; instead we use the formula given in Section 3.2 to calculate the airmass of the tiles.

As proposed in (Rana et al., 2016), one can adjust the final tile observation time to account for counterparts that may fade on timescales shorter than a day. We anticipate, however, that in most cases optical observatories tiling GW or GRB localization patches will continue observation until there is confirmed non-detection. As evidenced by the first joint GW-EM detection in 2017 (LIGO Scientific Collaboration & Virgo Collaboration et al., 2017), it will be difficult to determine exactly the timescale on which an optical counterpart will fade. However, the airmass-weighted  algorithm addresses this in part by allocating longer exposure times to tiles containing more distant galaxies.

Acknowledgements: This work made use of the Python libraries Numpy and Matplotlib. It also made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013, http://www.astropy.org). We would like to thank Patrick Brady, Leo Singer, Varun Bhalerao, G. C. Anupama, Om Sharan Salafia, and Shaon Ghosh for helpful discussions.We would also like to thank Michael Coughlin for carefully reading the manuscript and making useful comments in LIGO P&P review (https://dcc.ligo.org/LIGO-P1900019). This work was supported in part by the Navajbai Ratan Tata Grant.