Searches for Modulated {\gamma}-Ray Precursors to Compact Binary Mergers in Fermi-GBM Data

GW170817 is the only gravitational-wave (GW) event, for which a confirmed {\gamma}-ray counterpart, GRB 170817A, has been detected. Here we present a method to search for another type of {\gamma}-ray signal, a {\gamma}-ray burst precursor, associated with a compact binary merger. If emitted shortly before the coalescence, a high-energy electromagnetic (EM) flash travels through a highly dynamical and relativistic environment, created by the two compact objects orbiting each other. Thus, the EM signal arriving at an Earth observer could present a somewhat predictable time-dependent modulation. We describe a targeted search method for lightcurves exhibiting such a modulation, parameterized by the observer-frame component masses and binary merger time, using Fermi-GBM data. The sensitivity of the method is assessed based on simulated signals added to GBM data. The method is then applied to a selection of potentially interesting compact binary mergers detected during the second (O2) and third (O3) observing runs of Advanced LIGO and Advanced Virgo. We find no significant modulated {\gamma}-ray precursor signal associated with any of the considered events.


INTRODUCTION
Multimessenger astronomy started with the detection of the core-collapse supernova 1987A in both the electromagnetic (EM) and neutrino channels (Arnett et al. 1989). Since then, only one other unambiguous example of astrophysical event heralded by different messengers occurred: the simultaneous detection of the gravitational waves (GWs) from the binary neutron star (BNS) coalescence GW170817 (Abbott et al. 2017), and of several EM counterpart signals: the high energy photons of GRB 170817A (Goldstein et al. 2017;Savchenko et al. 2017), the ultra-violet, optical and infrared radiation of the kilonova AT 2017gfo (Coulter et al. 2017;Smartt et al. 2017) and the X-ray, optical and radio afterglow of the γ-ray burst (Lamb et al. 2019;D'Avanzo, P. et al. 2018). Additionally, there is a convincing claim for the coincident detection of the high energy neutrino IceCube-170922A and the multi-wavelength EM radiation coming from the γ-ray blazar TXS 0506+056 (Aartsen et al. 2018).
In the observable Universe, compact objects such as neutron stars and black holes are often found in pairs, forming binaries. During the inspiral, they lose angular momentum and binding energy, by emission of GWs (Taylor & Weisberg 1982). This implies a narrowing of the distance separating the binary components, leading in some cases to a merger in less than a Hubble time (Phinney 1991). Both the frequency and the amplitude of the GWs increase as the merger approaches; during the seconds prior to the merger, the frequency of the GWs sweeps from tens of Hz to above a kHz. If the binary is close enough to us, the GW strain is above the sensitivity threshold of and thereby detectable by GW detectors such as Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015).
One of the EM counterparts to compact binary mergers is the flash of γ-and X-rays from the short γ-ray burst (Berger 2014;Kochanek & Piran 1993), lasting less than 2 s (Kouveliotou et al. 1993) and possessing an isotropic-equivalent energy up to 10 51 erg (D'Avanzo 2015). This short, high-energy EM signal is followed by a longer-lasting, less energetic radiation, the γ-ray burst afterglow (Fox et al. 2005) covering a broad range of the EM spectrum, from X-ray (Vietri 1997), through optical (Meszaros & Rees 1997) to radio (Paczynski & Rhoads 1993). The kilonova comprises the ultraviolet, optical and infrared radiation associated with the radioactive decay of heavy elements (Lattimer & Schramm 1974;Li & Paczynski 1998;Metzger et al. 2010;Kasen et al. 2017); it has a quasi-thermal (Kasen & Barnes 2019) and quasi-isotropic emission (Darbha & Kasen 2020), very different from the γ-ray burst afterglows which display large anisotropies (Beloborodov et al. 2011), and consequently require nearly-aligned observer positions in order to be detectable.
While the association between the merger of compact objects and some γ-ray burst related messengers, such as the γ-ray prompt emission, the afterglow and the kilonova, has been unambiguously highlighted by GW170817, the presence of other kinds of EM emission is still debated. One such example is the precursor activity to short γ-ray bursts. Troja et al. (2010) claim that up to 10% of short γ-ray bursts, detected by Swift (Barthelmy et al. 2005), possess EM precursors, lasting less than 1 s and whose starting emission might go back about 100 s before the main γ-ray burst. The existence of such precursor signals is also supported by Zhu (2015). The investigation by Li et al. (2021) concludes that the precursors have shorter duration than the prompt γ-ray emission, but seem to be produced by similar central engine activity. 16 precursors to short γ-ray bursts in Fermi-GBM data have been found by Wang et al. (2020), who infer comparable duration for the main and precursor emissions, and possible fits with blackbody, non-thermal cutoff power law models. Precursor emission is also motivated on a theoretical basis. Two popular models are the resonant shattering of the crust of a neutron star during the inspiral (Neill et al. 2021;Tsang et al. 2012) and the "black hole battery" model (McWilliams & Levin 2011). Theoretical work has also been done to identify possible features of an EM signal emitted during the inspiral, namely a modulation induced by the orbital motion (Schnittman et al. 2018).
The prompt γ emission, afterglows and kilonovae described above are currently expected to follow from the disruption of a neutron star and the formation of an accretion disk during the merger of the binary. However, whether a neutron star will actually disrupt before merging with the companion strongly depends on the properties of the two objects, in particular on their masses, spins and structure (Foucart et al. 2018). Although BNS mergers are generally always expected to radiate the whole plethora of EM signals observed with GW170817, this is far from being guaranteed for neutron star -black hole (NSBH) mergers. In addition, the prompt γ emission is expected to be detectable only for very specific orientations of the binary with respect to the observer. Hence, the idea of a premerger, precursor emission that does not require disruption and may not be strongly anisotropic becomes particularly interesting: it may be the only EM signal systematically emitted by NSBH mergers, even those involving non-spinning black holes more massive than ≈ 10 M .
Motivated by the above consideration, and by the present ambiguity regarding the possibility of precursor emission, we propose a method to analyze archival γ-ray data in temporal proximity to GW events associated with BNS mergers. Following the idea of Schnittman et al. (2018), the analysis aims at detecting pulsations in the γ-ray data having the same orbital phase evolution like the GW inspiral signal, during the last several seconds before merger. The method is an extension of a more generic search for γ-ray transients based on a likelihood approach, introduced by Blackburn et al. (2015) and further optimized in Goldstein et al. (2016) and Goldstein et al. (2019). We demonstrate how the existing and proposed methods recover simulated signals added to archival data, and we apply the proposed method to data around a selection of compact binary merger events detected by Advanced LIGO and Advanced Virgo.
We employ data from the γ-ray Burst Monitor (GBM) instrument on the Fermi spacecraft (Meegan et al. 2009). GBM is an ideal instrument to study rapidly-evolving γray counterparts to compact binary mergers. It is able to perform time-resolved spectroscopy of high-energy EM radiation, thanks to its twelve sodium iodide (NaI) and two bismuth germanate (BGO) scintillation detectors, covering an energy range from 8 keV to 40 MeV (Atwood 1994), and distributed around the spacecraft. It has a temporal resolution of 2 µs, suitable to study the rapid modulations we are interested in, and can discriminate between 128 energy ranges/channels. Its field of view is only limited by the Earth, and its observing time is mostly limited by passages through the South Atlantic Anomaly, implying it can witness an astrophysical transient event 75% of the time.
The paper is structured as follows: in Section 2, we discuss the possible physical mechanisms responsible for the γ-ray precursor activity, we introduce a modulated EM waveform model and present how it can be used to simulate data from a γ-ray detector. The statistical search method is presented in Section 3. The sensitivity of the search is presented in Section 4, and the results of the search around some of the LIGO-Virgo events appear in Section 5. Finally, Section 6 offers the conclusion.
2. MODULATED γ-RAY PRECURSORS 2.1. Physical mechanisms and signal models While the physical mechanism responsible for such a γ-ray precursor is still in question, several possibilities have been suggested in the literature. Palenzuela et al. (2013) and Most & Philippov (2020) show that magnetosphere interaction in a BNS can power EM radiation, prior to the main emission, with luminosities reaching 10 45 erg/s. Copious EM radiation, prior to the merger of a magnetized neutron star-spinning black hole binary, can also be emitted through the unipolar inductor mechanism (D'Orazio et al. 2016;D'Orazio & Levin 2013;McWilliams & Levin 2011;Palenzuela et al. 2011): indeed an electric circuit is formed, where the roles of battery, resistor, and electrical wires are played by the black hole, the neutron star and its magnetosphere, and the magnetic field lines; in this way, black hole rotational energy is extracted by the magnetic field and sent far away by means of powerful, collimated Poynting flux.
During the last orbits of a quasicircular inspiral, or a periastron passage in an eccentric or hyperbolic close encounter, a neutron star can experience tidal forces that excite some of its many oscillation modes, e.g. core, shear, crustal discontinuity modes (Lai 1994;Shibata 1994). This process, especially when resonant, can cause quakes and/or shattering of the neutron star crust, followed by the release of a huge amount of energy in the form of EM radiation (Suvorov & Kokkotas 2020;Tsang et al. 2012;Reisenegger & Goldreich 1994;Tsang 2013).
While traditionally γ-ray burst emissions are detected by searching for an excess of photons with respect to the background (Meegan et al. 2009;Kocevski et al. 2018;Burgess et al. 2016), in the present work we attempt to increase the sensitivity of the pipeline by hunting for well-characterized EM waveforms. Schnittman et al. (2018) propose a method to calculate the EM precursor lightcurve from a NSBH, where the surface of the neutron star is uniformly shining. A similar method is proposed in Haiman (2017) for the case of a supermassive black hole binary, which are future LISA (Amaro-Seoane et al. 2017) signal progenitors. In the preceding papers, the lightcurve is modulated by physical processes such as relativistic beaming and gravitational lensing. There is a parallelism to the compact binary coalescence GW waveforms as the EM lightcurve could be locked in phase with the GW signal. As the merger time approaches, the amplitude and the frequency of the presumed EM signal increase. Radio precursor lightcurves to compact binaries containing at least one magnetized neutron star are also suggested in Sridhar et al. (2021).
The association of GWs and γ-ray burst precursors may represent an unique class of multimessenger events, in the coming years. Indeed, the theoretically-motivated signals presented above are expected to be visible only for nearby events. As GW detections will be nearby too, this is a huge gain for this kind of work. In addition, NSBH systems, for which the mass ratio is too high, are not expected to generate neither short γ-ray bursts nor kilonovae because the neutron star component is swallowed by the black hole companion before being disrupted by the tidal field. On the other hand, these systems might power precursors of γ-ray burst, during the inspiral phase. Moreover, lightcurve models like the one proposed in Schnittman et al. (2018) are interesting because they give information about the potential γ-ray signal we want to detect, based on the GW detection.

Orbital-modulation model
An important feature of an EM lightcurve is the evolution of the brightness with time. In this work, given a binary inspiral, we are interested in the orbital-phasedependence of the luminosity. For a compact binary, the first Post-Newtonian expansion (Blanchet 2014) term of the orbital angular frequency evolution can be written as where M = (m 1 m 2 ) 3/5 (m 1 + m 2 ) −1.5 , G and c are the observer-frame chirp mass, the gravitational constant and the speed of light in vacuum respectively. t c and t are the merger and variable time measured in the observer frame. m 1 and m 2 are the binary component masses. We define the orbital phase, Φ orbit (t) ∈ [0, 2π], with origin at t = −30 s, by Φ orbit (t) = t −30 s Ω(x)dx mod 2π.
We also introduce the simplified lightcurve, a lightcurve which, besides the chirp mass, depends on two other parameters, namely θ peak and θ width . The normalized luminosity has the following expression: This expression corresponds to a compact binary emitting only during the orbital phase window centered at θ peak , with width θ width . In Figure 1, there is an illus- tration of a simplified lightcurve. Additionally, for comparison purposes, we show a more realistic (albeit still approximate) lightcurve which attempts to explicitly model the relativistic beaming and gravitational lensing effects described in Schnittman et al. (2018). Although we will not use this model for the analyses described later, more details of its construction are presented in Appendix A. Figure 1 shows that the blue simplified lightcurve and the more realistic orange lightcurve are similar. The choice of the simplified lightcurve for the present paper is motivated by the pronounced modulation feature, limiting the spread of incoming photons in regions where the light flux is negligible. While the exacerbated modulation of the simplified lightcurve is most likely unrealistic, the preference for this model is mainly motivated by our attempt to check the validity of the search method in the case of signals with optimistically large modulation.

Simulating a Fermi-GBM observation
We use Fermi-GBM's Time-Tagged Event (TTE) data 1 , which consists of a list of photons characterized by their arrival time (to microsecond precision) and energy channel. In our analysis, for simplicity and to increase the photon statistic, we coalesce the 128 possible energy channels into 8 main channels, larger in energy.
In order to simulate a simplified-lightcurve-like signal, one needs to convert the spectra into detector counts. To this end, one needs the detector response to radiation, characterized by the energy of its photons, the arrival direction and the amplitude of the lightcurve. This function includes two components: the response due to the direct radiation (Kippen et al. 2007) as well as the response due to the scattering from both the Earth's atmosphere and the spacecraft (Pendleton et al. 1999). The spectral and directional dependence was validated experimentally by ground-based calibration (Bissaldi et al. 2009). As explained in Connaughton et al. (2015), based on the numerical values of the detector response taken at 41,168 grid points (accounting for 272 sky directions), the response to any arrival direction is constructed by interpolation among the three closest grid points, obtained by Delaunay triangulation. Regarding the radiation energetics, we make use of three spectral templates, which we refer to as: soft (lowest energy), normal and hard (highest energy). For the soft and normal templates, we use the Band parameterized functions (Band et al. 1993), introduced in Connaughton et al. (2015). Thus the flux of photons having an energy in with (α spec , β spec , E peak ) equal to (−1.9, −3.7, 70 keV) and (−1, −2.3, 230 keV) for the soft and normal spectrum. With respect to the hard spectrum, we employ the comptonized template proposed in Goldstein et al.
where E piv is a constant and (α peak , E peak ) equals (−0.5, 1.5 MeV). However, the γ-ray bursts detected to date reveal a richer set of spectra than the three tem-plates presented here. This fact motivates our choice of considering injected lightcurves parameterized by the variable κ ∈ [0, 2], such that κ = 0, 1 and 2 correspond to the hard, normal and soft spectrum. For all the other values, we consider linear combinations of the three templates with the following weights: where w 0 , w 1 and w 2 are the weights attributed to the hard, normal and soft spectrum templates.
The detector output to an input flux of high-energy photons is captured by the Response Matrix. More precisely, the Response Matrix is a 14 × 8 array, where each row designates one of the 14 detectors and each column stands for an energy channel. Finally, each element of the Response Matrix represents the photon rate as a function of time. As the TTE data counts are assigned arrival times, we convert the photon rate function into a time histogram using a Poisson distribution. Moreover, for technical reasons, the data is binned. More precisely, the time is divided in intervals of size equal to 0.5 ms, and the photons found in the same interval/bin and belonging to the same energy channel are summed. According to the Nyquist-Shannon sampling theorem, such a binning allows the preservation of signal frequency components up to 1000 Hz. According to Equation 1, the orbital frequency is always lower than this upper limit, for all binaries where the heavier compact object weighs more than 1.4 M and for a binary evolution up to the last millisecond before the merger. Finally a simplified lightcurve injection is completely parameterized by the tuple (t c , m 1 , m 2 , f amp , t start , ∆t dur , θ peak , θ width , θ ra , θ dec , κ), where the variables f amp , t start , ∆t dur , θ ra and θ dec designate the lightcurve amplitude factor, the EM signal start time, the signal duration, the right ascension and the declination associated with the sky location. It is worth mentioning that the photon flux is obtained by the multiplication of the amplitude factor f amp with the expression of the normalized simplified lightcurve proposed in Equation 2.

STATISTICAL SEARCH METHOD
We propose a statistical framework which adapts the search methods presented in Goldstein et al. (2016); Blackburn et al. (2015). Hereafter, we refer to the original search as the generic targeted search, since it aims at detecting a generic transient excess of highenergy photons above the detector background associated with a particular target time, regardless of its tem-poral morphology. We refer to our modified method as the chirp targeted search instead, as it aims at detecting excesses of photons that exhibit repetitions locked in phase with the orbital evolution of a compact binary that is about to merge. The main difference with respect to the generic targeted search is that we apply its statistical formalism in the orbital phase space, instead of the time space. The transformation from time space to orbital phase space is provided, approximately, by the coalescence time t c and the component masses m 1 and m 2 inferred from the GW signal. For a fixed number of bins N bins , the orbital phase is split into N bins equal intervals, i.e.
N bins , 2π . The top panel of Figure 2 illustrates the conversion from the time space to the orbital phase space. Once this transformation is defined, we can rebin the photons registered by GBM (the TTE data) into the orbital-phase intervals: using each photon's arrival time t, we identify its orbital phase interval I k such that Φ orbit (t) ∈ I k .
Next, we need an estimate of the background photon rate in each phase interval, i.e. the rate of photons registered by GBM in the absence of a modulated transient. To this end, we first estimate the background rate over time using the unbinned Poisson maximum likelihood technique introduced in Goldstein et al. (2016): at a given time t 0 , the background photon rate λ max (t 0 ) is defined as the ratio between the number of photons N photons contained in a large enough time window of duration T (in our case T = 100 s) and the width of the window, i.e. λ max (t 0 ) = N photons T . Assuming the background can be described as a stationary Poisson process over the search interval, the uncertainty in its rate can be written as σ 2 λmax (t 0 ) = λ max (t 0 )/T . The time window is slid over the time range of interest; thus a background photon rate (assigned with standard deviation) is calculated for an array of times, and finally the photon rates (and their uncertainties) are interpolated over time. A chi-squared statistic χ 2 is computed in order to evaluate the quality of the fit. If the fit is poor, i.e. χ 2 is too large, for one of the GBM energy channels, that channel is excluded from the search. The background rate is then transformed to the orbital phase space to predict the rate in each phase interval.
Once the foreground photon histogram (background photon fitting) is calculated (estimated) for the [0, 2π] orbital phase range, we aim to search for a subset of adjacent intervals I k presenting an excess of photons. Such a behaviour is equivalent to saying that the binary, during an orbit, emits the majority of the radiation in a specific orbital phase window. This feature is characteristic to the simplified lightcurves.
A quantity combining information about both source and noise is the likelihood ratio, defined as where d, H 1 and H 0 are the observed data, the signal presence hypothesis, and the hypothesis of noise alone. As in Blackburn et al. (2015), the assumption of uncorrelated Gaussian noise allows us to write the preceding probabilities in the following way: In those previous formulas,d i = d i − < n i >, where d i and < n i > are the foreground and the estimated background photons. σ 2 ni and σ 2 di represent the variances of the background and the expected data (background plus signal). The variances appearing in Equations 7 and 8 are computed in the orbital phase space, and are obtained by the summation of time space variances. Lastly the time space standard deviations are calculated as in Blackburn et al. (2015). r i stands for the detectorenergy response, which depends on both the EM source sky location and spectrum. Finally s is the amplitude (measured by the Earth) of the signal. Moreover, for all these quantities, the index i designates a pair (detector, energy channel). Given that σ di , σ ni andd i are measured quantities, while r i is calculated for a sample grid, accounting for all possible locations, as explained in Kocevski et al. (2018), then the amplitude parameter s is the only variable over which the marginalization needs to be done. Thus, the expression of the likelihood ratio becomes Λ(d) = P (d|H1,s) P (d|H0) P (s)ds, where P (s) is the amplitude signal prior. Maximizing the likelihood ratio is the same as maximizing the logarithm of it, L(d) = ln Λ(d). As explained in Blackburn et al. (2013), ln Λ(d|s) is almost a Gaussian function with variance σ 2 ln Λ(d|s) = . Therefore the maximum of ln Λ(d|s) is reached for s best , obtained by means of the iterative Newton's method. The (k + 1) th step in the Newton's method consists in refining the k th estimate by using the analytic second derivative, i.e. s k+1 ≈ s k − ∂L/∂s ∂ 2 L/∂s 2 , while the initial guess is . We assume the same well-behaved prior like in Kocevski et al. (2018), i.e.
where γ prior = 2.5 and β prior = 1. The value of γ prior ensures a prior almost constant over a range of σ L , while the value of β prior translates in a luminosity distribution independent of distance. The log-likelihood ratio becomes L(d) = ln σ ln Λ(d|s) +ln 1 + erf s best √ 2σ ln Λ(d|s) +ln Λ(d|s best )+ In Equation 10, erf(x) = 2 √ π x 0 exp −t 2 dt is the error function. And finally, as explained in Blackburn et al. (2015), we calibrate the log-likelihood ratio (hereafter LLR) by subtracting the quantity L ref = β prior ln γ prior + (1−β prior ) ln σ ref . Here we do not care about the value of σ ref because β prior = 1, and so the term (1−β prior ) ln σ ref cancels.
Various kinds of transients commonly appear in Fermi-GBM data and produce large LLR values, despite being certainly unrelated to GW events. One such class is represented by high-energy cosmic rays hitting the NaI detectors, being responsible for long-lived phosphorescent light emission. By means of the technique introduced in Blackburn et al. (2015), we remove most undesirable high-LLR triggers from this class. A second class is represented by sharp photon-rate changes due to Fermi approaching the South Atlantic Anomaly (SAA). These are discarded instead as explained in Goldstein et al. (2016).
Both the generic targeted search and our chirp targeted search take as input a GPS time t c , which is the compact binary merger time, measured at Fermi. In this work, the generic targeted search is performed over the following exact timescales: 0.064 s, 0.128 s, 0.256 s, 0.512 s, 1.024 s, 2.048 s, 4.096 s and 8.192 s. It uses a time displacement of 64 ms for the four shortest timescales and a time displacement of factor 8 for the four longest timescales (e.g., the 8.192 s search windows are separated by 1.024 s). Thus, in total there are 2270 such windows. The positions of some search windows on the time axis are shown in the bottom panel of Figure 2. For the chirp targeted search, we shorten/extend the same search windows by a small amount in such a way that each window contains an integer number of orbits. This procedure is important in order to avoid artificially unequal number of photons in different I k intervals, which would produce artificially large LLR values. Additionally, in the case of the chirp targeted search, for a fixed (t start , ∆t dur ), searches are done over 47 subsets of adjacent intervals I k . More precisely, we consider subsets of any length in {1, 2, ..., N bins } and we use a phase factor of 2 (e.g. subsets of length 2, 3, and 6 are separated by 1, 2 and 3 intervals I k ). Finally the most significant trigger, i.e. with the highest LLR, is reported.

Null distribution of LLR
It is worth mentioning that a statistically significant trigger is not necessarily a GW-related signal, and sometimes not even the effect of an EM radiation intercepted by the NaI and/or BGO detectors. Despite the filtering strategy discussed in Section 3, large LLR spurious signals survive. For this reason, an empirical measure of the false alarm probability (FAP) distribution is extremely useful. In this paper, for a given log-likelihood ratio LLR 0 , the value of FAP(LLR 0 ) represents the fraction of noise events in which at least one set member has a statistical significance higher than LLR 0 . In Figure 3, we show the FAP distribution of LLR, obtained by running the search on 1000 random times, spread over the period going from April 1, 2019 to April 1, 2020. The random times have been chosen such that they are at least 30 s away from the SAA entrance/exit. This choice is motivated by the exclusion of those situations where the 30 s time windows would otherwise analyze nonscience times. Based on Figure 3, at least two important remarks should be made: (i) even though similar, the distributions corresponding to the chirp targeted search have higher LLRs with respect to the generic-targetedsearch distribution; (ii) the noise output triggers have higher LLR as the number of bins N bins increases. The remark (i) is an expected behavior, because the generic targeted search is included in the chirp targeted search. In fact, the evaluation of the statistical significance of all I k intervals together is equivalent to performing the generic targeted search. Regarding (ii), there are at least two reasons favoring this behavior: (a) for two bins numbers N bins,1 and N bins,2 , with N bins,2 > N bins,1 , such that N bins,2 is a multiple of N bins,1 , the chirp targeted search with setting N bins,1 is included in the chirp targeted search with setting N bins,2 ; (b) the higher the number of bins N bins , the less true is the approximation of Gaussian background noise in the high energy detectors, on the scale of a I k interval.

SEARCH SENSITIVITY
In this section, we test the sensitivity of our chirp targeted search and we compare it to the sensitivity of the generic targeted search. In the case of modulated γ-ray signals, in phase with GWs, we expect better performance for the chirp targeted search. In this subsection, we describe the simplified lightcurve injections in the GBM data. We consider signals with durations ∆t dur log-uniformly distributed in . This choice was found empirically in order to respect the following requirements: the lower and upper limits impose for the majority of the injected signals to have statistical significance right above the background LLR distribution; the dependence on ∆t dur causes signals with different durations to have similar LLRs. Finally, all the modulated signals considered in this study are injected at random times spread over one year period starting at April 1, 2019 in such a way that the merger time is always at least ±30 s away from the closest SAA episode.

Comparison with the generic targeted search
In this subsection, we simulate signals with (m 1 , m 2 ) equal to (10 M , 1.4 M ) and (1.6 M , 1.4 M ), and recover them with the chirp targeted search, where (m 1 , m 2 ) is fixed at the same value. We also apply the generic targeted search to these injection sets. For each injection, the most significant trigger is selected and a FAP is derived according to the results summarized in Figure 3. The cumulative distribution functions (CDF) of the FAP for these injections sets are illustrated in Figure 4. The first remark one should draw from Fig-ure 4 is that the chirp targeted search is more sensitive than the generic targeted search in the case of simplified lightcurves, as long as N bins > 1. The case of the chirp targeted search with setting N bins = 1 should be equivalent to the case of the generic targeted search. This relation is verified here, the small discrepancy between the two distributions being due to the different technical implementations.
The N bins dependence of the chirp targeted search sensitivity warrants some discussion. Firstly, Figure 3 suggests that the statistical significance of the noise triggers increases with N bins , which has as effect the degradation of the sensitivity with the augmentation of N bins . Secondly, the performance of the pipeline is expected to depend on the relation between the width of intervals I k and the orbital phase length of the chirping signals we want to detect. Following this idea, the sensitivity should increase with N bins , as long as the width of I k is higher than the orbital phase length of the recovered signal. Indeed, for a given EM chirp radiation, if the signal is included in only one interval I k , the smaller the I k the higher the signal-to-noise ratio, because the higher the percentage of I k where the foreground photon rate is above the background photon rate. Therefore, one should expect a compromise between the two regimes described above: an increase of the sensitivity with the number of I k intervals at low N bins , then a saturation followed by a degradation of the performance at high N bins . This behavior is verified in Figure 4. For N bins equal to 1, 5 and 15, the width of the interval I k is 360 • , 72 • and 24 • . Given that θ width = 10 • , such a signal could a priori be included in one interval I k , unless it is situated at the border of two adjacent I k intervals. From Figure 4, one can note an increase of the sensitivity between the cases N bins = 1 and N bins = 5. However, the performances of the pipeline seem to be quite similar in the cases N bins = 5 and N bins = 15. Figure 5 shows the statistical significance dependence of the chirp targeted search with the parameters of the simplified-lightcurve injections. The top left panel of Figure 5 proves that our choice for the ∆t dur dependence of the amplitude factor f amp puts on an equal footing the signals with different durations. One should note that an injection with f amp = 1 corresponds to an energy flux of 1 erg/s/cm 2 in the 50−300 keV band. Thus, the same panel indicates that a signal spread over around 1 s and emitting during 10 • orbital phase window, with a flux reaching 10 erg/s/cm 2 in the Fermi-GBM band, might be recovered by our pipeline with an important statistical significance (LLR > 100). The top right panel shows that the pipeline sensitivity is independent of the signal position in the orbital phase space. From the bottom left panel, one can note better performance of the chirp targeted search when the signal has a spectrum either close to κ = 0 (hard ) or close to κ = 2 (soft). This might be explained by the fact that although for the injections we have considered κ anywhere in the interval [0, 2], the search recovery is realized with only the three spectrum templates corresponding to κ ∈ {0, 1, 2}. Therefore, an injection with κ not near to 0 or 2, represents an heterogeneous weighted addition of quite different spectral templates, and which is more difficult to be recovered with only one template in between {soft, normal, hard }. Finally, the bottom right panel of Figure 5 illustrates the independence of the sensitivity with the position in the sky of the injected signal.

Impact of orbital-phase uncertainty
The properties of a compact binary merger inferred from the GW data always carry some uncertainty. In particular, there are uncertainties on the chirp mass, the merger time and the sky location. The chirp mass uncertainty will impact our knowledge about the orbital phase evolution. The uncertainties on the geocentric merger time and the sky location together reflect on the uncertainty of the merger time at the position of the Fermi satellite, which in turn affects our knowledge of the positions of the I k intervals. Concerning the chirp mass detection errors, the O1, O2 and O3 LIGO-Virgo observing runs showed that the uncertainties are below 0.1 M for BNS or NSBH like objects, and of the order of a few M for binary black hole mergers. Concerning the geocenter merger time detection, the uncertainty is of the order of 1 ms. However the sky locations of GW mergers are often poorly constrained, especially in the case of single-interferometer detections, which means that the uncertainty of the merger time measured at Fermi is of the order of the photon flight time between Fermi and the center of the Earth, i.e. ≈ 23 ms. Before applying our method to LIGO-Virgo events, we have to evaluate the impact of imprecise GW measurements on the sensitivity of our search.
The variation of the sensitivity with the uncertainty on the Fermi location merger time and on the chirp mass is illustrated in Figure 6. In order to investigate the impact of the simultaneous inaccuracies of the merger time and chirp mass measurements on the recovery efficiency of the search, we use a set of 1000 injections. Given the total mass M total = m 1 + m 2 and the chirp mass M,  . This is equivalent to a GW trigger having a merger time uncertainty at Fermi of σ ∆t and a chirp mass uncertainty equal to σ ∆M . In the previous expressions the exponent i stands for the i th trigger. According to the results presented in Figure 6, when used with uncertainties σ ∆t ≤ 0.01 s and σ ∆M ≤ 0.01 M , the CDF chirp targeted search has a relative error with respect to the case of perfect measurements, i.e. (σ ∆t , σ ∆M ) = (0 s, 0 M ), lower than 10% over the range LLR ∈ [0,200]. This means that if we have a measurement with such low uncertainties, and if the output of the chirp targeted search, when used with setting N bins ≤ 10, indicates a trigger with LLR ≥ 100, then the EM candidate is promising, because it is above the LLR background range corresponding to −1σ lower limit.

RESULTS ON LIGO-VIRGO DETECTIONS
We finally apply the search method described previously to a few LIGO-Virgo detections which could plausibly be associated with a modulated γ-ray counterpart: the BNS mergers GW170817 (Abbott et al. 2017) and GW190425 (Abbott et al. 2020a); GW190814 (Abbott et al. 2020b), whose heavier object is a black hole, while the lighter object is either the heaviest neutron star or the lightest black hole observed to date; and the neutron star-black hole mergers GW200105 and GW200115 (Abbott et al. 2021a). Binary black hole mergers, although much more plentiful, are not considered as likely sources of detectable γ-ray counterparts in this work, and are not investigated. Moreover, the search for EM counterparts to such binary systems would impose a large trials factor due to the high rate of these mergers and the large uncertainties on both the chirp mass and the sky localizations (Abbott et al. 2021b). Hence, brighter EM counterparts may be required for a detection.
Observer-frame chirp masses and mission elapsed times (MET) since 2001.0 UTC (decimal), as well as the corresponding uncertainties, for the five GW events, are given in Table 1. Concerning the MET of GW170817 (respectively GW190814), 3 ms (respectively 20 ms) have been subtracted from the geocenter merger time posterior, in order to account for the angle between the direction to NGC 4993 (respectively the directions representing 90% credible regions of the GW190814 skymap) and the Fermi-Earth center baseline. For the other events, whose sky localizations are poorly constrained, we make the conservative simplifying assumption that the sky location is completely unknown. We then simply increase (decrease) the upper (lower) limits on the merger time by 23 ms, which is approximately the light travel time between the Earth center and the Fermi satellite.
To increase our chance of picking merger time and chirp mass values close enough to the true ones, we make use of a grid in the following way: if the MET (respectively the chirp mass) upper and lower limits are MET min and MET max (respectively M min and M max ), we consider 2D grid points (MET i , M j ), with MET i (respectively M j ) ranging from MET min to MET max (respectively from M min to M max ), such that MET i+1 − the upper and the lower limits, there is at least a grid point (MET i0 , M j0 ) such that |MET 0 − MET i0 |≤ 0.01 s and |M 0 − M j0 |≤ 0.01 M . This working method is motivated by the results obtained in the previous section and summarized in Figure 6. We create the grids corresponding to the five GW events and then we run both the generic and chirp targeted searches with the setting defined by these grid points.
The maximum LLR obtained for each run is converted into a FAP in the following way. For a given GW event with merger time t c , the chirp targeted search (and the generic targeted search) background distributions are obtained by running the search on 1000 random offsource times covering the interval [t c − 1 2 month, t c + 1 2 month]\[t c − 60 s, t c + 30 s]. In this way, the estimated noise distribution samples the high-energy EM activity in the few weeks around the time of interest, but we avoid the 30 s of on-source Fermi-GBM data preced-ing and following the merger. This prevents any possible candidate counterpart from contaminating the background. Of course, when constructing the off-source background of the chirp targeted search, mass parameters consistent with the particular on-source event of interest are used. The maximum LLR and FAP associated with each GW event are reported in Table 2. None of the FAP values have upper limit below 0.01, indicating that the LLR values associated with all events are compatible with the off-source background fluctuations. We conclude that there are neither statistically significant excesses of photons, nor significant modulated signals prior to the mergers we considered.

CONCLUSION
In this work we present a method, the chirp targeted search, to detect modulated γ-ray precursors to compact   Table 1. The MET, corresponding to the merger time at Fermi satellite, and the observer-frame chirp mass, for the GW detections explored in this work as possible sources of modulated γ-ray precursors. The values appearing here are the median (50th percentile), the upper limits (90th percentile) and the lower limits (10th percentile). A missing MET (respectively chirp mass) lower/upper limit means that the limit value is away from the median value by less than 1 ms (respectively 10 −3 M ). The parameter estimates are given in Romero-Shaw et al. (2020) for GW170817, Abbott et al. (2019) for GW190425 and GW190814, and Abbott et al. (2021a) for GW200115 and GW200115.
binary mergers in Fermi-GBM data. The existence of such signals is not confirmed so far. If they exist, there are several physical mechanisms which might be responsible for their emission. This fact makes the lightcurve amplitude dependence difficult to predict. Despite these difficulties, the presented method is very general. It aims to look for an excess of photons in the orbital phase space, while the GW frequency evolution is defined by the first order term in the post-Newtonian expansion. The sensitivity of the method has been tested on simplified lightcurves, for which the EM emission takes place only during the same orbital phase window. The performance of the chirp targeted search has been compared to that of an existing, more generic targeted search, which aims to detect an excess of photons in the time space. It has been displayed that the chirp targeted search has higher sensitivity than the generic targeted search, when the signal is modulated by the GW frequency of the binary. Finally, both pipelines have been used to search for EM precursors associated with confident GW events having a non-negligible probability to contain a neutron star, namely GW170817, GW190425, GW190814, GW200105 and GW200115. We found no significant candidate precursor signals associated with any of those events. However, given the potential of new physics provided by the presence of precursor signals, the proposed method, here demonstrated, will be applied to future BNS and NSBH observations by LIGO-Virgo-KAGRA, especially since an increased number of events is expected in the coming years. It is possible that this method could also be important for stellar-mass binary black hole mergers. As an example, these mergers could occur in gas-rich environments, such as active galactic nuclei disks, and could produce EM counterparts with orbital modulations, similarly to what has been proposed for more massive binaries in the LISA band (e.g. Tang et al. 2018). Potential EM counterparts have been reported for both GW150914 (Connaughton et al. 2016) and GW190521 (Graham et al. 2020). As stated above, however, in order to apply this method to the large number of binary black hole events detected by current ground-based interferometers, it will be necessary to address the statistical and computational challenges.
Although a sensitivity gain has been achieved by means of the chirp targeted search, improvements can be envisaged in the future. The actual method is trying to recover signals with spectra described by the related Band function. One can argue that this is not the most optimal choice. The Band functions are very appropriate to recover γ-ray burst prompt emission, at the ori-gin of which most likely the synchrotron radiation and the inverse Compton scattering are at play in producing photons. In the case of γ-ray precursors to compact binary mergers, one can imagine other physical emission mechanisms, like thermal emission. A higher sensitivity might then be obtained by including additional spectral templates in the search. Additional filtering strategies may also turn out to be effective. A cleaning of the search output is synonymous with decreasing the false alarm probability, and as a consequence, the sensitivity of the search gets higher.