Tito Dal Canton

Alexander H Nitz

Bhooshan Gadre

Gareth S Cabourn Davies

Veronica Villa-Ortega

Thomas Dent

Ian Harry

Liting Xiao

Verónica Villa-Ortega

Real-time Search for Compact Binary Mergers in Advanced LIGO and Virgo's Third Observing Run Using PyCBC Live

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The Advanced LIGO and Advanced Virgo gravitational-wave (GW) observatories conducted their first two observing runs, O1 and O2, between 2015 and 2017 (Aasi et al. 2015;Acernese et al. 2015;Abbott et al. 2018). From these two runs, over a dozen confident observations of binary black hole (BBH) mergers and one binary neutron star (BNS) merger were made (Abbott et al. 2017a(Abbott et al. , 2019a;;Nitz et al. 2019a,b;Venumadhav et al. 2019;Zackay et al. 2021a;Venumadhav et al. 2020;Nitz et al. 2020b). The BNS merger, GW170817, was observed within a minute of the data being recorded, was associated with a gamma ray burst (Abbott et al. 2017b) and was subsequently followed up by a large number of observatories spanning the whole electromagnetic (EM) spectrum (Abbott et al. 2017c). Without the realtime identification and localization of GW170817 as a merging pair of neutron stars, these observations would not have been possible.

The third observing run of Advanced LIGO and Advanced Virgo, O3, began in April 2019 and ended in March 2020. Thanks to hardware improvements in the detectors, the ranges for compact binary mergers were expected to increase by 6-85% with respect to O2, depending on the source type and detector, with the largest improvement expected for Virgo (Abbott et al. 2018). The sensitive volume of the run was then predicted to be 3.3 × 10 6 Mpc 3 yr for BNS mergers, and 3.4 × 10 8 Mpc 3 yr for BBH mergers. With this unprecedented sensitivity, we expected to observe up to ≈ 40 BBH mergers, up to ≈ 10 BNS mergers, and potentially make the first observation of a neutron-star-blackhole (NSBH) merger. Indeed, dozens of candidates from O3 have been uploaded to GraceDB1 and announced on the Gamma-ray Coordinates Network2 (GCN). Four such candidates have been published as notable compact binary mergers so far (Abbott et al. 2020a,b,c,d). Many more confident detections from the first half of O3 have also been presented in the GWTC-2 (Abbott et al. 2021a) and 3-OGC (Nitz et al. 2021) catalogs. O3 was the first observing run in which three observatories operated for the full duration of the run, increasing the observing duty cycle of the network, reducing the uncertainty in the sky location of observed events and therefore maximizing the chance of making a multi-messenger observation (Abbott et al. 2018). Rapid processing of data from all three observatories was therefore a crucial requirement.

PyCBC Live, first introduced by Nitz et al. (2018), is one of several search codes analysing the data in realtime to observe compact binary mergers, other analyses being GstLAL (Messick et al. 2017;Sachdev et al. 2019), MBTA (Adams et al. 2016;Aubin et al. 2021), SPIIR (Hooper et al. 2012;Chu 2017;Chu et al. 2020) and the generic-transient search cWB (Klimenko et al. 2016).

Having multiple analyses provides for redundancy, both in terms of the possibility that one of the analysis fails for any reason, and in terms of the independent methodology that each of these analyses applies to identify compact binary mergers. PyCBC Live is based on the more general PyCBC software package (Nitz et al. 2020a) and uses a precalculated bank of compact binary merger waveform models combined with matched filtering in the frequency domain (Allen et al. 2012;Babak et al. 2013). PyCBC Live has been instrumental in many of the GW observations to date, both in O2 (Abbott et al. 2017d,e,f,a) and O3 (Abbott et al. 2020b,a,d).

In this paper we describe the improvements that have been made to PyCBC Live in preparation for, and during, O3. Specifically, we discuss the new techniques that (i) enabled the simultaneous and symmetric analysis of data from three observatories, and the reliable assessment of the statistical significance of observed candidates, (ii) improvements in the handling of instrumental transients, (iii) an updated technique to detect signals in data from a single detector and (iv) a method to rapidly infer the nature of the compact objects involved in a candidate merger. We then evaluate the effect of these improvements in terms of search sensitivity by simulating compact binary signals in Gaussian data at the design sensitivity of advanced GW detectors, as well as in real O3 data from Advanced LIGO and Advanced Virgo. We also evaluate the accuracy of the source classification method and the effect of these improvements on the latency of the produced candidates using O2 data.

NEW METHODS FOR THE THIRD OBSERVING RUN

This section describes the new methodology that was used or developed in PyCBC Live during O3. Each subsection describes changes to a specific aspect of the analysis, and subsections are ordered so as to follow the data flow through the pipeline as much as possible.

Search space and template bank

Although it is an improvement to the configuration of PyCBC Live rather than the code itself, we begin by describing the search space and template bank adopted during O3, as a reference for future work.

The bank covers the same mass, spin and waveform duration space as that proposed by Dal Canton & Harry (2017) and previously adopted during O2. The same waveform models are also employed: TaylorF2 (Buonanno et al. 2009;Bohé et al. 2013) for BNS and lowmass NSBH templates, and a reduced-order frequencydomain version of SEOBNRv4 (Bohé et al. 2017;Pürrer 2014) for BBH and heavy NSBH templates.

However, the O3 bank utilises a template placement method based on an optimized hybrid geometric-random approach (Roy et al. 2017[START_REF] Roy | [END_REF][START_REF] Roy | 42nd COSPAR Scientific Assembly[END_REF] which is more efficient than the "manual" combination of geometric and stochastic methods used for the O2 bank (Dal Canton & Harry 2017) in the sense that the obtained bank is ∼ 25% smaller and can be produced much faster. As a result, the O3 bank is only 13% larger than the O2 bank, despite the increased sensitivity of the detectors.

Improved rejection of instrumental transients

Loud instrumental transients in GW data (glitches) which last much less than a second can corrupt the results of transient searches on a time scale much longer than the glitch itself (Dal Canton et al. 2014). In particular, due to the relatively long-lasting impulse response of the various filters involved in the SNR calculation, loud glitches can cause the SNR time series for a given template to cross the trigger-generation threshold many times over several seconds. The resulting high-SNR triggers then dominate over quieter triggers from the underlying stationary noise and possible astrophysical signals, effectively blinding the search for the entire duration of the impulse response of the filters. In early O3, this issue appeared prominently in the results of PyCBC Live as occasional gaps in the production of triggers from a given detector, lasting from a few seconds to several tens of seconds, depending on the glitch.

A simple and widely-used solution to this problem, called gating, consists of windowing out the GW strain data for a short window centered on the glitch prior to matched filtering (Abbott et al. 2016;Usman et al. 2016;Sachdev et al. 2019). PyCBC offline searches already implement this method by detecting glitches as loud excursions in the whitened strain, and then multiplying the data with the complement of a Tukey window centered on the glitch time. We adopted the same algorithm in PyCBC Live during O3. We used a threshold of 50σ on the absolute value of the whitened strain time series as a glitch detector. Each detected glitch was gated with a symmetric complemented Tukey window, configured to have 0.5 s of central zeroes and 0.25 s of smooth taper on both sides. This approach significantly reduced the impact of high-SNR non-Gaussian transient noise with no visible impact on the latency of the analysis. The chosen gating duration is justified because many high-SNR glitches tend to be shorter than 1 s (Davis et al. 2021) and significantly longer gates might cause more damage to downstream stages of the analysis than the glitch itself. A fixed duration also removes the need to estimate the duration of the glitch, which is nontrivial due to the impulse response of the whitening filter, often longer than the glitch itself. Nevertheless, a fast gating procedure which more accurately identifies the time-frequency extent of a glitch could be beneficial in the future.

Another improvement inherited from PyCBC's offline search during O3 was the inclusion of the highfrequency sine-Gaussian discriminator in the ranking of single-detector triggers. The discriminator, described in Nitz (2018), exploits the fact that most compact binary mergers induce a negligible amount of signal power at frequencies higher than the ringdown of the dominant quadrupole mode. By measuring the excess power at the time of peak signal amplitude, and at various frequencies above the ringdown, a χ 2 statistic is constructed and used to reweight the single-detector trigger ranking statistic. The discriminator is most effective at preventing glitches from triggering high-mass templates with final frequencies of ≈ 100 Hz or less, hence it increases the search sensitivity to high-mass black hole mergers. We adopted the same implementation of the discriminator used by PyCBC's offline search with a negligible impact on PyCBC Live's latency.

Inclusion of Virgo in the coincident search

Advanced Virgo began operating in conjunction with the LIGO observatories in the last few months of O2, when its sensitivity was typically 1/4 to 1/3 that of the LIGO instruments. Despite the smaller sensitivity, the inclusion of Virgo markedly improved the localization of important candidates, such as GW170817 (Abbott et al. 2017a). However, Virgo's contribution to the overall network sensitivity was limited, as quantified by its integrated BNS observed time-volume of 4×103 Mpc 3 yr compared to 5×10 5 Mpc 3 yr for LIGO-Hanford 3 . Hence, in order to produce candidate events, the PyCBC Live analysis introduced in Nitz et al. (2018) relied on a simple coincident detection between the two LIGO observatories only. Additional detectors were analyzed by Py-CBC Live for the purpose of improving the rapid spatial localization, but they did not contribute to the significance of candidates, and they could not produce candidates in coincidence with one of the LIGO detectors.

However, as the relative sensitivities of the instruments within the global GW network become comparable, each instrument's contribution to the overall sensitivity also increases. In the coming years, the Virgo observatory may approach 60-80% of the LIGO detector sensitivities, limited primarily by its shorter arm length (Abbott et al. 2018). In addition, new instruments will be joining the global network: KAGRA (Somiya 2012;Aso et al. 2013;Akutsu et al. 2020), which conducted its first observing period shortly after the end of O3, and LIGO-India (Iyer et al. 2011), scheduled to begin operation in the mid 2020's (Abbott et al. 2018). The higher network sensitivity comes about from two sources: (1) improvement in overall network uptime due to overlap between instruments live time, and (2) improvement in detection confidence (reduction of the false-alarm rate) from additional detectors. Hence, in order to start exploiting the benefits of a larger and more symmetric network, the PyCBC Live analysis has been modified for the O3 run.

In its O3 configuration, PyCBC Live correlates the full bank of template waveforms with data from all operating detectors. Then for each detector pair, we independently perform the same double-coincident analysis used for LIGO-Hanford and LIGO-Livingston during O2.4 When a pair of detectors labelled A and B observe a coincident candidate, a false-alarm rate F AB is computed using time-shifted analysis of these two detectors's data, as done in O2. The F AB estimate is local, using only the last 5 hr of data, and therefore tracks possible slow changes in the properties of the detector noise. At times when A and B are the only operating detectors, any candidate events are then submitted directly to GraceDB with false-alarm rate F AB , the analysis being effectively identical to O2. However, if additional detectors are observing at the time of the candidate (C, D. . . ), a combined false-alarm rate is computed as follows.

First, we correct the double-coincident false-alarm rate to include the trials factor associated with the choice between possible detector pairs,

F tf AB = N 2 F AB (1)
where N is the number of detectors that can generate double coincidences at the time of the candidate.

If multiple instrument pairs generate coincident candidates from the same transient, we select the candidate having the lowest false-alarm rate. If tied, we select the candidate with the largest ranking statistic value. Then, using the template of the selected candidate, we calculate the signal-to-noise ratio (SNR) time series for the remaining detectors, which we call followup detectors, and then use this time series to obtain a p-value for each followup detector. Assuming a plane GW traveling at the speed of light, and using the arrival times estimated at detectors A and B, an on-source time interval of possible arrival times is determined at each followup detector. The maximum SNR within the on-source interval is identified. Next, 150 s of data immediately before the on-source window, called off-source data, are segmented into N off intervals of the same duration as the on-source data (Figure 1). The maximum SNR in each off-source interval is calculated, and the number M off of off-source intervals having a SNR larger than the on-source SNR is used to compute the p-value,

p C = 1 + M off,C 1 + N off , (2)
where C labels the followup detector(s). This is the probability of producing a SNR as large as the one observed under the assumption that detector C's data contains only noise. Such a p-value is produced for each followup detector (p C , p D . . . ). In addition, we obtain the p-value for the original double-coincident candidate as

p AB = 1 -exp -T live F tf AB (3)
where T live = 4.38 hr (0.005 yr) is a fiducial livetime used to convert false-alarm rates to p-values. Its value is close to the amount of single-detector data stored for background computation, and also close to the minimum inverse false-alarm rate required for uploading a candidate to GraceDB. The two-detector p-value p AB is then combined with the followup detector p-value p C using an adaptation of Fisher's method [START_REF] Fisher | Statistical methods for research workers[END_REF]:

p ABC = P (2, -ln(p AB p C )) = (4) = p AB p C [1 -ln(p AB p C )] (5)
where P(•, •) is the regularized gamma function. (For more than one followup detector, our algorithm performs this combination iteratively.) Next, we convert back to a false-alarm rate,

F ABC = -T -1 live ln(1 -p ABC ). (6)
For the final significance, we additionally select the minimum of the original two-detector and multi-detector false-alarm rates,

F = 2 min{F tf AB , F ABC }, (7)
where the trials factor of 2 accounts for this additional choice. This procedure produces a self-consistent rate of false alarms under the null hypothesis, and assuming time-invariance of the noise for the collected background. The latter assumption may be violated if one or more detectors rapidly change to a new operating state, hence the desire is to use as short a background collection time as possible to track these changes.

We now apply the above equations to two limiting examples to illustrate possible results in practice. We first consider the case of a nearby source and 3 detectors with similar sensitivities. The source produces a very loud coincidence in detectors A and B as well as a very high SNR in detector C, where by "very loud" we mean that the SNR is higher than any background it is compared to. In this case we obtain the lower bound (F AB ) -1 > 100 yr, determined by the amount of data chosen to estimate the double-coincident background. Assuming an on-source window of 40 ms for detector C, we also have the upper bound p C 2.7 × 10 -4 . Then, using the above formalism, we obtain F -1 3500 yr. Consider now a similar situation, but the sensitivity of detector C is so low that the signal is completely undetectable in its data, leading to p C ≈ 0.5. This results in the bound F -1 17 yr, weaker than the original bound on F AB due to the combined effect of the trials factors and a detector with low sensitivity. Nevertheless, the limit is still well below the threshold required for a public alert and for considering the candidate worthwile of additional astronomical observations. In these examples, we emphasized that the resulting false-alarm rates are to be understood as upper limits, as they are limited by the amound of data chosen to estimate the background, and not by the SNR of the candidate. The limits can be lowered by using more data, but this comes at the risk of increasing our sensitivity to sudden changes in the noise properties of the detectors.

Note that the on-source SNR in the followup detectors is not required to cross any threshold. Hence, this method naturally allows even weak signals in the followup detectors to increase the significance of the original coincident candidates. However, a potential limitation of the method is represented by the implicit assumption that all detector pairs are equally likely to The window of possible signal arrival times at Virgo (on-source region, dark gray in the bottom-right panel) is calculated using the Hanford and Livingston triggers, as indicated by the horizontal arrows. The Virgo SNR time series is calculated (violet, thick curve) and searched for its maximum within the on-source region. Once the maximum is identified, its statistical significance is determined by comparing it to the maxima occurring in the off-source intervals (vertical stripes of alternating color in the top and bottom-left panels).

produce an initial double coincidence in response to a signal. This leads to the trials factor in Eq. 1 being potentially very large: for a network of 3 detectors, such as during O3, the total trials factor can be as high as 6. As additional observatories join the network, the trials factor will grow rapidly and this method may not offer a sensitivity close to optimal. One way to reduce the trials factor would be to use the local sensitivity of each instrument to predict the instrument pair(s) most likely to produce a detection, and only consider those pair(s) in the calculation of the combined false-alarm rate. An alternative strategy for a heterogeneous network could be modifying Eq. 1 to weight each detector pair in a way that accounts for different sensitivities. Finally, the iterative application of Fisher's method de-scribed above could be replaced by a single application of the method to all available p-values. Although the two approaches produce p-values well within a factor of 2 in the case of 2 followup detectors (i.e. a LIGO-Virgo-KAGRA network), a single combination may lead to a higher sensitivity in the case of a 5-detector network.

Identification of signals in a single detector

During an observation run, there will undoubtedly be times during which a detector cannot operate, or is affected by severe nonstationary or transient noise. Moreover, all detectors are blind to signals coming from certain positions in the sky, and these positions are detector-dependent. Therefore, a non-negligible fraction of signals are unable to produce coincident trig-gers in two or more detectors (Callister et al. 2017;Nitz et al. 2020b). By chance, two events that were very promising targets for EM followup observation were first identified as single-detector triggers: GW170817, initially seen as a single-detector trigger in LIGO Hanford due to a large glitch affecting LIGO Livingston (Abbott et al. 2017a) and GW190425, observed by LIGO Livingston and Virgo, but too weak to be detectable in Virgo (Abbott et al. 2020a). In cases like these, we have a single-detector trigger and no coincidence. We cannot rely on the usual robust time-slide approach to establish the false-alarm rate of the trigger, and an alternate method is required. We note and caution that formally, the false alarm rate can only be assigned with confidence to be less than once per livetime for singledetector candidates. Beyond this point, extrapolation is used in order to generate low-latency alerts. If the detector noise changes in unexpected ways, the extrapolation may become invalid and the rate of false positives for single-detector candidates may no longer match an extrapolated false alarm rate.

During O2, PyCBC Live relied on a simple algorithm to identify single-detector candidates based on a pre-determined ranking statistic threshold, and a set of signal-consistency cuts, which were chosen based on early detector data. The algorithm was only active when a single detector was observing. We improve on this method by implementing a complete ranking statistic and procedure for extrapolating the false-alarm rate. We further increase the coverage of our single-detector search by allowing it to operate during times when multiple detectors are observing, as the relative sensitivities of two detectors may imply that a signal can only be seen in one of them.

Our single-detector trigger ranking statistic is the usual reweighted SNR (Usman et al. 2016)

, ρ = ρ (1 + (χ 2 r ) p/2 )/2 -1/p if χ 2 r > 1 ρ if χ 2 r ≤ 1 , ( 8
)
where ρ is the matched-filter SNR, χ 2 r is the timefrequency discriminator described by Allen (2005), and p is an index set to the usual value of 6. Note that the same ρ, calculated for each detector and complemented by other terms, also defines the ranking statistic of coincident triggers.

As we do not have a coincident trigger to corroborate the evidence of a signal, we must apply strict cuts in order to remove triggers which are likely to be glitches. Therefore we remove any trigger with ρ < 9 or χ 2 r > 2, which comes at the cost of possibly removing some signals that are not particularly well matched by the templates. We further restrict the calculation to triggers coming from a template describing a system with a nonnegligible mass remaining outside of the final black hole. By doing so, we ignore many of the templates which match best to common types of glitch in the LIGO detectors, as well as focus on the templates corresponding to sources which are most interesting for potential generation of EM counterpart emission. Applying the remnant mass equation from [START_REF] Foucart | [END_REF] to all templates in our bank, we find that this mass is negligible for templates with duration shorter than ∼ 7 s. Therefore, we prevent triggers associated with shorter templates from generating single-detector candidates.

The probability distribution of ρ for triggers associated with detector noise is expected to follow a falling exponential, as described in Davies et al. (2020),

p(ρ|noise) = α exp [-α(ρ -ρth )] if ρ > ρth 0 if ρ ≤ ρth , (9)
where ρth is a threshold on the reweighted SNR, and α is the fit coefficient. Given the selection cuts defined above, we find empirically that this model describes the detector noise very well.

For each detector, we fit Eq. 9 to the triggers from a day's worth of data and record the fit coefficient, combining these over a longer period of time. Then, during PyCBC Live's operation, we can use this coefficient to estimate the false-alarm rate for each single-detector trigger that passes the selection cuts described above. In general, this fit could be performed separately for each template, but we instead choose to perform the fitting in five bins which are spaced logarithmically in template duration. This choice allows us to group many templates which behave similarly in the presence of noise, and hence increase the number of triggers available for the fit.

Noise in ground-based GW detectors is affected by slowly-varying environmental processes, like weather, and by commissioning activities that change the detector properties over time. Therefore, the statistical properties of the noise are time-dependent, and the fit coefficients for each bin will also vary over time. To account for this variation, we combine the daily fit values in one of two ways. The maximum likelihood choice of α is proportional to the inverse of the mean ρ for each template (Davies et al. 2020), and so we can take the mean of α -1 over the different days, weighted by the number of triggers from each day in order to combine these fits. Alternatively, we could take a low percentile of the α distribution over different days. This would lead to a much more conservative estimate of the false alarm rate. In our case we use the 5 th percentile value, and call this the conservative fit coefficient. This conservative choice would be used for issuing alerts during future observing runs, as it would reduce the number of false alarms compared to the use of the mean coefficient.

We see the variation and combination of the trigger fit coefficient in Figure 2, where the fit coefficient in each template duration bin is plotted for each day of July 2017, along with the mean and conservative combinations.

To calculate the expected rate of louder triggers, as well as the trigger distribution of Eq. 9, we must estimate the overall rate of triggers. This is done by simply counting the triggers which pass the specified thresholds in the daily fits, and for the mean fit, this is simply an addition of the daily trigger counts. For the conservative fit choice, we choose the 95 th percentile daily trigger count and multiply by the number of days over which the fit smoothing is performed.

Using the originally recorded trigger from LIGO Hanford, and the conservative fit values calculated from the O2 data from July 2017, we estimate the single-detector trigger false alarm rate in LIGO Hanford of GW170817 to be 1 per O(10 9 ) yr. If we choose the mean fit combination, then the estimated false alarm rate would be 1 per O(10 11 ) yr.

Source classification

Source classification between different possible types of coalescing compact binary is an important element in generating followup alerts for EM or other counterpart signals. For the O3 run, four astrophysical source types: BNS, BBH, NSBH and MassGap, were considered (LIGO Scientific and Virgo Collaborations 2019). The desired classification designated every object with a mass below 3 M as a neutron star, every object with mass above 5 M as a black hole, and every object between these two thresholds as a MassGap object; any binary containing one or more MassGap components was then considered as MassGap. Accurate classification in low latency is a considerable challenge: in general for lower-mass binaries, only the chirp mass

M = (m 1 m 2 ) 3/5 (m 1 + m 2 ) 1/5 (10)
can be precisely measured, while the mass ratio has large uncertainty. In addition, typically search pipelines report only a point estimate of redshifted component masses, and these template masses may be subject to systematic biases relative to the true source-frame component masses (see e.g. Biscoveanu et al. 2019). During O3, astrophysical source classification for PyCBC Live candidates was performed by the LIGO/Virgo rapid alert infrastructure using a "hard cuts" method, which Figure 2. Fit coefficients α calculated daily for triggers from a month of O2 LIGO Hanford data which meet the ρ, χ 2 r and template duration cuts as described in the text. One plot is given for each template duration bin, and longer templates generally have fewer triggers at high ρ, so the fit coefficients are larger. The dashed and dotted lines show the mean and conservative combinations of the fit coefficients, used to estimate the false-alarm rates of future single-detector candidates, and the black lines are for comparison to the fit values if they were not separated into different duration bins.

assigns Boolean weights (either 1 or 0) to the different source types based on component mass cuts applied to the reported search template. The fact that this method neglects uncertainties in component masses and does not account for any systematic error in mass recovery suggests the potential for improvement (compare Kapadia et al. 2020).

A new approach developed during the later part of the O3 run, which will be described in detail in Villa-Ortega et al. (2022, in preparation), uses the chirp mass recovered by the search pipeline as input, and implicitly assumes a uniform density of candidate signals over the plane of component masses m 1 , m 2 . These sourceframe masses are constrained to the interval 1 M < m < 45 M , where the lower bound is the lower limit on the template bank mass space described in the Section 2.1, and the upper bound is chosen based on BBH population studies up to the first half of the third observing run, O3a (Abbott et al. 2019b(Abbott et al. , 2021b)). Constraining the chirp mass to be within a confidence interval around a point estimate derived from the search template determines an allowed region in the m 1 -m 2 plane; we then estimate the probability that a candidate belongs to each source category to be directly proportional to the area of the allowed region satisfying the criteria for the given category. 5 The output of the method for a given event is a list of probabilities {P BNS , P NSBH , P MG , P BBH } summing to unity.

To compute these areas we require accurate estimates of the candidate chirp mass in the source frame, thus we apply a correction for the bias caused by cosmological redshift. For this correction we use estimates of the luminosity distance and its uncertainty derived from the effective distances (Allen et al. 2012) of the trigger(s) comprising the event:

DL = const • min(D eff ), ( 11
) σDL = DL • const • ρ -p coinc , (12)
where min(D eff ) is the minimum effective distance over all triggers obtained for a given coincident or singledetector event, and the constants of proportionality and power law exponent p may be derived from a fit to previously obtained three-dimensional localizations produced by the BAYESTAR pipeline (Singer & Price 2016) for PyCBC Live events. Although in most cases the uncertainty in the source chirp mass is expected to be dominated by the redshift (distance) uncertainty, we combine this with a nominal, small uncertainty of 1% in the detector-frame chirp mass relative to the value recovered by the search pipeline.

Our approach may be compared to the 'ellipsoidbased' method outlined in Chatterjee et al. (2020). The ellipsoidal error region considered there accounts for expected uncertainties in the recovered masses and spins in the limit of Gaussian detector noise and high SNR (though without attempting to correct for the source redshift). It was also noted there that the actual recovery of parameters other than M by search pipelines was subject to significantly larger systematic error, motivating an alternative machine learning based method. Our approach is simpler in that we effectively consider the uncertainties in such parameters to be infinitely large: we leave improvements to this approximation to future work.

Search for maximum-SNR template

The rapid spatial localization of GW alerts performed by BAYESTAR (Singer & Price 2016) relies on the knowledge of the template that maximizes the likelihood assuming stationary Gaussian noise (or equivalently, the template that maximizes the network matched-filter SNR) for a given candidate. The template is usually taken directly from the candidate produced by the lowlatency search. However, when a search like PyCBC Live generates a candidate in response to an astrophysical signal, it includes both astrophysical priors and the presence of nonstationary and non-Gaussian noise features into the significance of the candidate. Hence, the template immediately associated with a candidate will not necessarily maximize the network SNR under the assumption of stationary Gaussian noise. The sparseness of the template bank will generally also drive the reported template parameters away from the maximum SNR. This could in principle introduce biases in the rapid spatial localization, and more generally affect any low-latency result which uses the mass and/or spin parameters of the search template, such as the source classification we described above in Section 2.5. See Biscoveanu et al. (2019); Chatterjee et al. (2020) for investigations of such possible biases.

In order to remove some of the sources of these biases, we developed a followup process that starts after PyCBC Live reports a candidate to GraceDB. The process reanalyzes a short amount of strain data around the candidate, and uses differential evolution (Storn & Price 1997) to find the template parameters that maximize the network SNR. The maximization explores the mass and spin parameter space in a continuous fashion, regardless of the placement of the search templates. Once the optimization converges, or a predefined timeout of 400 s is reached (whichever comes first) a new candidate is uploaded to GraceDB using the best tem-plate found by the optimization. The new candidate can then be used to generate new spatial localization and source classification results, free of potential biases from the initially-reported template.

EVALUATING THE IMPROVED SEARCH TECHNIQUE

In this section we evaluate the impact and performance of the techniques described in Section 2 using both simulated and real data.

Sensitivity in simulated data

We first characterize the search sensitivity by simulating a population of astrophysical signals, adding the signals to a portion of simulated Gaussian noise, analyzing the data with PyCBC Live and counting how many signals are recovered at a given false-alarm rate. The noise models correspond to final design sensitivities of Advanced LIGO and Advanced Virgo (Abbott et al. 2018). We focus on evaluating the significance calculations described in Sections 2.3 and 2.4. To this end, we compare the sensitivity of the search under different network configurations: HL, HLV, H, L and V, where H, L and V indicate respectively LIGO-Hanford, LIGO-Livingston and Virgo. In the HL and HLV configurations, all observatories are assumed to be observing at the same time.

We construct a population of BNS signals with component masses distributed uniformly between 1.35 M and 1.45 M . Spins are assumed to be aligned with the orbital angular momenta, and spin magnitudes are distributed uniformly between 0 and 0.05. The signals are simulated using a waveform model based on post-Newtonian theory. The sources have isotropic orientation and sky location. In order to increase the number of detected signals, we distribute the sources uniformly in chirp distance (Allen et al. 2012) up to a maximum value. When computing the sensitive volume, we then weight each source such that the effective population has a uniform spatial distribution, as described in Usman et al. (2016).

The result of the simulation is shown in Figure 3 in terms of sensitivity distance as well as relative search volume between the HLV and HL configurations. We can see that adding Virgo to the LIGO network increases the detection rate of BNS systems by a few tens of percent under ideal noise conditions. The single-detector distances are approximately half of what is achieved by a multidetector network.

Sensitivity in real data

Here we repeat a similar test as presented in Subsection 3.1, with the difference that we consider broader ranges of masses and spins (therefore effectively including BBH and NSBH systems) and we add their signals to real data from the third observing run of Advanced LIGO and Advanced Virgo, as opposed to simulated stationary Gaussian detector noise. We use ∼ 8 days of O3 data starting from 2019-05-04 13:15:32 UTC. In the simulated binaries, neutron stars have masses distributed uniformly between 1 M and 3 M , and spin magnitudes distributed uniformly between 0 and 0.05. Black holes have masses distributed uniformly between 3 M and 100 M , and spin magnitudes distributed uni- formly between 0 and 0.985. Spins are aligned with the orbital angular momenta in all cases. BNS signals are simulated using a post-Newtonian waveform model, while NSBH and BBH signals use the SEOBNRv4 opt inspiral-merger-ringdown model (Devine et al. 2016;Bohé et al. 2017).

The resulting sensitivity for BNS mergers is shown in Figure 4. We can see that the improvement in detection rate is around 10% at false-alarm rate thresholds relevant for public alerts. This estimate is consistent with the earlier estimate from simulated data at design sensitivities. NSBH and BBH mergers, albeit arguably less interesting for rapid alerts, show similar relative improvements at the same false-alarm rate threshold, and are not shown here.

Redetection of GW170814 and GW170817

PyCBC Live in the O3 configuration detects GW170814 in O2 data as a coincidence between the two LIGO detectors, with a LIGO-Virgo network SNR of ≈ 15. With the method described in Section 2.3, and using 5 hours of lookback background, the candidate is assigned a false-alarm rate of 1 in ≈ 200 yr.

The case of GW170817 is more interesting because of the glitch affecting the LIGO-Livingston data seconds before merger (Abbott et al. 2017a). Although the glitch is automatically gated by PyCBC Live, the surrounding data are still flagged as affected by a glitch by the low-latency data-quality flags, preventing a LIGO double coincidence from taking place. In addition, the small Virgo SNR also prevents a double coincidence between LIGO-Hanford and Virgo. Nevertheless, using the method described in Section 2.4, GW170817 is reported as a LIGO-Hanford single-detector candidate with falsealarm rate of 1 in ≈ 10 9 yr.

PyCBC Live can also be configured to ignore dataquality flags. Under this configuration, GW170817 is detected instead as a LIGO double coincidence and is assigned a false-alarm rate lower than 1 in ≈ 17 yr by the method of Section 2.3. The much higher value with respect to the single-detector candidate is due to the upper limit on the double-coincident false-alarm rate imposed by the duration of the lookback background (1 in 100 yr), combined with a trials factor of 6 caused by having 3 observing detectors at the time of the candidate. The absence of a detectable signal in Virgo produces a relatively large followup p-value (see Eq. 2), which cannot overcome the penalty of the trials factor. In fact, this situation matches our second numerical example in Section 2.3. Hence, when comparing these GW170817 false-alarm rates, one has to bear in mind that the LIGO-Hanford-only rate is an extrapolation from months of data, while the Hanford-Livingston-Virgo rate is an upper limit based on just 5 hr of data.

In both configurations, however, GW170817 is reported with a false-alarm rate well beyond what is required to issue a public alert on the GCN and to consider the candidate worthwhile of followup observations.

Latency

We measure the latency of the analysis by repeatedly replaying a week of O2 data and analyzing it with Py-CBC Live, thus simulating an actual observing run with a realistic computing configuration. The test amounts to a total wall-clock time of approximately 50 days. For each candidate uploaded to GraceDB, we calculate the latency as the difference between the upload time recorded by GraceDB, and the merger time estimated . Cumulative distribution of PyCBC Live's latency, for a period of replayed O2 data analyzed using the O3 configuration. The latency is defined here as the time elapsed between the estimated coalescence time of a candidate and the time of creation of the corresponding GraceDB event. The shaded region is the range of expected latency due to PyCBC Live alone, as described in Nitz et al. (2018).

by PyCBC Live. This quantity includes the processing time in PyCBC Live, as well as the latency due to the generation and distribution of the replayed O2 data. We expect the former to be the dominant contribution.

The cumulative latency distribution is shown in Figure 5. Most candidates are available in GraceDB within a few tens of seconds after coalescence, as expected. The tail extending to ≈ 100 s is due to occasional and temporary issues with the computing infrastructure running the test, typically starvation of the available computational resources or interruptions of network connections. A similar tail is also found in a typical production analysis.

Astrophysical classification

In order to test the chirp-mass-based classification of Section 2.5 we used the same population of simulated signals described in Section 3.2, with an additional constraint. Given the 1 M < m < 45 M limits imposed on component masses, the results for asymmetric highmass NSBH systems outside these limits are not representative of its accuracy of the method: here, we restrict the black hole components of simulated NSBH events to be below 50 M . For each simulation recovered by the search, we find the estimated source probabilities {P BNS , P NSBH , P MG , P BBH }. We then consider the figures of merit shown in Figure 6.

The first figure of merit is the distribution of probabilities for correct classification, plotted for each category as determined by the true source masses. The great ma-jority of BNS and BBH simulations are assigned high or very high correct class probabilities, as expected given the positions of the target classes in the m 1 -m 2 plane. No NSBH simulations are assigned very high probabilities of NSBH origin, since their contours of constant chirp mass always overlap other source target regions to some extent, but the majority are assigned P NSBH > 0.5. In contrast, the majority of MassGap simulations are assigned P MG 0.5; this again is expected due to the narrow extent of the MassGap region and its very high overlap with other source target regions.

The second figure of merit is the correctness of the highest estimated probability for each simulation, i.e. the most likely source category as assigned by our method. Comparing this highest probability to the true (target) classification determined by the simulated source masses, we construct the confusion matrix shown in the right panel of Figure 6. Simulations of all source types except MassGap are assigned most likely classifications that are correct in a large majority of cases; Mass-Gap simulations are, though, preferentially assigned as most likely to be NSBH. Given the very high uncertainties on the rates and masses of actual NSBH and Mass-Gap sources (e.g. Abbott et al. 2020c), this bias can be argued to be acceptable, and will also yield a conservative outcome as the method will err on the side of recommending EM followup for signals consistent with NSBH origin even if the true source class is MassGap.

DISCUSSION

We described how the PyCBC Live analysis was improved with respect to the O2 configuration between the end of O2 and the end of O3. The most significant changes are the inclusion of more than two detectors in the significance calculation, which allowed Virgo to play a prominent role in the generation of O3 candidates, and the ability to assign significances based on data from a single detector. The single-detector significance calculation and source classification methods were not used during O3 due to its premature end, but they are ready to be utilised in future runs. We evaluated these improvements in multiple ways: first by recovering simulated signals added into ideal detector noise, then by recovering simulated signals added into real O3 data, and finally by reanalyzing a segment of O2 data containing the GW170814 and GW170817 events, and discussing how these events are detected by the new analysis. The improvements we introduced do not impact the latency.

During O3, 56 alerts were issued on the GCN without retraction. PyCBC Live contributed to 34 of them. Only 1 of the 24 retracted alerts was produced by Py-CBC Live. The next observing run will include a brand new detector, KAGRA. Unless Virgo and KAGRA reach a sensitivity comparable to LIGO, the larger resulting detector network will probably warrant further development of the multidetector significance calculation in order to limit the impact of trials factors, as discussed earlier.

In preparation for the next observing run, we plan to investigate further improvements in the handling of instrumental transients. In particular, the inpainting method presented by Zackay et al. (2021b) could potentially broaden the applicability of gating to a wider class of glitches, if it is compatible with the latency requirements of PyCBC Live. We also plan to study the impact of applying data-quality flags to the online analysis, and to characterize the effect of the SNR maximization implemented during O3.

Our proposed rapid source classification method could also be extended to consider template parameters other than the chirp mass, although these are subject to higher statistical and systematic errors (e.g. Biscoveanu et al. 2019). A possible implementation using a larger set of triggers associated with a given astrophysical event to quantify parameter errors was shown in Stachie et al. (2021).

Finally, as the latency of the entire LIGO-Virgo public alert system keeps being reduced, further development is also under way to reduce the latency of PyCBC Live. In particular, the so-called "early-warning" detection of BNS mergers (Cannon et al. 2012;Sachdev et al. 2020;Magee et al. 2021) has been implemented in PyCBC Live as well (Nitz et al. 2020c), and will be optimized and characterized in a future study.

Figure 1 .

 1 Figure1. Visualization of the process used by PyCBC Live to generate a three-detector coincident candidate. In this example, the LIGO-Hanford and LIGO-Livingston SNR time series (red and blue curves) have two nearby peaks above the triggergeneration threshold (dashed horizontal line). The coincident peaks then produce a Hanford-Livingston coincident candidate. The window of possible signal arrival times at Virgo (on-source region, dark gray in the bottom-right panel) is calculated using the Hanford and Livingston triggers, as indicated by the horizontal arrows. The Virgo SNR time series is calculated (violet, thick curve) and searched for its maximum within the on-source region. Once the maximum is identified, its statistical significance is determined by comparing it to the maxima occurring in the off-source intervals (vertical stripes of alternating color in the top and bottom-left panels).

Figure 3 .

 3 Figure 3. Sensitivity of PyCBC Live, with the O3 configuration, for a population of simulated BNS signals added to simulated Gaussian noise at design sensitivity. The "HL" configuration corresponds to a detector network formed by LIGO-Hanford and LIGO-Livingston only. The "HLV" configuration includes Virgo. The top panel shows the relative search volume of the HLV and HL configurations. The bottom panel shows sensitivity distances for the multidetector coincidence in the HL and HLV configurations (solid lines), and for the single-detector triggering (dashed lines). The lighter bands represent the 1σ uncertainties from the Monte Carlo sampling.

Figure 4 .

 4 Figure 4. Comparison of PyCBC Live's sensitivity to BNS mergers using an O2-like ("HL") and the final-O3 ("HLV") search configurations. The simulated signals are added to real data from Advanced LIGO and Advanced Virgo's third run. The lighter bands represent the 1σ uncertainties from the Monte Carlo sampling.

  Figure5. Cumulative distribution of PyCBC Live's latency, for a period of replayed O2 data analyzed using the O3 configuration. The latency is defined here as the time elapsed between the estimated coalescence time of a candidate and the time of creation of the corresponding GraceDB event. The shaded region is the range of expected latency due to PyCBC Live alone, as described inNitz et al. (2018).

Figure 6 .

 6 Figure 6. Left: For simulated signals of each source type, we show the probability density distribution of assigned probabilities for the recovered event to be of the same (i.e. correct) type. Right: Confusion matrix for the source classification assigned the highest probability for each simulated signal.

https://gracedb.ligo.org/superevents/public/O3

https://gcn.gsfc.nasa.gov

See https://www.virgo-gw.eu/O2.html and https://www.gwopenscience.org/detector status/O2.

Note that the ranking of double-coincident triggers in each detector pair depends on the pair, since it involves the distribution of expected relative signal phases, times and amplitudes appropriate for the two detectors(Nitz et al. 2017).

Systems with chirp masses outside the considered limits are assumed to be BNS if they lie below the minimum value, and BBH if they are above the maximum.

We are grateful to Stuart Anderson, Juan Barayoga, Patrick Brockill and Sharon Brunett for computing assistance during the operation of PyCBC Live in O3. We also thank Gregory Mendell and Alan Weinstein for reviewing PyCBC Live's O3 configuration, and Francesco Pannarale for comments on the manuscript.

TDC was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by Universities Space Research Association under contract with NASA, during part of this work. GSCD, VVO and TD acknowledge financial support from the Maria de Maeztu Units of Excellence program MDM-2016-0692, from Xunta de Galicia (Centro singular de investigación de Galicia accreditation 2019-2022), and from European Union ERDF. GSCD and IWH acknowledge the STFC for funding through grant ST/T000333/1.

We acknowledge the Max Planck Gesellschaft and the Atlas cluster computing team at AEI Hannover for computing support in preparation of this draft. We are also grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.

This paper has LIGO document number LIGO-P2000296 and Virgo document number VIR-0693E-20.


Anne Hosmalin
email: anne.hosmalin@inserm.fr

Bénédicte Manoury
email: benedicte.manoury@inserm.fr

A tribute to Nilabh Shastri and a special issue on Antigen Processing and Presentation in

The 10 th Antigen Processing and Presentation Workshop was a glorious meeting, when we were all still young and unconscious. Unconscious of the Covid-19 crisis looming, we wandered happily without masks in an exciting flurry of unpublished news and informal exchanges on our preferred topic in a garden in the middle of Paris across the street from Institut Cochin. We heard short presentations from good-humored very famous scientists generously giving their time to presentations from younger scientists. We were inspired by a great scientific EMBO women in science lecture. We discussed posters in depth, and played a game on Immunology in Paris. We had very hot discussions on the relative roles of spliced peptides and other neoepitopes in the induction of T cell responses, and more [START_REF]Antigen Processing and Presentation 10 (APP 10)[END_REF]. We were proud to really feature 50% invited women speakers. We were also unconscious of Nilabh Shastri's life coming to an end so soon after, whereas he had given a mysterious title to his Keynote and had promised a future contribution to this special issue 'Immune surveillance of immune surveillance' (Figure 1). After this loss struck our APP community, the special issue linked to the Workshop turned into a Tribute to Nilabh. His trainees at the Johns Hopkins University School of Medicine give us a testimony on their privileged experience 'chasing the rainbow' with his guidance, Jon Yewdell a recollection of testimonies from 'Hoodlums', UC Berkeley students and post-docs, UC Berkeley Faculty and APP colleagues, and Robin Fahraeus, Sebastien Apcher and colleagues focus on Nilabh's pioneering discovery of non-canonical mRNA translation initiation and all its consequences in antigen presentation. We all are deeply sad from this loss of a true groundbreaking genius in our field, who also displayed all his life exceptional elegance, kindness and mentorship.

In this special issue, Julie Magarian Blander reviews the role of the Transporter Associated with antigen processing (TAP) in antigen presentation and cross-presentation, and she details non-canonical cross-presentation which relies on endoplasmic reticulum-Golgiintermediate Compartment (ERGIC)-derived Class I Molecular Histocompatibility Molecules (MHC-I), without Toll-like receptor control, thus maintaining CD8 + T cell priming even when TAP is blocked by viral escape mechanisms. Bénédicte Manoury reviews the role of endoplasmic reticulum stress in the MHC class I antigen presentation pathway of dendritic cells (DC), highlighting the role of the unfolded protein response in this pathway. Louise Boyle presents primary data where she compared MHC-I when it lacks a glycan, like in recombinant molecules derived from bacteria used for in-vitro functional assays, to physiologically glycosylated MHC-I. She found that when MHC-I lacks a glycan, it has a stronger affinity for the peptide editor TAP-binding protein related (TAPBR), thus peptides are more easily dissociated by TAPBR, impacting peptide exchange. MHC-I exerts a strong immune surveillance on tumors, which escape by multiple mutations at all possible molecular and compartment levels, the consequences of these mutations are reviewed by Marlieke Jongsma, Jacques Neefjes and Robbert M. Spaapen. Finally, as DC remain the 'professional antigen presenting cells' as coined by Ralph Steinman, Justine Mintern and Jose Villadangos discuss how the triggering of different receptors for the constant fragment of immunoglobulins (FcR), especially neonatal FcR vs FcgR, regulates the immune functions of the different DC populations, and Elodie Segura proposes a model in which conventional DC and monocyte-derived DC play complementary roles for the induction of T cell responses, by presenting antigens either in lymphoid organs or in inflamed tissues, respectively. The traditional and generous sharing of many more unpublished data of Antigen Processing and Presentation Workshop was vivid and exciting, as testified by the abstracts published in Molecular Immunology [START_REF]Antigen Processing and Presentation 10 (APP 10)[END_REF] and by articles published since.

Nilabh was at the origin of the creation of these Workshops with Betsy Mellins and Eli Sercarz in 1995, a successful series of workshops across the world with a small crowd of afficionados (Table 1). The future Antigen Processing and Presentation Workshop is yet to come, but this rendez-vous is recurrent: now from France to Australia and (unlike a boomerang) in a different location somewhere in the world in the future.

Figure 1 .

 1 Figure 1. Nilabh Shastri at APP10: 'Immune surveillance of immune surveillance'

Table 1 . List of APP workshops

 1 Institut Cochin, Inserm U1016, CNRS UMR 8104, Université de Paris, F-75014 Paris, France b Institut Necker Enfants Malades, INSERM U1151, CNRS UMR 8253, Université de Paris, F-75015 Paris, France" This is inaccordance with my co-author, Dr Bénédicte Manoury.

			Berkeley,
	1	1995	Oxnard (USA)	Nilabh Shastri, Elizabeth D Mellins, Eli Sercarz
			Bar Harbor	John D. Lich, Janice S. Blum
	2	1999	(USA)
			Vaux de Cernay	Sebastian Amigorena, Alexander V Chervonsky,
	3	2002	(France)	Alexander Rudensky, Peter Van Endert
			Bar Harbor	Alexander Rudensky, Alexander V Chervonsky
	4	2004	(USA)
			Dunk Island	Jose Villadangos, James McCluskey, William
	5	2007	(Australia)	Heath
			Cargese	Philippe Pierre, Kenneth Rock
	6	2010	(France)
			Amsterdam	Jacques Neefjes, Emil Unanue (EMBO)
			(The
	7	2012	Netherlands)
			Philadelphia	Lawrence Eisenlohr, Lisa Denzin, Paul Roche,
	8	2015	(USA)	Tim Elliott
			Salamanca	Margarita Del Val, Luis C Antón, Jack Bennink,
	9	2017	(Spain)	Vincenzo Cerundolo (EMBO)
			Paris	Anne Hosmalin, Bénédicte Manoury, Peter
	10	2019	(France)	Cresswell, Sebastian Amigorena (EMBO)

Acknowledgments

We greatly acknowledge funding by EMBO (European Molecular Biology Organization (EMBO), Université de Paris, European Federation of Immunological Societies (EFIS), Institut Cochin, DCBiol Labex, Centre National de la Recherche Scientifique (CNRS), Institut National de la Recherche Médicale (Inserm), Fondation ARC (Association de Recherche contre le Cancer), Biolegend, Becton-Dickinson, 10x Genomics, Tetramer Shop, and Molecular Immunology for publishing the abstracts and this Special Issue.


Élisabeth Remy

Paul Ruet

Average sensitivity of nested canalizing multivalued functions

We prove that the average sensitivity of nested canalizing multivalued functions is bounded above by a constant. In doing so, we introduce a generalization of nested canalizing multivalued functions, which we call weakly nested canalizing, for which this upper bound holds.

Introduction

Boolean canalizing functions are Boolean functions f : (Z/2Z) n → R such that at least one input variable, say x i (1 i n), has a value a = 0 or 1 which determines the value of f (x). Nested canalizing (NC) functions are a "recursive" version of canalizing functions: an NC function f is canalizing as above and moreover its restriction f ↾ x i =a is itself NC.

These classes of Boolean functions have been introduced by Kauffman [START_REF] Kauffman | The origins of order: Self organization and selection in evolution[END_REF][START_REF] Kauffman | Random Boolean network models and the yeast transcriptional network[END_REF] to formalize the "canalizing" behaviour observed in some discrete systems. This idea is also at the basis of Waddington's work in embryology: he described an epigenetic landscape guiding embryogenesis by canalizing configurations [START_REF] Waddington | Canalization of development and the inheritance of acquired characters[END_REF].

NC functions are particularly interesting because they have "low complexity". The average sensitivity AS(f ) (also called influence or total influence) of a Boolean function f is a measure of its complexity. It can be defined in several ways, in particular via Fourier-Walsh analysis. It is related to spectral concentration, learning properties, decision tree complexity [START_REF] Donnell | Analysis of Boolean functions[END_REF]. For arbitrary Boolean functions, AS(f ) = O(n), but some functions have significantly smaller average sensitivity. For NC functions, AS(f ) is bounded above by a constant [START_REF] Li | Boolean nested canalizing functions: A comprehensive analysis[END_REF][START_REF] Klotz | Bounds on the average sensitivity of nested canalizing functions[END_REF].

NC functions are notably used as appropriate rules in Boolean models of gene regulatory networks [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF]. The dynamical systems in biological networks are far from random, and it has been shown that NC functions ensure expected stability properties [START_REF] Kauffman | The origins of order: Self organization and selection in evolution[END_REF][START_REF] Kauffman | Genetic networks with canalyzing Boolean rules are always stable[END_REF] and are indeed predominant in large databases of Boolean gene networks [START_REF] Subbaroyan | Minimum complexity drives regulatory logic in Boolean models of living systems[END_REF].

In most cases Boolean variables are sufficient, but for some situations this description is too crude, and it may be necessary to consider other levels. To model such a situation correctly, multivalued variables have been introduced [START_REF] Thomas | Regulatory networks seen as asynchronous automata: a logical description[END_REF]. Then it is necessary to consider multivalued functions f : (Z/kZ) n → Z/kZ for some k 2. The notion of average sensitivity generalizes to the multivalued setting [START_REF] Donnell | Analysis of Boolean functions[END_REF], and multivalued NC functions are defined in [START_REF] Murrugarra | Regulatory patterns in molecular interaction networks[END_REF][START_REF] Murrugarra | The number of multistate nested canalyzing functions[END_REF]. Very little is known about their spectral properties. In [START_REF] Kadelka | Multistate nested canalizing functions and their networks[END_REF], a variant of average sensitivity, the normalized average c-sensitivity, is defined for multivalued functions, and used to measure the stability of networks made of NC functions.

A natural question is whether the average sensitivity of NC multivalued functions is bounded above by a constant, too. We prove in Theorem 3 that this is the case. We actually show that the upper bound holds for a more general class of functions, which we call weakly nested canalizing, and at the same time this enables us to establish the upper bound in a simpler way than in [START_REF] Li | Boolean nested canalizing functions: A comprehensive analysis[END_REF] for Boolean NC functions.

Nested canalizing multivalued functions

Let k, n be positive integers, k 2. Z/kZ is the ring of integers modulo k.

Following [8, 9, 1], we shall say that f : (Z/kZ) n → Z/kZ is canalizing with respect to coordinate i and (a, b) ∈ Z/kZ×Z/kZ if there exists a function g : (Z/kZ) n → Z/kZ different from the constant b such that

f (x) = b if x i = a g(x) if x i = a.
We shall simply say that f is canalizing if it is canalizing with respect to some i, a, b.

A segment is a subset of Z/kZ of the form {0, . . . , i} or {i, . . . , k -1}. Let σ ∈ S n be a permutation, A 1 , . . . , A n be segments, and c 1 , . . . , c n+1 ∈ Z/kZ be such that c n = c n+1 . Then f is said to be nested canalizing (NC) with respect to σ, A 1 , . . . , A n , c 1 , . . . , c n+1 if

f (x) =                  c 1 if x σ(1) ∈ A 1 c 2 if x σ(1) / ∈ A 1 , x σ(2) ∈ A 2 . . . . . . c n if x σ(1) / ∈ A 1 , . . . , x σ(n-1) / ∈ A n-1 , x σ(n) ∈ A n c n+1 if x σ(1) / ∈ A 1 , . . . , x σ(n-1) / ∈ A n-1 , x σ(n) / ∈ A n .
We shall simply say that f is NC if it is NC with respect to some σ, A 1 , . . . , A n , c 1 , . . . , c n+1 .

Weakly nested canalizing multivalued functions

In Theorem 3, we shall give an upper bound on average sensitivity which holds not only for NC functions, but for the more general class of weakly nested canalizing functions, which we define now.

Let n be a positive integer. For each i ∈ {1, . . . , n}, Ω i is a finite set of cardinality k i > 0, Ω = i Ω i , and f : Ω → R. Note that we do not require k i 2 for all i. If k j = 1 for some j, f could be viewed as a function with one less variable, i.e. as a function on i =j Ω i , but we still consider it as a function defined on i Ω i .

We shall say that f is weakly canalizing with respect to coordinate i and

(a, b) ∈ Ω i × R if f (x) = b whenever x i = a,
and simply that it is weakly canalizing if it is weakly canalizing with respect to some i, a, b.

Note that this definition differs slightly from the usual definition by the absence of condition on the values of f for x i = a: we do not require the existence of some x such that x i = a and f (x) = b. In particular, constant functions are weakly canalizing, though not canalizing.

If f is canalizing with respect to i, a, b and k i 2, we shall consider

f ↾ x i =a : Ω ∩ {x | x i = a} → R, the restriction of f to the set of x ∈ Ω such that x i = a.
The class of weakly nested canalizing on Ω = i Ω i is then defined by induction on the cardinality |Ω| = i k i of Ω:

• If |Ω| = 1, i.e. k i = 1 for all i, any f : Ω → R is weakly nested canalizing (WNC) on Ω. • If |Ω| > 1, f : Ω → R is WNC on Ω if it is weakly canalizing with respect to some i, a, b such that k i 2 and f ↾ x i =a is WNC on Ω ∩ {x | x i = a}, a strict subset of Ω.
Intuitively, a function f : Ω → R is WNC if its domain Ω can be "peeled" by successively removing coordinate hyperplanes (defined by equations of the form x i = a) whose points are mapped by f to the same value, whence the following characterization: Proposition 1. Letting K = i k i , f is WNC if and only if there exist a function v : {1, . . . , K} → {1, . . . , n} and numbers a i ∈ Ω v(i) and b i ∈ R for each i ∈ {1, . . . , K} such that:

f (x) =            b 1 if x v(1) = a 1 b 2 if x v(1) = a 1 , x v(2) = a 2 . . . . . . b K if x v(1) = a 1 , . . . , x v(K-1) = a K-1 , x v(K) = a K .
In decomposing an NC function f : (Z/kZ) n → Z/kZ, each coordinate i ∈ {1, . . . , n} is considered exactly once (in some order prescribed by a permutation σ) and the value of f is fixed for x σ(i) in some segment A i . This can be realized by successively fixing the value of f for each α ∈ A i , and therefore, the class of WNC functions contains the class of NC functions, as stated in the following Proposition:

Proposition 2. If f : (Z/kZ) n → Z/kZ is NC, then it is WNC.
Proof. Assume f is NC with respect to σ, A 1 , . . . , A n , c 1 , . . . , c n+1 . For each i ∈ {1, . . . , n}, let

A i = {α 1 i , . . . , α |A i | i } (Z/kZ) \ A i = {α 1+|A i | i , . . . , α k i } with α 1 i < • • • < α |A i | i and α 1+|A i | i < • • • < α k i . This defines K = nk numbers α j i ∈ Z/kZ.
For each i ∈ {1, . . . , n} and j ∈ {1, . . . , k}, let

β j i = c i if j |A i | c n+1 otherwise.
To comply with the characterization of WNC functions (Proposition 1), we relabel the numbers α j i , β j i by identifying the list

α 1 1 , . . . , α |A 1 | 1 , . . . , α 1 n , . . . , α |A i | n , α 1+|A

Examples

• As we have already observed, constant functions from (Z/kZ) n to Z/kZ are WNC but not NC.

• In decomposing a WNC function f : (Z/kZ) n → Z/kZ, it is possible to "peel" a coordinate hyperplane defined on some coordinate i (i.e. by some equation x i = a), then a coordinate hyperplane defined on j, and later a coordinate hyperplane defined on i again. This is because of the recursive definition of WNC functions, and gives more freedom in the construction of WNC functions than in the construction of NC functions.

For instance, the functions min and max : (Z/kZ) 2 → Z/kZ are not NC, as observed in [START_REF] Kadelka | Multistate nested canalizing functions and their networks[END_REF]. However, an easy induction on k shows that they are WNC. For instance, min = min k : {0, . . . , k -1} 2 → {0, . . . , k -1} is weakly canalizing with respect to 1, 0, 0, min k ↾ x 1 =0 is weakly canalizing with respect to 2, 0, 0, and min k ↾ x 1 =0,x 2 =0 is identical to the function min k-1 : {1, . . . , k -1} 2 → {1, . . . , k -1}, which is WNC.

• Also, in constructing a WNC function f : (Z/kZ) n → Z/kZ, the values a used to define f (x) for x i = a need not be extremal values (initially 0 or k -1), they can be intermediate values: 0 < a < k -1.

For instance, the function from Z/3Z to Z/3Z defined by 0 → 0, 1 → 1, 2 → 0 is not NC because it is canalizing with respect to either the intermediate value 1 (for its unique variable), or the values 0 and 2 (which do not form a segment). But any function from Z/kZ to Z/kZ is WNC.

NC functions and genetic networks

We have already mentioned that canalizing functions are significantly predominant in gene network modellings [START_REF] Subbaroyan | Minimum complexity drives regulatory logic in Boolean models of living systems[END_REF], and that networks with nested canalyzing rules are stable [START_REF] Kauffman | Genetic networks with canalyzing Boolean rules are always stable[END_REF]. In this context, one is interested in the discrete-time asynchronous evolution of the expression levels of n genes, where "asynchronous" means that at each time step, the level of at most one gene can change. Moreover the expression level of each gene belongs to a finite set, typically {0, 1} or {0, 1, 2}.

The following example is inspired from the logical modelling of the phage lambda, a biological model widely studied to understand the decision between lysis and lysogenization [START_REF] Ptachne | A genetic switch. Phage lambda and higher organisms[END_REF][START_REF] Thieffry | Dynamical behaviour of biological regulatory networks II. Immunity control in bacteriophage lambda[END_REF][START_REF] Remy | From minimal signed circuits to the dynamics of Boolean regulatory networks[END_REF]. It involves two genes, CI and Cro. CI is either expressed or not, and its expression level is therefore modelled by a Boolean variable, Cro can take 3 values {0, 1, 2}. This simple model is sufficient to display both multistability (representing lysis and lysogeny fates) and oscillations (lysogeny state) [START_REF] Remy | Graphic requirements for multistability and attractive cycles in a Boolean dynamical framework[END_REF][START_REF] Ruet | Local cycles and dynamical properties of Boolean networks[END_REF].

In state x = (x CI , x Cro ) ∈ Z/2Z × Z/3Z, the next value of CI is given by the following function f CI : Z/2Z × Z/3Z → Z/2Z:

f CI (x) = 0 if x Cro 1 1 otherwise.
For instance, in state (1, 2), the next value of CI can be 0 because f CI (1, 2) = 0, and in state (0, 2), the value of CI cannot change because f CI (0, 2) = 0. Similarly, the next value of Cro is given by a function f Cro : Z/2Z × Z/3Z → Z/3Z. However, the following two choices for f Cro :

f 1 Cro (x) =      0 if x CI = 1 1 if x CI = 0 and x Cro = 2 2 otherwise f 2 Cro (x) =      1 if x Cro = 2 0 if x Cro = 2 and x CI = 1 2 otherwise
give rise to the same asynchronous trajectories, represented by the following graph, where vertices are states of the system (x CI , x Cro ) and arrows link two consecutive states: (

These two functions do not have the same canalizing property:

f 1 Cro is NC, f 2
Cro is only WNC. Thus, in this example two functions that represent the same asynchronous dynamics do not have the same canalizing properties. The observations of [START_REF] Kauffman | Genetic networks with canalyzing Boolean rules are always stable[END_REF] can thus be an assistance for modelling the biological system, a task known to be difficult, as the number of network-compatible functions is enormous.

Average sensitivity of WNC multivalued functions

Following [10, Chapter 8], we shall take the following definition of average sensitivity. First, Fourier decomposition is generalized to non Boolean domains. Let Ω = n i=1 Ω i be as above, with |Ω i | = k i . On the vector space of real-valued functions defined on Ω, an inner product is given by f, g = E x [f (x)g(x)], where E denotes the expectation. Here, x ∈ Ω and we assume independent uniform probability distributions on the Ω i . A Fourier basis is an orthonormal basis (ϕ α ) α∈ i {0,...,k i } such that ϕ (0,...,0) = 1. It is not difficult to see that a Fourier basis always exists, although it is not unique.

Then, fix a Fourier basis (ϕ α ). The Fourier coefficients of f : Ω → R are f (α) = f, ϕ α , and E i f = α i =0 f (α)ϕ α turns out to be independent of the basis. for all i ∈ {1, . . . , n}, let the ith coordinate Laplacian operator L i be the linear operator defined by

L i f = f -E i f .
Finally, the influence of coordinate i on f is defined by Inf i [f ] = f, L i f , and the average sensitivity (also called influence or total influence) of f is then AS

[f ] = i Inf i [f ].
By Plancherel's theorem (see [START_REF] Donnell | Analysis of Boolean functions[END_REF]), we have

Inf i [f ] = α i =0 f (α) 2 = E x [Var y i [f (x 1 , . . . , x i-1 , y i , x i+1 , . . . , x n )]],
where Var denotes the variance (Var[g] = E[g 2 ] -E[g] 2 ) and y i ∈ Ω i . The above equality makes clear that the definition of influence of multivalued functions generalizes the Boolean case, as expected. Now, for an arbitrary f : Ω → [0, M ], we have

Var i [f ](x) (M/2) 2 for all i, therefore Inf i [f ] M 2 /4 for all i and AS[f ] n • M 2 /4 = O(n).
For WNC functions, this upper bound can be greatly improved. In the Boolean case, [START_REF] Li | Boolean nested canalizing functions: A comprehensive analysis[END_REF] proves (by a different method from ours) that AS[f ] 2 for NC {-1, +1}-valued functions. This bound is improved in [START_REF] Klotz | Bounds on the average sensitivity of nested canalizing functions[END_REF], where it is proved that AS[f ] 4/3. For NC functions f : {0, 1} n → {0, 1}, the result in [START_REF] Li | Boolean nested canalizing functions: A comprehensive analysis[END_REF] means

AS[f ] 1/2.
Theorem 3 generalizes this result, by establishing that, in the more general multivalued case, the average sensitivity of WNC functions is bounded above by a constant.

Theorem 3. Let Ω = n i=1 Ω i where each Ω i has cardinality k i > 0. Let f : Ω → [0, M ] and κ = max i (k i -1)/k i < 1. If f is WNC (in particular if it is NC), then AS[f ] M 2 4(1 -κ) .
Proof. We prove this by induction on i k i . If k i = 1 for all i, the inequality holds trivially: actually AS[f ] = 0. Now assume f is canalizing with respect to j, a, b, with k j 2, and let f ′ = f ↾ x j =a . Let Ω ′ be the set of x ∈ Ω such that x j = a, so that f ′ : Ω

′ → R + . The induction hypothesis is AS[f ′ ] M ′2 /(4(1 -κ ′ )), with M ′ = max x∈Ω ′ f ′ (x) = max x∈Ω ′ f (x) M κ ′ = max k j -2 k j -1 , max i =j k i -1 k i κ.
Note that the induction hypothesis implies AS

[f ′ ] M 2 /(4(1 -κ)). We shall use the notation Var i [f ](x) = Var y i [f (x 1 , . . . , x i-1 , y i , x i+1 , . . . , x n )]. Then AS[f ] = E x [ i Var i [f ](x)] and AS[f ] • i k i = x i Var i [f ](x) = x j =a
Var j [f ](x) +

x j =a

Var j [f ](x) + i =j

Var i [f ](x) since f (x) is constant when x j = a, so that Var i [f ](x) = 0 for i = j. Furthermore, Var j [f ](x) is independent of x j , and on the other hand, Var i [f ](x) = Var i [f ′ ](x) if x j = a and i = j. Thus

AS[f ] • i k i = k j • x j =a Var j [f ](x) + x∈Ω ′ i =j Var i [f ′ ](x).
Since 0 f (x) M for all x, we have Var j [f ](x) M 2 /4. Therefore

AS[f ] • i k i k j • i =j k i • M 2 /4 + x∈Ω ′ n i=1 Var i [f ′ ](x) = i k i • M 2 /4 + AS[f ′ ] • (k j -1) • i =j k i and AS[f ] M 2 4 + AS[f ′ ] • k j -1 k j M 2 4 + κ • AS[f ′ ].
To conclude the proof, it suffices to observe that AS[f ′ ] M 2 /(4(1 -κ)) implies AS[f ] M 2 /(4(1 -κ)).

In the Boolean case, κ = 1/2 and M = 1, so that the upper bound M 2 /(4(1 -κ)) equals 1/2 and the above result is a generalization of the result in [START_REF] Li | Boolean nested canalizing functions: A comprehensive analysis[END_REF]. The proof is also significantly simpler than the one in [START_REF] Li | Boolean nested canalizing functions: A comprehensive analysis[END_REF]. It can be easily checked that in the Boolean case, our argument on variance essentially amounts to compute the fraction of edges in the Hamming cube {0, 1} n which are boundary edges (i.e. edges (x, y) with f (x) = f (y)).

An obvious question is whether the bound M 2 /(4(1-κ)) can be improved for multivalued WNC, or at least NC, functions, along the lines of [START_REF] Klotz | Bounds on the average sensitivity of nested canalizing functions[END_REF].

1 | 1 ,

 11 . . . , α k 1 , . . . , α 1+|An| n , . . . , α k n as the list a 1 , . . . , a K , and by identifying similarly the list b 1 , . . . , b K . Call ϕ this relabeling, mapping r ∈ {1, . . . , K} to the pair ϕ(r) = (i, j) such that a r = α j i and b r = β j i . For instance, ϕ(1) = (1, 1) and ϕ(K) = (n, k). Then finally, a function v : {1, . . . , K} → {1, . . . , n} is defined by v(r) = σ(i) if ϕ(r) = (i, j). Then f clearly enjoys the characterization of WNC functions, with the choice of function v and numbers a r , b r .