{'fraction_non_alphanumeric': 0.03112230618527848, 'fraction_numerical': 0.011456292564604908, 'mean_word_length': 4.602759774200293, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Long-duration transient, gravitational-wave search pipeline

As the sensitivity and observing time of gravitational-wave detectors increase, a more diverse range of signals is expected to be observed from a variety of sources. Especially, long-lived gravitational-wave transients have received interest in the last decade. Because most of long-duration signals are poorly modeled, detection must rely on generic search algorithms, which make few or no assumption on the nature of the signal. However, the computational cost of those searches remains a limiting factor, which leads to sub-optimal sensitivity. Several detection algorithms have been developed to cope with this issue. In this paper, we present a new data analysis pipeline to search for un-modeled long-lived transient gravitational-wave signals with duration between 10 and 1000 s, based on an excess cross-power statistic in a network of detectors. The pipeline implements several new features that are intended to reduce computational cost and increase detection sensitivity for a wide range of signal morphologies. The method is generalized to a network of an arbitrary number of detectors and aims to provide a stable interface for further improvements. Comparisons with a previous implementation of a similar method on simulated and real gravitational-wave data show an overall increase in detection efficiency depending on the signal morphology, and a computing time reduced by at least a factor 10.


I. INTRODUCTION
A new era in astronomy began in September 2015 with the observation of gravitational waves (GWs) from the merger of two stellar mass black holes [1]. Since then, the Advanced LIGO [2] and Advanced Virgo [3] have regularly observed a larger volume of the Universe leading, among major discoveries, to the observation of the merger of two neutron stars [4] in August 2017 associated to gamma ray burst GRB190817A [5] followed up kilonova AT2017gfo in NGC4993 [6]. As of mid of 2021, Advanced LIGO and Advanced Virgo have reported ∼ 50 confirmed mergers of compact objects, black holes and/or neutron stars [7].
Yet many GW sources have not yet been observed: core collapse supernova, isolated neutron stars, magnetars, cosmic strings, and the resulting stochastic background of GWs [8]. The diversity of the GW signal expected from these sources require different detection algorithms. When the GW signal waveform is predicted analytically, matched filter techniques can be used. In practice, this concerns mainly compact objects binary coalescence, cosmic strings signals [9] and GWs from pulsars [10,11].When the GW emission is poorly modeled, detection will rely on unconstrained searches that make few assumptions about the characteristics of the signal. In the last twenty years, several search algorithms have been developed, mostly focusing on GW signals of duration less than a few seconds [12][13][14][15][16]. More recently, transient GW signals of longer duration have received attention, bridging the gap between short-duration transient and continuous emission of GWs, and dedicated search algorithms have been developed [17][18][19][20][21][22][23].
Several astrophysical processes could be at the origin of long-duration transient GWs emission, for example, those related to core collapse supernova, compact object binary mergers and isolated neutron stars [18]. Some of them are associated to the most energetic phenomena observed in the Universe. There is evidence [24,25] that core collapse supernovae and long gamma-ray bursts (GRB) are connected to the death of massive stars where the iron core collapses under its own gravity, forming either a black hole or a highly magnetized neutron star, releasing an incredible amount of energy (10 53 erg) mainly through neutrino emission, while ∼ 1% goes into the kinetic energy of the explosion [26]. Once the collapse is triggered, very powerful non-spherical flows develop in the outer region of the proto-neutron star that are expected to generate GWs energetically bounded to 10 44 − 10 47 ergs [27]. The GW emission will last until the onset of the explosion or until a black hole is formed. The signal is expected to be no longer than 1 − 2 s. In the collapsar model, massive stars collapse to black holes either without an initial supernova explosion or via fallback accretion after a successful but weak explosion [28]; a rotating black hole is formed while the inner layers of the star lacks momentum to eject all the matter. Over a period of minutes to hours, 0.1 − 5 M falls back onto the collapsed remnant, turning it into a black hole and establishing an accretion disk. GWs may be emitted by disk turbulence and disk instabilities that may lead to clumping or disk fragmentation [29,30]. The GW signal expected from accretion disk fragmentation would last O(10 − 100) s with a characteristic strain h ∼ 10 −22 at 100 Hz for a source at 100 Mpc [29]. When the core collapse explosion is successful, a magnetar is formed. Convective currents and dynamical and secular nonaxisymmetric rotational instabilities in the proto-magnetar develop and may emit GWs [31]. In both scenarios, a GRB jet is launched either thanks to magnetohydrodynamical processes and neutrino pair annihilation powered by accretion or by the high Lorentz factor outflow that follows the birth of the proto-magnetar.
When a magnetar is formed, gravitational wave emis-sion from viscosity-driven "spin-flip" instability may last hours to days, with a detection horizon of 3 − 4 Mpc for Advanced LIGO/Advanced Virgo detectors and unmodelled searches [32][33][34]. The merger of two neutron stars will form a hot supermassive neutron star; depending on the component masses, the centrifugal forces induced by differential rotation and the stiffness of the nuclear equation of state may allow it survive for hundreds of milliseconds before collapsing to a black hole or form a massive neutron star [35,36]. It is very likely that the remnant is surrounded by an accretion disk that may endeavor instabilities like in the collapsar scenario. If the newly formed neutron star survives more than a few seconds, it could emit long-lived GW through magnetic field-induced ellipticity [37,38] or r-mode instabilities [39], although the precise amplitude of such signals remains unclear. So far, no post-merger GW signals have been detected for any of the binary neutron star mergers found in LIGO and Virgo data [40][41][42].
Isolated neutron stars are another potential source of long-duration GW signals. Sudden speed-ups of the rotation of pulsars observed in radio and X-ray data are followed by a period of relaxation (weeks long) during which the pulsar slows down. GWs may be emitted during this period but the amplitude is expected to be low as the rotational energy changes remain below 10 43 erg [43][44][45][46]. Seismic phenomena in the crust of magnetars are thought to be at the origin of soft gamma repeaters and anomalous X-ray pulsars. Soft gamma repeaters giant flares are associated to huge emission of electromagnetic energy, up to 10 46 erg, followed by long duration quasi periodic oscillations which may be associated to GW emission over the same time scale [47][48][49]. The recent observation of GRB 200415a, suggesting that magnetar giant flare may be a distinct class of short GRB, with a substantially higher volumetric rate than compact object mergers [50], is re-enforcing the interest for magnetar giant flare events in nearby galaxies.
The diversity of long transient expected GW waveforms has lead to the development of algorithms that do not rely on a signal model. Coherent waveburst [13,51] and X-pipeline [15], used for short-duration searches, have been adapted to search for transients with duration up to a few hundred of seconds, while the STAMP excess cross-power algorithm [18] has been developed to target specifically long and very long transient signals lasting up to several weeks. It has been used to search for long duration GW transients associated to GRBs [52], for postmerger GW signals associated to GW170817 [40,42] and adapted to perform an all-sky/all-time search for long duration GW transient in LIGO and Virgo data [53][54][55].
An enhanced version of the STAMP algorithm is presented in this article. It is a complete rewrite in python of the all-sky/all-time STAMP-AS pipeline that was built using the STAMP algorithm library written in Matlab. As such, it has been optimized to search for GW signals of duration in the range 10 − 10 3 s in a large data set at a lesser computing cost than STAMP. It especially implements a hierarchical strategy, similar to the algorithm proposed in [19] to select the most interesting periods of the data without loosing detection efficiency. This paper is organized as follows. In section II, we present the formalism of the analysis and the methods used to generate candidate events in the framework of a 2 detector search. We describe the implementation of the pipeline and the methods used for background and efficiency estimation in Section III. Section IV summarizes the performance of the pipeline over simulated Gaussian noise and real data from the second LIGO-Virgo observation campaign (O2). Finally, we summarize those results in section V and propose several improvements to increase the pipeline sensitivity in the future.

II. OVERVIEW OF A CROSS-CORRELATION GW
TRANSIENT SEARCH ALGORITHM

A. Definitions and conventions
We are considering a network of GW detectors whose strain data time-series s(t) = n(t) + h(t) is a linear sum of independent detector noise n(t) and the detector's response to a GW strain amplitude given by h(t). The detector noise is itself the sum of random noise and non-Gaussian noise transients, or "glitches". The GW signal is assumed to be described by two polarization modes, h + (t) and h × (t) and originates from a pointlike source whose sky-position is given by the right ascension and declination (α, δ). We defineΩ as the direction to the source andh(f ) the Fourier transform of any h(t) time-series. The detector's response to a GW strain is the linear combination of the two polarisations weighted by the detector antenna factors We consider an interval of duration T of GW strain data that are discrete measurements sampled at f s . In the following, the variable t; refers to the time segment start time.
The STAMP algorithm is an extension of the radiometer method developed to detect point-like sources of stochastic background GWs [56]. To estimate the GW strain power spectrum of a transient signal, excess power is searched in frequency-time maps (f t-maps) formed by cross-correlating the data of two spatially separated gravitational wave detectors I and J. Following [18] an estimator of the GW power in a single f t-pixel is given byŶ where Q IJ (t; f,Ω) = 1 is a filter function that takes into account the arrival time delay of the signal in the two detectors whose distance is given by ∆ x IJ and the pair efficiency (3) which weights the GW strain cross-power according to the alignment of the detectors. To normalize the crosscorrelation, we compute the variance ofŶ for which an estimator iŝ where P I (t; f ) is the noise one-sided auto-power spectrum.
We then define the signal-to-noise ratio SNR(t; f,Ω) for a single pixel SNR(t; f,Ω) depends only on the single-detector whitened statisticỹ and the time delay τ ≡Ω · ∆ x IJ /c of the signal in the two detectors.
In the context of an all-sky search, the source direc-tionΩ, and therefore τ , are unknown. An error in the time delay induces a dephasing in the computation of Y (t; f,Ω) that can cause an underestimation of the SNR of coherent signals. A solution is to span all sky-positionŝ Ω and retain the one that gives the largest SNR. That was the strategy implemented in STAMP-AS used to search for long duration transient GW signals in initial LIGO data [53,57] and advanced LIGO data [54,55]. However, the computational time required to process numerous sky positions was a limitation of the pipeline. Besides, background estimation requires repeating, a large number of times, the same coherent cross-correlation of the data streams for each sky position tested using complete f tmaps while a large fraction of the pixels do not contain relevant information. As a consequence the amount of simulated background was restricted to ∼ 100 years, and the number of sky positions tested was limited to a few. All these sub-optimal features resulted in a loss of sensitivity of ∼ 10 − 20% [53].
The PySTAMPAS pipeline addresses these limitations by implementing the hierarchical approach proposed in [19] which consists of first identifying the most interesting clusters of pixels in single-detector auto-power f t-maps. In a second stage, a coherent detection statistic is computed using only the pixels that have been selected in the first stage. The computationally intensive calculations are carried out only once, allowing rapid background estimation without sacrificing the sensitivity gained by the use of coherence and spanning the whole sky positions. The gain in computational performance has also allowed the introduction of the use of several time-frequency resolutions to gain sensitivity to GW signals that may have time-varying frequency evolution. In the following sections, we describe the different computations that are performed at each stage.
B. Single detector stage

Single detector ft-map
The simplest time-frequency representation of the GW time series s I (t) is a spectrogram obtained using onesided Fourier transform of short segments of duration ∆t. The short segments are first Hann-windowed and overlap by 50% with each other such that the pixels resolution is respectively (∆t/2) × (1/∆t) in time and frequencythe factor 1 2 comes from the 50% overlap between short segments.
The spectrograms are whitened by the one-sided power spectral density of each segment P I (t; f ). Two methods to compute the auto-power have been implemented. The first one, inherited from STAMP, takes the average of |s I (t; f )| 2 over time-neighboring pixels in a similar way to the Welch's method. The other one considers the median over the frequency-neighboring pixels. The pros and cons of the two methods are discussed in section IV A 1. For each time-frequency resolution, f t-maps of the whitened statisticỹ I (t; f ) are built.
The duration ∆t of the Fourier transformed segments is an arbitrary choice that depends of the type of signal searched. Long-duration GW searches generally use Fourier transformed segments of duration 1s which are suited to reconstruct signals lasting ∼ 10 1 − 10 3 s. However when the frequency evolution of the signal is changing with time, parts of it can be better reconstructed using different resolutions. In order to improve signal reconstruction as demonstrated by coherent waveburst [13], we opt for a multi-resolution approach which consists in building several f t-maps of different resolutions and combine them into a single, multi-resolution f t-map.

Clustering
The long-duration GW signal signature in f t-maps appears as a cluster of pixels that a pattern recognition algorithm must be able to reconstruct without assuming a model. A year-long data set is typically used. The unknown morphology assumption leads us to consider a seed-based clustering algorithm. The principle is to group high-energy pixels together by proximity, without imposing any preferred morphology for the cluster. For PySTAMPAS, we have adapted the burstegard algorithm, developed for STAMP [57] to multi-resolution f t-maps.
We consider all pixelsỹ I (t; f ) from every map with individual resolution ∆t i × ∆f i . Pixels for which |ỹ I (t; f )| exceeds a given threshold are kept to form a set of pixels for which we keep the time and frequency of the bottom left corner, ∆t i , ∆f i , andỹ I (t; f ). The clustering algorithm starts with a seed pixel, the first pixel in the list, as the order does not matter. All pixels that are above threshold and within a given distance (in time and frequency) of the seed become part of the same cluster, whatever their resolution. Each new pixel added to the cluster is then becoming the seed pixel and the same process is repeated recursively until no more pixels can be added. The next remaining unclustered pixel becomes the seed of the next cluster and the operations are applied again until all isolated pixels have been clustered. To eliminate clusters composed of only a few pixels, we select clusters that have a user-determined minimal number of pixels. The different parameters of the clustering (pixel threshold, radius and minimal number of pixels per cluster) are free parameters that can be tuned considering that the number of operations is proportional to O(N log(N ))), where N is the number of pixels above threshold. As the GW signal energy is spread over many pixels, the threshold on |ỹ I (t; f )| should not be too selective, and the distance between 2 pixels should not be too strict as well.

C. Coherent analysis
Considering all possible detector pairs, clusters from one detector are cross-correlated with the other detector's pixels. At this stage, the clusters can be composed of pixels of different time-frequency resolution. To be able to cross-correlate pixels of different time-frequency resolution, we define virtual pixels that have resolution min ∆t i × min ∆f i . Each of these pixels is assigned a value that is the largest |ỹ I (t; f )| value of all pixels that overlap the virtual pixel. As ∆t i ×∆f i is constant over all resolution maps, the virtual pixels assigned values have the same weight. The same construction of virtual pixels is performed for the pixels of the other detector's f t-map.
The cross-correlation SNR given by Eq. (5) is then computed considering the virtual pixels. As already mentioned, pixel SNR depends only on the time delay between two detectors τ =Ω · ∆ x IJ /c. In an all-sky search, the direction to the source is not known a priori, and an error on the time delay can cause to underestimate the SNR of coherent signals. A solution is to span the time delay parameter space over all possible values for a given pair. The maximal SNR loss due to an error of d Ω corresponding to dτ is The time delay bin size dτ is determined such that the maximal SNR loss is lower than ∈ [0 1] for the maximal frequency considered in each cluster. dτ = arccos( ) 2πf max (8) with f max being the maximal frequency of all pixels of the cluster, which can be much lower than the maximal frequency of the search, reducing the number of time delays to test. In the most general case, we would need to test N τ time delay values between 0 and ∆x IJ /c by steps of dτ to recover a signal accurately. However, because we are considering a phase factor, a degeneracy appears: = cos(2πf dτ ) = cos(2πf (dτ +1/f )). As a consequence, for a pixel at frequency f , it is sufficient to test time delays in the interval [0, 1/f ] instead of [0, ∆x IJ /c] to get the correct phase factor. In the case of a broadband cluster with pixel frequencies between f min and f max , this interval is the largest for f = f min , so we need to test time delay values between 0 and 1/f min by steps of dτ . Finally, the number of time delays to test to recover a signal with an accuracy is In the end, N τ remains rather small (less than a few hundred), especially compared to the thousands of sky positions that need to be tested using a regular grid in sky coordinates α, cos(δ). Since the computations are done only over the small subset of pixels that constitute the cluster, it is possible to test hundreds of time delays in a reasonable time and therefore limit the loss of SNR to = 0.95 regardless of the signal morphology as shown in section IV A 2. The time delay τ 0 that maximizes the sum of all pixel SNRs provides a detection statistic that reflects the significance of the cluster This detection statistic is used to test the hypothesis of a GW signal or the null hypothesis. However, hierarchical processing methods such as PySTAMPAS applied on real GW data tend to bias the selection of triggers because of the presence of noise outliers in one detector. When combined with noise fluctuation in the second detector, such triggers may have large SNR Γ values despite being incoherent. In order to mitigate this effect we estimate the residual noise energy that is left in one detector's data after subtracting the sum of |ỹ I (t; f )| 2 over all pixels belonging to the cluster. We define the quantity For a coherent GW signal, recovered with the right time delay, this residual energy is expected to be much smaller than both SNR Γ and the auto-power energy E I On the contrary, for a cluster due to a noise outlier in one of the detectors E res I may become large in the second detector. We can then define a second discriminant variable in addition to SNR Γ , Finally, we combine these two variables into a single detection statistic Λ defined by Λ should tend to 1 in presence of a coherent GW signal and take << 1 values in case of noise outliers. For convenience we define such that the detection statistic increases with the significance of the trigger. Note that p Λ is not the only possibility to combine SNR Γ and Σ res . We show in Section IV C that p Λ is robust to loud noise triggers using a sample of real data from GW detectors, but other combinations may be relevant depending on the distribution of background noise.

III. DETAILS OF THE PIPELINE IMPLEMENTATION
In the following sections, we describe the implementation of PySTAMPAS in the case of a 2-detectors network, and we propose a generalization to the case of network of more than 2 detectors. In practice, the pipeline is implemented using Python 3 and relies on the GWpy package [58].

A. Data conditioning
The GW detectors' data streams are first searched individually to reveal clusters of energy which may contain coherent GW signals. Real GW detector data are available as an ensemble of time-series of different lengths. For a given pair of GW detectors, only coincident times are analyzed. This reduces the data set to a list of coincident segments of time. For each of the coincident segments, the data are split into windows of duration T win that overlap by 50%. The duration of the data window is a free parameter that can be adjusted to the typical duration of the GW signal that is being investigated. In this study, we use T win 500 s, as it is done in previous longduration searches [54,55]. STAMP was originally designed to search for signals with duration up to several weeks [18]. Although there is no fundamental limitation to extending PySTAMPAS to longer signals, we limit ourselves to signals in the range 10 − 10 3 s in this paper. Working with very large windows leads to dropping up to T win /2 s of data at the end of each coincident segments, and increases the computing cost of clustering.
The data are first high-pass filtered to suppress energy outside the analysis frequency bandwidth whose lower boundary is adapted to the GW detectors' noise spectrum of each data set. Real GW detector data often contain non-Gaussian, short duration spikes ("glitches") [59,60]. When the magnitude of the glitch is large, an excess of energy is present in the f t-maps and generates single-detector clusters with very large energy (orders of magnitude larger than what a real GW signal would generate). The coherent step is usually not able to eliminate them completely and a better strategy consists in gating the data time-series before computing the f t-maps. PySTAMPAS is mitigating the effect of the loud glitches by applying a Planck window on the h I (t) samples that exceeds a fixed threshold 1 . This threshold is a free parameter that should be tuned for each analysis in order to remove most of the glitches without penalizing signal recovery. After this pre-processing step, f t-maps of y I (t; f ) are built.
As shown in all GW detectors' noise spectra [61,62], real GW data contain many spectral artefacts corresponding to mechanical resonances, power lines and pump or fan-like machines surrounding the detectors [63,64]. Most of these spectral lines are of low amplitude and relatively constant over time, while some have a timevarying frequency. These artifacts can generate false long duration tracks in f t-maps. To attenuate the impact of lines, we set to zero ("frequency notch")ỹ I (t; f ) pixels corresponding to a list of frequencies that are constructed following two steps: 1. For each f t-map built, we compute the median valueȳ(f ) over time of |ỹ(t; f )|. Frequencies for whichȳ(f ) is higher than a given threshold are flagged.
2. If a frequency is flagged in more than a given fraction of the total f t-maps, it is added to the list.
This last condition reduces the risk that a monochromatic GW transient of duration T win is mistaken for an instrumental line and notched. One should note, however, that a very long monochromatic transient signal (on the order of weeks or months) could still be flagged if it is spread over a fraction of the total f t-maps higher than the threshold chosen. If a signal crosses a notched frequency, it may be divided into several parts which will be reconstructed by burstegard as separate clusters, reducing the significance of the signal. To reconnect these parts, we implement the findtrack algorithm [57]. If the minimal distance between the corners of two clusters is smaller than a given radius, these clusters are connected and treated as one single cluster.

B. Coincident search
The coincident search is the proper analysis during which true GW signals are searched in the data. The individual detector's f t-maps are searched for clusters of excess energy following the procedure described in section II B. Two lists of clusters are extracted from a pair of detectors. They are saved along with the f t-maps to be analyzed in the coherent stage following the procedure described in section II C. The pipeline produces a list of coherent triggers that are ranked according to p Λ .

C. Background estimation
In order to assess the significance of triggers found in coincidence, one has to estimate the accidental rate of noise triggers caused by instrumental and environmental effects. Like almost all GW transient search pipelines, to encompass any particular effect in the data and augment the total volume of data, we use the time-slides technique to estimate our background [65]. One data stream is time-shifted with respect to the other one by an amount of time greater than the light travelling time between the detectors. Assuming the number of detectable GW signals is small, this assures that the cross-correlated data does not contain a coherent GW signal. In the meantime, non-Gaussian and non-stationary features of the detectors' noise are preserved. By repeating the analysis for many time-shift values one simulates multiple instances of the noise.
In PySTAMPAS, time-shifts are performed considering data streams split over N win windows that are time ordered on a circle. Data are shifted by a multiple of windows (lags) and for each lag by a multiple of ∆t max the maximal time resolution (mini-lags). For example, considering only lags, at the n-th lag, clusters from detector I that have been extracted in window i are matched with detector J data from window (i+n). With this technique, the maximal number of time-shifts is The total background lifetime simulated T bkg is the number of time-shifts performed times the duration of data available for a pair of detectors. The cumulative background trigger rate gives an estimation of the false-alarm rate (FAR) as a function of the detection statistic which is used to rank the triggers.

D. Sensitivity studies
PySTAMPAS performs sensitivity studies by injecting simulated signals into the data. A simulated signal consists primarily of a waveform which describes the two polarizations modes h + (t) and h × (t) of a GW. Waveforms are stored in files in the form of two time series sampled at f s , as well as metadata (duration, frequency range, physical model, etc). A bank of waveforms with various properties is available to sample the rather large parameter space of long-duration transient GW signals with representative signal morphologies.
To compute the detector's response h I (t) to a given GW signal one has to specify a waveform and the following parameters: • the time of arrival t 0 at the center of the Earth; • the directionΩ to the source; • the inclination and polarization angles (ι, ψ) that characterize the orientation of the source's reference frame with respect to the Earth equatorial frame; • a scaling amplitude factor α to modulate the strength of the signal.
Source frame GW polarizations are then rotated to be expressed in the Earth equatorial frame where a + ≡ 1+cos ι 2 2 and a × ≡ cos ι 2 . The polarizations are then time-shifted by the delay of arrival between the detector's position r I and the center of the Earth and rescaled by the amplitude factor α such that finally, where F + I (t;Ω) and F × I (t;Ω) are the detector's sensitivity to + and × polarizations (expressed in the Earth equatorial frame) of a GW signal coming from direction Ω at time t. The computed response is resampled and interpolated to match with the detector's sampling, and the first and last seconds of the time series are tapered with a Hann window to avoid numerical artifacts when the signal starts or stops abruptly. The signal is injected in the data, which are then analyzed the same way as in a coincident search (restricted to the windows that overlap the injection to gain time). An injection is considered detected if the search produces a trigger within the time and frequency boundaries of the simulated signal, and with a detection statistic p Λ larger than a given threshold.
To estimate the detection sensitivity to a given waveform at a given amplitude, a statistically significant number of injections are performed with random starting time, sky position, polarization angle and cosine of the inclination. Starting times are selected in such a way that they always fall within a coincident data segment. By computing the fraction of recovered injections for different signal amplitudes, it is possible to characterize the detection efficiency as a function of signal's strength, which is usually expressed with the root-sum-squared amplitude h rss given by

E. Generalization to a network of detectors
The search algorithm can be generalized in a straightforward manner to a network of N detectors (I, J, K...), constituting p(N ) = N (N − 1)/2 pairs. For a given timefrequency pixel (t; f ) and sky directionΩ, we define the total coherent SNR as the sum of cross-correlated SNRs from all detector pairs: with SNR IJ (t; f,Ω) the coherent SNR computed from Eq. (5) corresponding to the pair IJ. This allows us to generalize the definitions of E res I and SNR Γ for a cluster of pixels Γ, and finally the definition of Λ remains unchanged: The pipeline's implementation does not fundamentally change with N ≥ 3 detectors. The clustering step is performed independently over each individual detector's f t-maps, following the hierarchical method of [19] . For each cluster, cross-correlation is computed for the p(N ) pairs to compute its ranking statistic p Λ . However, as the degeneracy between sky direction and time delay between detectors is broken for N ≥ 3, it is necessary in this case to test all sky positions by choosing uniformly α and cos(δ) and select that position that maximizes SNR Γ . Therefore, a full-scale study of the pipeline's performances over a network of 3 or 4 detectors will be necessary in the future, considering realistic detector's sensitivity curves.

IV. PERFORMANCES AND COMPARISONS
To test the pipeline and demonstrate its performance, we consider 13 waveforms commonly used in longduration searches [54,55] whose main characteristics are listed in Table I. Most of the waveforms are based on astrophysical models and fall into three categories: eccentric inspiral-merger-ringdown nonspinning compact binary coalescence (ECBC) [66], broadband chirps from innermost stable circular orbit waves around rotating black holes (ISCOchirp) [67,68] and accretion disk instability models (ADI) [30]. We include two ad hoc waveforms to better cover the parameter space; a 250 s long sine Gaussian signal (SG-C) with a decay time of 50 s and a 20 s long band-limited white noise burst (WNB-A). These signals of different morphology cover the time-frequency space with durations within 9-290 s and frequencies in the 10-2048 Hz range. In the following we consider the case of a 2-detector search to compare performance with STAMP-AS. If not stated differently, we are using simulated Gaussian noise following LIGO's best sensitivity during the second observing run (O2) [69] to simulate the data from the two LIGO detectors at Hanford (H1) [70] and Livingston (L1) [71].

A. Signal reconstruction
We investigate the effects of several parameters of the pipeline on the detection capability and the signal reconstruction in order to find a set of parameters that maximize the detection of a wide range of different morphology signals, while keeping the computational costs affordable.

Power spectral density estimation
The accuracy of the noise power spectral density (PSD) estimation is playing a central role to reconstruct GW signals efficiently. Yet, this task is complicated in the case of GW detectors as the noise contains non-Gaussian and non-stationary features such as glitches, spectral lines and slow drifts of the noise amplitude.
Consider a detector's strain time series given by s I (t) = h I (t) + n I (t), where h I (t) is a deterministic GW signal and n I (t) random noise. A good estimator of the onesided PSD of the noise is given by the squared modulus of its Fourier transform, This value is not directly accessible because in case of an unknown GW waveform it is not possible to disentangle a priori signal from noise. One has to rely on the observable |s I (t; f )| 2 , that may contain GW signal. Assuming signal and noise are not correlated,  Therefore, an assumption over the nature of the signal h I (t; f ) must be made in order to build an unbiased estimator of P I (t; f ). PySTAMPAS implements two methods to estimate the PSD that are suited for different signal morphologies.
The first method consists in taking the average of |s I (t; f )| 2 over n t symmetrically chosen neighbouring Fourier transformed segments. The underlying assumptions are that (1) the noise is stationary over the time window considered, and (2) no signal is present in the adjacent pixels. As discussed above, (1) is often wrong because of the presence of short glitches in the data, which are therefore not factored in the PSD and appear as signal. Conversely, (2) is wrong when a monochromatic or quasi-monochromatic signal is present in the data, leading these to be mistakenly included in the PSD. Degraded sensitivity to monochromatic signals is a known weakness of STAMP [53].
To address these issues, we propose to estimate the PSD by taking the median of |s I (t; f )| 2 over n f adjacent frequency bins. The pros and cons of this method are opposite to the first one: short glitches are well taken into account and monochromatic signals are better reconstructed. However, signals whose frequency evolution is rapid tend to be less well reconstructed. In case of noise only, both methods provide similar PSD estimates, except that spectral narrow features are better reconstructed with the method averaging the neighboring time-segments pixels, as shown in Fig. 1. We use the median as it is more robust that the average to extreme values. Because of instrumental lines, it is likely that one of the neighbouring frequency bins has pixels with a very high value of s(t;f), which would spoil the PSD estimation.
The effect of the PSD estimation on the signal reconstruction in PySTAMPAS is illustrated in Fig. 2. Two signals with very different spectral morphologies, a broadband ISCO chirp (ISCOchirp-C) [68] and a monochromatic sine Gaussian (SG-C), are injected in

Gaussian noise.
By taking the median over adjacent frequency bins (hereafter referred as frequencymedian PSD) instead of averaging over neighbouring Fourier transformed segments (time-average PSD), the sine Gaussian signal is better reconstructed, but the fast frequency evoluting part of the ISCOchirp is blurred out. The optimal choice of a PSD estimation method depends on the type of signals targeted and the characteristics of the noise, especially spectral lines and/or non-stationary features. Another way to restore the sensitivity to monochromatic triggers would be to consider a very long (∼ 10 3 s) time period to estimate the PSD in the case of the time adjacent pixels method. However, noise from GW detectors tend to become non-stationary over such time intervals at low frequencies (below ∼ 100 Hz) [63].

Source sky location determination
The number of sky locations tested in the coherent step (which reduces to a single time delay parameter in the case of a two detector network) is currently a limiting factor of all-sky searches, and illustrates the necessary tradeoff between detection sensitivity and computational cost [53]. The hierarchical processing implemented in PyS-TAMPAS allows for scanning many positions at a low cost.
In Section II C, we have seen that N τ , the number of time delays between detectors to be tested, depends on the ratio between the maximal and the minimal frequency of the trigger. Here, we investigate empirically the pipeline sensitivity loss as function of the number of time delays for different waveform families.
Signal waveforms are injected coherently into Gaussian noise, simulating data from LIGO Hanford and LIGO Livingston, from a given sky directionΩ 0 , and are recovered by PySTAMPAS. We vary the number of time delays and keep the maximal SNR Γ obtained which is compared to SNR Γ (Ω 0 ), the SNR value corresponding to the true source positionΩ 0 .
The ratio SNR Γ to SNR Γ (Ω 0 ) as function of the number of time delays between detectors is shown in Fig. 3 for sine Gaussians of different central frequency and for a selection of waveforms of different morphology/durations. We compare the number N τ of delays tested to get = 0.95 to the theoretical prediction from Eq. (9) given in Table II. We see that the optimal value of N τ does not depend on the signal frequency, but mainly on its frequency range f max /f min . Monochromatic sine Gaussians are recovered equally rapidly no matter their frequency, and faster than signals of broader band. The empirical values are overall lower than the theoretical ones. This discrepancy comes from the fact that the clustering algorithm does not always reconstruct the entirety of the waveform, leading to a lower effective value of f max /f min . To optimize the detection efficiency while keeping the number of tested sky positions minimal, we fix such that the maximal SNR loss parameter to 5% and N τ is determined for each cluster following Eq. (9).

Multi-resolution and clustering
The energy of long-duration GW signals is spread over a potentially large number of pixels. This would mean it is necessary to lower the threshold on the individual pixel's energy |ỹ I (t; f )| and rely on the clustering algorithm to group all pixels belonging to the cluster. Clustering a large number N of pixels is computationally expensive since burstegard's time complexity is O (N log N ). However, because of the hierarchical implementation, that step is computed only once per f tmap and is therefore no longer a bottleneck for analyzing long periods of data. Yet, the risk is to include pixels due to noise fluctuations and generate clusters that are only composed of noise pixels. By increasing the minimal number of pixels per cluster, one can control the rate of noise clusters that are generated.
Another way to collect, as best as possible, all the energy in the f t-maps is to process the data with a range fmax/fmin Nτ ADI-A of different time-frequency resolutions that match well all the different GW signal shapes. The choice of timefrequency resolutions depends on the waveform, but we have seen that for the diversity of signals we are targeting, a limited number of time-resolutions is enough to improve the detection efficiency of non-monochromatic GW signals. Using a set of 4 resolutions ranging from 4 s × 0.25 Hz to 0.5 s × 2 Hz, we report an efficiency increase by 5 − 40% for the waveforms tested (at constant FAR), compared to 1 s × 1 Hz pixels. We report the rela- TABLE III. Relative increase on the distance at 50% detection efficiency between a single-resolution approach (pixels of 1 s × 1 Hz) and the multi-resolution approach implemented in PySTAMPAS (4 different resolutions from 0.5 s × 2 Hz to 4 s × 0.25 Hz), all other parameters being equal, for the astrophysical waveforms described in Table IV tive increase in detection efficiency for each astrophysical waveform in Table III. It is not possible to perform a fine optimization of all PySTAMPAS parameters for a generic all-sky/all-time search because of the large parameter space, but we present in the next Sections the pipeline performance for both Gaussian simulated noise and real GW data to detect long duration GW signals using the set of parameters given in Table IV.

B. Test on simulated data
We carry out a study with simulated Gaussian noise to test the pipeline as a whole and evaluate its performance.
First, we generate two sets of 14 days of stationary Gaussian GW noise following LIGO's O2 sensitivity to simulate the data from the two LIGO detectors at Hanford and Livingston. We analyze these data with PySTAMPAS, using parameters given in Table IV. Background triggers are generated following the method described in III C. We perform 128, 000 timeslides, simulating ∼ 4, 900 years of background noise accounting for 34 days of CPU time on a dual-core modern processor. As a comparison, the previous version of STAMP-AS took 95 days of CPU time to perform 1, 000 time-slides over the same data, meaning that PySTAMPAS is faster by at least one order of magnitude. On Fig. 4 showing the cumulative false alarm rate (FAR) as function of p Λ , the blue curves correspond to the distribution of simulated Gaussian noise triggers for the two PSD estimation methods; the shape of the two curves is similar, but the median-frequency PSD method produces ∼ 60 more triggers than the time-average PSD. This has little effect on the pipeline sensitivity as the tail of the p Λ distribution are similar.
For each waveform described in Table I, we estimate the detection efficiency as function of h rss following the method described in section III D. We fix a detection threshold corresponding to a FAR of 1/50 yr −1 and determine the value h 50% rss of h rss for which 50% of the injections are recovered. To provide a comparison, we perform the same search with STAMP-AS over the same simulated Gaussian noise. We use the quantity h 50% rss to estimate the detection efficiency of the search. It is inversely proportional to the typical detection range. In Fig. 5, we show the ratio of h 50% rss between STAMP-AS and PySTAMPAS for each waveform and each PSD estimator.
For a majority of the waveforms tested, PySTAMPAS is more sensitive than STAMP-AS, up to a factor 2, with the exception of the ISCOchirp family for which detection efficiencies are worse by down a factor 0.8 − 1 in the best case with the time-average PSD. For this specific family, the single-detector clustering algorithm reconstructs low amplitude signals poorly because the energy is spread over too many pixels. Down to a certain amplitude, most pixels fall below the clustering threshold and the signal is not reconstructed at all. A finer tuning of burstegard could be done to address this limitation, but this type of signal would certainly be better reconstructed by seedless clustering algorithms. This also illustrates the difficulty of tuning the pipeline to maximize sensitivity to a wide variety of waveforms.
The ad-hoc waveforms illustrate the most extreme cases. Detection efficiency is multiplied by ∼ 6 for the monochromatic sine Gaussian signal (SG-C) when the PSD is computed over adjacent frequency bins as compared to STAMP-AS. On the other hand, the large band white noise burst (WNB-A) is not recovered at all with this method, and recovered almost equally well with the time-average PSD. We note that the sine Gaussian is also better recovered using the time-average PSD. This is due to the fact that we consider a wider time window to com-pute the PSD (32 pixels from each side instead of 8).

C. Tests on real data
Real data from GW detectors have non-Gaussian and non-stationary features that challenge pipelines. To understand the behaviour of PySTAMPAS on real GW noise, we analyze LIGO data from the Advanced LIGO and Advanced Virgo O2 observing run downloaded from the Gravitational Wave Open Science Center [72,73]. The chosen period runs from 2017-08-01 00:00:00 UTC to 2017-08-15 00:00:00 UTC and contains 9.21 days of coincident data from H1 and L1. We keep the pipeline's parameters given in Table IV, but switch on the spectral lines removal algorithm described in Section III A. About 5% of the total frequency bins are flagged as spectral lines and notched for each detector. As in the simulated data study, we consider both PSD estimation methods. Cumulative FAR distributions for O2 data are compared to simulated Gaussian noise FAR distributions in Fig. 4. For both PSD estimation methods, an excess of triggers is present compared to the simulated Gaussian noise distributions, meaning that the FAR of the search for a given value of p Λ is higher than with Gaussian noise.
For the frequency-median PSD, the excess of triggers (∼ 20% more triggers in real data than in the Monte Carlo study with Gaussian noise) consist of long-duration (> 50 s), quasi-monochromatic events that correspond to instrumental lines being punctually excited. These lines are too low amplitude and are not excited regularly enough to be flagged by the spectral lines removal algorithm. However, that excess becomes marginal for large value of p Λ and thus is not affecting the overall pipeline sensitivity for this set of data.
Using the time-average PSD method, the excess of triggers compared to Gaussian noise is much larger, by at least 1.5 orders of magnitude. It is dominated by short glitches with frequencies between 20-100 Hz which have passed the gating procedure. They generate triggers with high p Λ that populate the tail of the distribution. To discriminate those triggers, we implement a veto, Rveto, based on the ratio of incoherent energy between the detectors R = E I /E J , similar to what is done for STAMP-AS in [55]. Fig. 6 shows the cumulative distributions of R for background triggers and for triggers recovered for a GW waveform (ADI-A). Vetoing triggers with R > 4 allows to reduce by a factor 5 the number of triggers but more interestingly, the tail of the distribution of p Λ is drastically reduced to approach the Gaussian noise triggers estimation, while no more than 5% of GW signal triggers are vetoed. In this paper, we are just illustrating that the pipeline behavior changes considerably in presence of non Gaussian and non stationary data. We also show that simple post-processing selection criteria can be easily developed and applied with a relatively small penalty for the overall pipeline sensitivity.
As we have done for the study with simulated Gaus- sian noise, we now estimate the detection efficiency of this search with the two PSD estimators and compare it to results obtained by STAMP-AS during the second Advanced LIGO observing run [55] for a FAR of 1/50 yr −1 .
For the time-average PSD, signals with R > 4 are rejected like is done in the background study. Best results obtained for each waveform among the two PSD methods are presented in Table V and compared to STAMP-AS.
The relative detection efficiency depends on the waveform, but the overall PySTAMPAS pipeline efficiency increase observed with real data is very similar to what was obtained on simulated Gaussian data.

V. CONCLUSION
In this paper, we have presented PySTAMPAS, a new data analysis pipeline designed to search for GW of duration ∼ 10−10 3 s in a network of detectors with minimal FIG. 5. Ratio between the hrss at 50% detection efficiency obtained with STAMP-AS and with PySTAMPAS for a FAR = 1/50 yr −1 for both PSD estimation methods. The white noise burst waveform WNB-A was not recovered at all using the frequency-median PSD. A ratio above 1 means that PySTAM-PAS recovers the signal better than STAMP-AS.
FIG. 6. Distribution of the incoherent energy ratio R obtained for background triggers (in blue) and GW signal triggers from the ADI-A waveform (in orange) using the time-average PSD. Rejecting triggers with R > 4 allows for reducing the excess of large pΛ background triggers, while marginally affecting the pipeline efficiency to recover GW signals.
assumptions on the nature and origin of the signal. The search algorithm relies on a hierarchical method, initially designed for a seedless clustering algorithm [19], where candidate events are first identified in single-detectors f t-maps, and a coherent detection statistic is then computed by cross-correlating data streams from each pair of detector. This method provides a significant gain in computational efficiency compared to the initial implementation of STAMP-AS with seed-based clustering, while still benefiting from the increased sensitivity of coherent searches. This is especially critical for all-sky/all-time searches for which both the data set and the parameter space can be very large.  The reduced computational cost allows us to implement several new features to improve the overall sensitivity of the pipeline. The use of multi-resolution f t-maps enables the better reconstruction of signals with fast frequency evolution. An alternative method to estimate the noise PSD is proposed that is best suited for monochromatic and quasi-monochromatic signals. We also introduce a new detection statistic that compares the coherent SNR of an event to the incoherent auto-power in single detectors in order to discriminate coherent GW signals from loud noise events. Additionally, it is now feasible to scan hundreds of sky positions during the coherence stage, and therefore to reduce the loss of SNR due to an error in the sky position to less than 5%. The combination of these features results in a detection efficiency increased by a factor ∼ 1.5 on average compared to the previous version of STAMP-AS with seed-based clustering for the different waveforms tested, which have durations between 8 − 291 s, frequencies between 10 − 2048 Hz, and various spectral morphologies. We note that the changes in detection efficiency are dependant on the type of waveform, with PySTAMPAS performing slightly less well on waveforms from the ISCO chirp family and better for the remaining waveforms. We plan to improve the tuning of the clustering algorithm to address this issue.
PySTAMPAS is able to perform all-sky or targeted searches over a full observing run and a network of detectors, and provides a basis for further developments. For example, the burstegard algorithm has been used here to identify clusters of excess power pixels, but other detection algorithms could be considered, such as seedless clustering [74] or more complex pattern recognition algorithms. This will be need to be done in order for the pipeline to be fully competitive, as show by the example of the ISCO chirp waveforms family, which are currently slightly less well recovered by PySTAMPAS. We have shown that real GW data search requires to develop specific trigger selection to cope with non Gaussian and non stationary features of GW detectors data, but another possibility of improvement could consist in implementing a better identification and subtraction of non-Gaussian features of the GW detectors noise, as well as better discriminant variables.
Acknowledgments -N.C. is supported by NSF grant PHY-1806990. M.C. is supported by NSF grant PHY-2010970. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gwopenscience.org/), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain [72]. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.