Extending the Event-weighted Pulsation Search to Very Faint Gamma-ray Sources

TL;DR

This paper introduces two novel methods to enhance pulsation detection in very faint gamma-ray sources by scanning spectral parameters, thereby expanding our understanding of gamma-ray pulsar populations despite source confusion and limited spectral information.

Contribution

The paper presents two new techniques that enable pulsation searches for very faint gamma-ray sources without prior spectral knowledge, improving detection sensitivity.

Findings

Both methods successfully detect pulsations in all tested pulsars.

The methods work even for pulsars with no significant phase-averaged emission.

They reduce the need for detailed spectral models in pulsation searches.

Abstract

Because of the relatively broad angular resolution of current gamma-ray instruments in the MeV-GeV energy range, the photons of a given source are mixed with those coming from nearby sources or diffuse background. This source confusion seriously hampers the search for pulsation from faint sources. Statistical tests for pulsation can be made significantly more sensitive when the probability that a photon comes from the pulsar is used as a weight. However, the computation of this probability requires knowing the spectral model of all sources in the region of interest, including the pulsar itself. This is not possible for very faint pulsars that are not detected as gamma-ray sources or whose spectrum is not measured precisely enough. Extending the event-weighted pulsation search to such very faint gamma-ray sources would allow improving our knowledge of the gamma-ray pulsar population. We present two methods that overcome this limitation by scanning the spectral parameter space, while minimizing the number of trials. The first one approximates the source/background ratio yielding a simple estimate of the weight while the second one makes use of the full spatial and spectral information of the region of interest around the pulsar. We test these new methods on a sample of 144 gamma-ray pulsars already detected by the Fermi Large Area Telescope data. Both methods detect pulsation from all pulsars of the sample, including the ones for which no significant phase-averaged gamma-ray emission is detected.

Figures

Tables

Table 1. Table 1: Results of the TS < 25 TS 25 \mathrm{TS}<25 pulsars: TS of the spectral fit, model and simple weight pulsation significances, P m subscript 𝑃 𝑚 P_{m} and P s subscript 𝑃 𝑠 P_{s} .

Pulsar	$(l,b)$	TS	$P_m$	$P_s$
J0922+0638	(225.4,36.4)	15.2	9.9	8.4
J1224$-$6407	(300.0,-1.4)	24.3	16.7	13.8
J1455$-$3330	(330.7,22.6)	24.8	15.4	12.7
J1513$-$5908	(320.3,-1.2)	2.7	140.5	144.0
J1646$-$4346	(341.1,1.0)	1.0	22.8	15.0
J1739$-$3023	(358.1,0.3)	2.1	7.9	6.1
J1801$-$2451	(5.3,-0.9)	1.0	153.4	125.4
J1831$-$0952	(21.9,-0.1)	8.8	15.6	10.7
J1832$-$0836	(23.1,0.3)	20.5	7.7	8.0
J1856+0113	(34.6,-0.5)	22.9	66.7	61.3
J1909$-$3744	(359.7,-19.6)	18.2	9.8	6.8
J2317+1439	(91.4,-42.4)	22.3	5.4	4.5

Table 2. Table 2: Pulsars with TS > 25 TS 25 \mathrm{TS}>25 and P m / P fit > 1.2 subscript 𝑃 𝑚 subscript 𝑃 fit 1.2 P_{m}/P_{\mathrm{fit}}>1.2 . The reported spectral index corresponds to the power-law fit, expect for the two σ curv > 4 subscript 𝜎 curv 4 \sigma_{\mathrm{curv}}>4 pulsars (indicated with a †), for which the spectral index is the result of the fit with the power law with an exponential cutoff.

Pulsar	TS	$P_m$	$P_{\mathrm{fit}}$	$\gamma$	$S_f$
J0729$-$1448	62.8	24.1	19.7	$2.5\pm 0.1$	4.7
J0742$-$2822†	556.2	273.8	219.0	$1.7\pm 0.2$	3.4
J1151$-$6108	150.0	39.3	31.6	$2.4\pm 0.1$	4.6
J1410$-$6132	116.9	33.1	26.5	$2.7\pm 0.1$	4.4
J1431$-$4715	76.4	10.2	8.1	$2.7\pm 0.1$	4.3
J1531$-$5610	155.9	40.8	32.3	$2.6\pm 0.1$	4.1
J1828$-$1101†	251.4	9.0	6.1	$2.1\pm 0.2$	4.0
J1935+2025	89.2	15.4	9.2	$2.7\pm 0.1$	4.3

Equations

$$H_{mw}=\max \left[Z_{iw}^2-c\times (i-1)\right],1\le i\le m$$

$$Z_{mw}^2=\frac{2}{{\sum }_{i=1}^N{w}_i^2}\times \sum _{k=1}^m({\alpha }_{wk}^2+{\beta }_{wk}^2)$$

$${\alpha }_{wk}=\sum _{i=1}^N{w}_i\cos (2\pi k{\varphi }_i)$$

$${\beta }_{wk}=\sum _{i=1}^N{w}_i\sin (2\pi k{\varphi }_i),$$

$$w(E,\text{@vec}\Omega )=\frac{r_{\mathrm{psr}}(E,\text{@vec}\Omega )}{r_{\mathrm{psr}}(E,\text{@vec}\Omega )+r_{\mathrm{bkg}}(E,\text{@vec}\Omega )}.$$

$$w(E,\text{@vec}\Omega )=r_{\mathrm{psr}}(E,\text{@vec}\Omega )/r_{\mathrm{bkg}}(E).$$

$$w(E,\text{@vec}\Omega )=f(E)\times g(E,\text{@vec}\Omega )$$

$$K(x,s)=\frac{1}{4\pi s^2}{\left(1+\frac{x^2}{4s^2}\right)}^{-2}$$

$$g(E,x)={\left(1+\frac{9x^2}{4{\sigma }_{\mathrm{psf}}^2(E)}\right)}^{-2}$$

$${\sigma }_{\mathrm{psf}}(E)=p_0(E/100{)}^{p_1}\oplus p_2$$

$$\mathrm{d}N/\mathrm{d}E\propto (E/E_0{)}^{-\gamma }{e}^{-(E/E_{\mathrm{c}}{)}^{\beta }}$$

$$f(E)=e^{-0.5({\log }_{10}E-{\mu }_w{)}^2/{\sigma }_w^2}.$$

$$w(E,x,{\mu }_w)=e^{-2({\log }_{10}E-{\mu }_w{)}^2}{\left(1+\frac{9x^2}{4{\sigma }_{\mathrm{psf}}^2(E)}\right)}^{-2}$$

$$P_s=\underset{{\mu }_w\mathrm{scan}}{\max }{P}_w-{\log }_{10}6$$

$$\mathrm{d}N/\mathrm{d}E\propto (E/E_0{)}^{-\gamma }{e}^{a(E_0^{\beta }-E^{\beta })}$$

$$P_m=\underset{L_m\mathrm{scan}}{\max }{P}_w-{\log }_{10}6$$

$$S_f=0.576(\gamma +\delta \gamma )+0.817{\log }_{10}(E_c+\delta E_c).$$

$${\lambda }_1(N)={\lambda }_0+0.0525796e^{-N/215.170}+0.086406e^{-N/35.5709}$$

$$\log_{10}P(H_{20}>x)\sim\left\{\begin{array}[]{ll}\lambda_{0}x&\text{if }x<15,\\ 15\lambda_{0}+\frac{\lambda_{0}+\lambda_{1}(N)}{2}(x-15)&\text{if }15<x<29\\ 22\lambda_{0}+\lambda_{1}(N)(x-22)&\text{if }x>29.\end{array}\right.$$

$$\log_{10}P(H_{20w}>x)\sim\left\{\begin{array}[]{ll}\lambda_{0}x&\text{if }x<15,\\ 15\lambda_{0}+\frac{\lambda_{0}+\lambda_{1}(W+5)}{2}(x-15)&\text{if }15<x<29\\ 22\lambda_{0}+\lambda_{1}(W+5)(x-22)&\text{if }x>29.\end{array}\right.$$

Full text

*Aims. *Statistical tests for pulsation can be made significantly more sensitive when the probability that a photon comes from the pulsar is used as a weight. However, the computation of this probability requires knowing the spectral model of all sources in the region of interest, including the pulsar itself. This is not possible for very faint pulsars that are not detected as gamma-ray sources or whose spectrum is not measured precisely enough. Extending the event-weighted pulsation search to such very faint gamma-ray sources would allow improving our knowledge of the gamma-ray pulsar population.

*Methods. *We present two methods that overcome this limitation by scanning the spectral parameter space, while minimizing the number of trials. The first one approximates the source/background ratio yielding a simple estimate of the weight while the second one makes use of the full spatial and spectral information of the region of interest around the pulsar.

*Results. *We test these new methods on a sample of 144 gamma-ray pulsars already detected by the Fermi Large Area Telescope data. Both methods detect pulsation from all pulsars of the sample, including the ones for which no significant phase-averaged gamma-ray emission is detected.

1 Introduction

One of the important contributions of the Fermi Large Area Telescope (LAT) (Atwood et al., 2009) to gamma-ray astronomy is the detection of many pulsars. Prior to the Fermi launch in 2008, at most 10 gamma-ray pulsars were known (Thompson, 2008). Ten years later, more than 200 have been detected111https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars, representing the most abundant Galactic source class at GeV energies (Acero et al., 2015). Detecting more pulsars, especially faint ones, is very important to characterizing as well as possible the gamma-ray pulsar population.

One way of finding gamma-ray pulsars is to search for periodic emission after phase-folding gamma-ray data with rotation ephemerides, obtained from contemporaneous timing observations conducted in radio or X-rays. Prior to Fermi’s launch, a Pulsar Timing Consortium (PTC) was organized to support LAT observations of pulsars by providing accurate pulsar ephemerides (Smith et al., 2008). Another way of finding gamma-ray pulsars is to look at LAT catalog unassociated sources with a pulsar-like spectrum (see e.g Clark et al., 2017; Pleunis et al, 2017) but, by construction, it does not target very faint gamma-ray sources. We thus focus on the former approach in the rest of the paper.

The sensitivity to pulsation depends on the purity of the event sample, which depends on the level of background and on the instrument performance, especially its angular resolution. In the MeV-GeV energy range, the Galactic diffuse emission is bright and is the main source of background for candidate pulsars near the Galactic plane. Neighboring sources can also significantly increase the level of background because of the finite angular resolution of the instrument. The LAT Point Spread Function (PSF) at low energy is driven by multiple scattering of the electron and positron created by pair-conversion of the photon in the tracker. The 68% containment angle ranges from 5.5 degrees at 100 MeV to 0.1 degrees at 30 GeV (Atwood et al., 2009).

As shown in Bickel et al. (2008) and Kerr (2011), it is possible to efficiently mitigate the background problem by weighting the pulsation statistical test with the probability of each event to come from the pulsar position. To compute the exact probability, one needs a complete description of the Region of Interest (RoI) around the pulsar (large enough to fully take into account the effects of the PSF), including the spectral model of all gamma-ray sources in the RoI. So one prerequisite is to be able to measure the spectrum of the candidate pulsar, which is usually done by performing a maximum-likelihood fit of the RoI. When the fit leads to a clear detection of the pulsar, its spectrum is well measured and the photon probabilities can be computed. On the contrary, a weak detection of the pulsar or obviously its non detection precludes a good measurement of its spectrum. As a consequence, current event-weighted methods to search for pulsation cannot be applied to very faint gamma-ray sources.

The only possible way to overcome this limitation is to scan over the pulsar spectral parameters in order to find the ones maximizing the statistical test used to search for periodic emission. The obvious drawback is that the resulting significance must be corrected for the number of trials performed during the scan. Consequently, the challenge of extending the event-weighted methods to very faint pulsars consists of minimizing the number of trials when scanning the spectral parameter space, while maximizing the pulsation sensitivity, i.e. getting as close as possible to the spectral parameters yielding the optimal weights.

After presenting in Section 2 the statistical test we use for pulsation searches, we present in Section 3 a method using a simple weight definition that does not involve a detailed spatial nor spectral modelling of the RoI. In Section 4, we describe a more complex method that takes advantage of the complete knowledge of the RoI and the full capabilities of the instrument.

Both methods are tested on a sample of 144 known LAT pulsars. This sample comprises the 117 pulsars from the LAT Second Gamma-ray Pulsar Catalog (2PC) (Abdo et al., 2013), among which only two are not significant gamma-ray sources when analysing 8 years of LAT data. Adding 27 post-2PC detected pulsars (Hou et al., 2014; Laffon et al., 2014; Smith et al., 2017) to our sample increases the number of faint pulsars to 12. The results, based on the analysis of 8 year of Pass 8 SOURCE class data (Atwood et al., 2013; Bruel et al., 2018), are presented and discussed in Section 5.

2 The event-weighted H-test

To test the signal periodicity, we use the weighted H-test derived by Kerr (2011) from the original H-test proposed by de Jager et al. (1989):

[TABLE]

The weighted H-test definition involves the weighted $Z_{mw}^{2}$ test:

[TABLE]

where

[TABLE]

and

[TABLE]

where the sums run over the list of $N$ photons with pulsar rotational phase $\phi_{i}$ and weight $w_{i}$.

In the following, we adopt the H-test parameter recommended values $m=20$ and $c=4$ (since they provide an omnibus test) and write $H_{w}$ instead of $H_{20w}$. We note that an obvious property of $Z^{2}_{w}$ and $H_{w}$ is that they are not sensitive to a global scaling of the photon weights.

The analytic, asymptotic null distribution of $H_{w}$ has been derived by Kerr (2011). Because we want to use $H_{w}$ even with small event samples, we derive this distribution directly with Monte Carlo simulations, allowing us to know, for any $x>0$, the probability $P(H_{w}>x)$, that $H_{w}$ evaluated on a non-periodic sample can fluctuate to values larger than $x$. The calibration procedure is described in Appendix A and Appendix B. We show that $P(H_{w}>x)$ can be parameterized as a function of the sum of the weights, computed under the prescription that the maximum weight is 1.

Rather than using directly $H_{w}^{\mathrm{meas}}$, the measured value of $H_{w}$, when searching for pulsation, we use instead $P_{w}=-\log_{10}{P(H_{w}>H_{w}^{\mathrm{meas}})}$, that we name pulsation significance for simplicity’s sake. We note that the relation between $P_{w}$ and $n_{\sigma}$, the significance expressed in number of $\sigma$, is given by $P_{w}=-\log_{10}\left(1-\mathrm{erf}(n_{\sigma}/\sqrt{2})\right)$. The 3, 4 and 5$\sigma$ levels correspond to $P_{w}\sim 2.57,4.20$ and 6.24, respectively.

3 Simple weights

Let us consider an RoI centered on a pulsar. The probability that a photon originates from the pulsar depends on the position and spectrum of all gamma-ray sources in the ROI and the response functions of the instrument. We assume that all sources are steady and that the pulsar spectrum does not depend on phase. For a photon at a position $\@vec{\Omega}$ with an energy $E$, the probability that it comes from the pulsar is the ratio of the pulsar differential rate, $r_{\mathrm{psr}}$, and the total differential rate at this energy and position. Defining $r_{\mathrm{bkg}}$ as the sum of all but the pulsar differential rates, we have:

[TABLE]

3.1 Simple weight definition

Since we are interested in very faint pulsars, we assume that the pulsar differential rate is negligible compared to $r_{\mathrm{bkg}}$. Let us also assume that the background only comprises an isotropic diffuse emission. The background differential rate thus depends only on the energy $E$ and we have:

[TABLE]

We note that the resulting weight $w$ is directly proportional to the pulsar absolute flux. Since the event-weighted statistical tests are unchanged when all weights are scaled by a constant, in the limit of very faint pulsars only the shape of the pulsar spectrum matters and not its absolute flux normalization. We also note that the weight position dependence only comes from the pulsar differential rate. Since the pulsar is a point source, the position dependence is described by the instrument PSF.

We can rewrite the weight as the product of two functions:

[TABLE]

where $g(E,\@vec{\Omega})$ contains all the position dependence and is defined such that it is 1 at the pulsar position at all energies. $f(E)$ is the weight at the pulsar position and depends on the pulsar and background spectral shapes and on the energy dependent part of the PSF normalization. The maximum of $f(E)$ is set to 1 to follow the maximum weight prescription mentioned in Section 2.

At a given energy, the LAT PSF is well approximated by a Moffat profile (Moffat, 1969)222Formerly referred to as King profile in https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_PSF.html:

[TABLE]

where $x$ is the angular distance to the source and $s$ a scale parameter.

The integral $\int_{0}^{\infty}\int_{0}^{2\pi}K(x,s)x\mathrm{d}x\mathrm{d}\theta=1$ and we note that the integral up to a distance of $3s$ is about 0.68, so that the parameter $s$ corresponds to a third of the PSF 68% containment angle. We can thus write the function $g(E,\@vec{\Omega})$ as

[TABLE]

where $\sigma_{\mathrm{psf}}(E)$ is the energy dependent PSF 68% containment angle, which for the LAT Pass 8 SOURCE event class can be parameterized as follows (Ackermann et al., 2013):

[TABLE]

with $p_{0}=5.11$ degrees, $p_{1}=-0.76$, $p_{2}=0.082$ degrees and $E$ in MeV and the addition is in quadrature.

The weight position-independent part, $f(E)$, corresponds to the weight at the position of the pulsar. Its computation is made simple by the fact that both the pulsar and background rates involve the same instrument effective area, which cancels out in Equation 6. As a consequence $f(E)$ can be computed directly using the pulsar and background spectral shapes, as well as the energy dependent part of the PSF normalization ($1/4\pi s^{2}$ in Equation 7).

As shown in Abdo et al. (2013), the spectral shape of pulsars between 60 MeV and 60 GeV energy range can be modeled with a power law with an exponential cutoff:

[TABLE]

with, for most of pulsars, an index $\gamma$ ranging from 0.5 to 2, an energy cutoff $E_{\mathrm{c}}$ between 0.6 and 6 GeV (and $\beta$ set to 1 in 2PC). We consider three $(\gamma,E_{\mathrm{c}})$ test-cases: (2.0, 0.6 GeV), (0.5, 0.6 GeV) and (2.0, 6 GeV). The last two encompass the bulk of the pulsars, while the first one corresponds to pulsars whose pulsation is driven by low energy photons.

Regarding the background spectral shape, we use the spectrum of the Galactic diffuse emission towards the Galactic center. It can be modeled with a smoothly broken power law, with spectral indices of $\sim 1.6$ and $\sim 2.5$ below and above $\sim 3$ GeV, respectively333The spectral indices are estimated using the Fermi-LAT Galactic diffuse emission model gll_iem_v06.fits available at https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html..

Figure 1 shows $f(E)$ for the three pulsar test-cases. Expressed as a function of $\log_{10}E$, the function $f$ is gaussian-like. The pulsar energy cutoff is responsible for the decrease of $f$ at high energy. At low energy, its decrease is mainly driven by the increase of the instrument PSF. The center position parameter $\mu_{w}$ is located between 2.5 and 4.5 (corresponding to 30 MeV and 30 GeV, respectively) while the half 68% width parameter $\sigma_{w}$ ranges from 0.3 to 0.45.

This result does not change qualitatively when using the spectrum of the Galactic diffuse emission at high Galactic latitude: for a given pulsar spectrum, changing the background spectral index moves the peak position of $f(E)$ but its shape remains gaussian-like, with a width in the same [0.3,4.5] range. As a consequence, we can use the same weight definition for young pulsars, that are mostly near the Galactic ridge, and millisecond pulsars, that tend to lie at high Galactic latitude.

We approximate $f(E)$ with a Gaussian in $\log_{10}E$:

[TABLE]

We note that the examples in Figure 1 are not exactly gaussian. Using other functional forms that better model the negative tail of $f(E)$ does not improve the performance of the pulsation search, that is, better matching of $f(E)$ does not compensate for the simple assumptions that we use (very faint source on top of an isotropic background). Similarly, in Section 4.1 we show that refining the PSF information fails to improve the performance.

The gaussian definition of $f(E)$ introduces two parameters, $\mu_{w}$ and $\sigma_{w}$. Formally, one would have to scan over these two parameters in order to fully explore the pulsar spectral parameter space. However, $\sigma_{w}$ varies much less than $\mu_{w}$, allowing us to fix $\sigma_{w}$ and to vary only $\mu_{w}$. We note that the $\sigma_{w}$ range was obtained under the very faint pulsar hypothesis. For brighter pulsars, the energy range over which the pulsar emission is significant is larger, leading to a wider shape for $f(E)$, which corresponds to larger values of $\sigma_{w}$. For this reason, we choose to fix $\sigma_{w}$ to 0.5. We check that, on average, this choice is valid for all pulsars in Section 5.

Thus, the simple weight definition that we propose is:

[TABLE]

3.2 Simple weight scan

For a given value of $\mu_{w}$, we can compute $H_{w}$ and the corresponding pulsation significance $P_{w}$ using the weights $w(E,x,\mu_{w})$. Figure 2 shows $P_{w}$ as a function of $\mu_{w}$ for PSR J1646$-$4346 (Smith et al., 2017). We chose this pulsar to illustrate the method because it is one of the pulsars in our sample that are not detected as a gamma-ray source with 8 years of LAT data and, as a consequence, for which the original Kerr (2011) method can not be applied.

As can be seen in this example, $P_{w}(\mu_{w})$ has a gaussian shape around its maximum. This shape is simply the result of the relative matching between the weights $w(E,x,\mu_{w})$ and the optimal weights, when $\mu_{w}$ varies. We name $\mu_{P}$ and $\sigma_{P}$ the maximum position and the 68% width of $P_{w}(\mu_{w})$. Figure 3 shows the distribution of $\sigma_{P}$ versus $\mu_{P}$ for the 144 pulsar sample. The bulk of the pulsars have $3<\mu_{P}<4$ while a few have a low $\mu_{P}$, corresponding to a pulsation signal driven by low energy photons (e.g. PSR J1513$-$5908, aka PSR B1509$-$58 (Abdo et al., 2010; Kuiper et al., 2017), whose energy cutoff is below 30 MeV). The width $\sigma_{P}$ is between 0.35 and 0.9.

Searching for pulsation consists in finding the maximum of $P_{w}(\mu_{w})$. The easiest way to do so is to perform a very fine scan over $\mu_{w}$ but that would lead to a too large number of trials. Taking advantage of the gaussian-shape of $P_{w}(\mu_{w})$ around its maximum, we perform the following algorithm (illustrated in Figure 2):

Test three values of $\mu_{w}=(2,3,4)$. Let $\mu_{0}$ be the one giving the maximum $P_{w}$. 2. 2.

Test two more values of $\mu_{w}=(\mu_{0}-0.5,\mu_{0}+0.5)$. The 0.5 distance is chosen such that it is of the order of the minimum of the 68% width of $P_{w}(\mu_{w})$. Let $\mu_{1}$ be the one giving the maximum of $P_{w}$ among all the trials with $2\leq\mu\leq 4$; 3. 3.

Test $\mu_{w}=\mu_{g}$, the position of the maximum of the gaussian passing through the three points $\mu_{1}-0.5,\mu_{1},\mu_{1}+0.5$ (all tested in the previous steps);

The resulting pulsation significance of the simple method, corrected for the 6 trials, is:

[TABLE]

We note that the 6 trials are partially correlated, especially the ones with a small difference in $\mu_{w}$. Unfortunately there is no simple way to estimate and take into account these correlations when correcting the pulsation significance and so we simply ignore them and adopt the conservative definition of $P_{s}$ given by Equation 12.

4 Model weights

To increase the sensitivity of the pulsation search, one has to go beyond the simple weight definition of Equation 11. There are two possible approaches:

• •

use the spatial and spectral description of all sources in the RoI;

• •

take into account the full PSF information.

While the former requires precise RoI modelling including spectral fits, the latter can be achieved by a minimal change of the weight definition, as it is a consequence of a refinement of the PSF definition.

4.1 PSF event types

One of the new features of Pass 8 (Atwood et al., 2013), the latest version of the reconstruction and selection of LAT data, is the introduction of the PSF event types. Prior to Pass 8, the events of a given class (e.g. SOURCE class) were divided into two subclasses, FRONT and BACK, depending on the location of the photon pair-conversion in the LAT tracker. Because of the different thicknesses of the converters in these two sections of the tracker, the PSF of the two subclasses are significantly different: at 1 GeV the FRONT and BACK PSF 68% containment angles are 1.2 and 0.61 degrees, respectively. In Pass 8, thanks to a multivariate analysis, the events are divided into four event types PSF0, PSF1, PSF2 and PSF3, ordered in decreasing quality of the reconstructed direction, with a PSF 68% containment angle at 1 GeV of 1.8, 1.0, 0.66 and 0.42 degrees, respectively.

It is straightforward to take into account the PSF events types in the context of the simple weight definition of Equation 11: one has just to properly modify the $f(E)$ part of the weight induced by the $1/4\pi s^{2}$ factor of the normalization of the Moffat profile (Equation 7). We note that this modification is made simple by both the faint-source and uniform-background hypotheses.

Tested on LAT data, this modification failed to clearly increase the pulsation significance of known pulsars. This failure can be explained by two reasons. First, as soon as the faint-source hypothesis does not hold, the correction induced by the $1/4\pi s^{2}$ factor of the Moffat profile normalization is not correct: in the extreme situation in which the background is negligible, all photons should have a weight $\sim 1$, regardless of their PSF event type.

Moreover, taking into account the better PSF description given by the PSF event types is in fact equivalent to a better spatial description of the RoI. This improvement might not be large enough to compensate the crude approximation of the uniform-background hypothesis. In other words, one has to use a precise spatial and spectral description of the RoI to fully take advantage of the additional information provided by the PSF event types.

As a consequence, going beyond the simple weights presented in the previous section requires the full spatial and spectral description of the RoI around the pulsar.

4.2 Pulsar RoI description

To compute the weights for each pulsar, we perform a binned maximum-likelihood fit, following the standard Fermi-LAT analysis procedure, of the $10.05\times 10.05$ degrees RoI centered on the pulsar, with a pixel size of 0.05 degrees. We use 37 bins in $\log_{10}E$, between 63.1 MeV and 316.2 GeV.

We include in the source model the following components:

• •

the Galactic diffuse emission and the isotropic template (that accounts for the isotropic diffuse emission as well as the residual background) 444https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html;

• •

all point-like and extended sources from the preliminary 8-year Fermi-LAT source list (FL8Y) 555https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/, within 5 degrees of the RoI border.

If the closest point source to the pulsar is within a distance of 1.5 times the source 95% error radius, we consider the source to be the pulsar and set its position to the pulsar position. Otherwise we add a new source at the pulsar position.

We perform the maximum-likelihood fit with SOURCE class events (combining all PSF event types) above 160 MeV whose zenith angle is less than 90 degrees to avoid Earth limb contamination, collected between 2008 August 4 and 2016 August 2. The significance of each source in the model is estimated using the Test Statistic, TS, defined as twice the difference in log-likelihood obtained with and without the source. A $\mathrm{TS}=25$ corresponds to $\sim 4\sigma$ significance (Mattox et al., 1996). The 160 MeV energy threshold for the maximum-likelihood fit has been chosen to mitigate the effect of systematic errors due to our imperfect knowledge of the Galactic diffuse emission.

Using the model derived by the fit, we build for each PSF event type the 3-D map (sky position and energy) of the number of predicted events coming from the pulsar as well as the 3-D map of the total number of predicted events. We are then able to compute the weights by simply dividing the 2 maps.

When computing $H_{w}$, we assign to each event the weight of the bin corresponding to the event position and energy of the 3-D map corresponding to the PSF event type of the event.

In Kerr (2011), the weights are computed using the Fermi-LAT Science Tool gtsrcprob666https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/overview.html that assigns to each event the probability that the event belongs to a given source of the RoI. Our binned approach is less CPU intensive and the chosen binning is fine enough to not induce any significant loss of sensitivity.

4.3 Spectral parameter scan

For bright gamma-ray pulsars, the maximum-likelihood fit gives a large TS for the pulsar and its spectral parameters are well estimated. We can thus use these spectral parameters to compute the weights, as done in Kerr (2011). On the contrary, for faint pulsars (TS¡25) or pulsars just above the TS=25 threshold for which the spectrum is not precisely estimated, we have to scan over the spectral parameters.

Instead of using Equation 9 to model the power law with an exponential cutoff, we use the following expression:

[TABLE]

as used in FL8Y and the forthcoming 4FGL catalog, with $\beta$ given by the fit for bright pulsars or fixed to $2/3$ for faint or not-detected ones. The formal change between the energy cut off and parameter $a$ is convenient to scan the spectral parameter space, as will be shown later.

For a given set of $(a,\gamma)$, we fix the corresponding spectral parameters and perform the maximum-likelihood fit with the pulsar normalization being the only free parameter. If $\mathrm{TS}<4$, we set the normalization to the 68% confidence limit. Using the resulting pulsar spectrum, we compute the four PSF event type weight maps and compute the pulsation significance, taking into account the PSF event type information by using for each event the PSF event type weight map corresponding to the event. Figure 4 shows how the pulsation significance varies in the $(a,\gamma)$ plane for PSR J1646$-$4346.

The iso-level contours have an ellipse-shape. Let us consider the high $P_{w}$ region corresponding to $P_{w}$ larger than 90% of the maximum, whose ellipse-like contour is named ${\cal E}_{90}$. To characterize ${\cal E}_{90}$, we use the lowest and highest $a$ extremities of ${\cal E}_{90}$. These points are shown in Figure 5 for the 144 pulsar sample. The important result is that the lowest and highest $a$ points of ${\cal E}_{90}$ are not mixed and it is possible to define a line separating them.

The fact that ${\cal E}_{90}$ is clearly elongated along a direction that tends to go from the upper-left to the lower-right of the $(a,\gamma)$ plane is the consequence of how the weights vary with the spectral parameters. For a given background spectrum, when $a$ increases, the effective position of the spectrum cutoff decreases and the energy position of the maximum weight at the pulsar position shifts to lower energy. On the other hand, the maximum weight energy position shifts to higher energy when $\gamma$ decreases. So when $a$ increases and $\gamma$ decreases, the two effects partly counterbalance each other and the maximum weight energy position does not vary very much. The situation is reversed along the minor axis of ${\cal E}_{90}$: the increases of $a$ and $\gamma$ both shift the maximum weight energy position to lower energy. This is why ${\cal E}_{90}$ is very eccentric.

We note that the convenient elliptical shape of the high $P_{w}$ region is the result of the change from Equation 9 to Equation 13 to model the power law with an exponential cutoff.

4.4 Model weight scan

As for the simple weight method, the goal is to find the maximum of $P_{w}$ in the lowest possible number of trials. The characterization of the high $P_{w}$ region obtained in the previous section naturally suggests a 2-step procedure: first finding the major-axis of ${\cal E}_{90}$ then scanning along the major-axis to find the maximum $P_{w}$.

To find the major-axis of ${\cal E}_{90}$ we perform a scan along the line $L_{m}$, going from ($a=0,\gamma=-1$) to ($a=0.03,\gamma=3$), that crosses almost all ${\cal E}_{90}$, as shown in Figure 5. The choice of $L_{m}$ is motivated by the fact that increasing both $a$ and $\gamma$ at the same time allows a fast variation of the maximum weight energy position and, therefore, an efficient exploration of the pulsar spectrum parameter space. This choice might be slightly refined in the future thanks to the analysis of a larger sample of gamma-ray pulsars.

An example of the $L_{m}$ scan is shown in Figure 6. Because the $P_{w}$ variation along $L_{m}$ is gaussian-like around the maximum, we can perform a similar 6-trial algorithm as for the simple weight method, with the following test positions:

Test three values of $a=(0.005,0.015,0.025)$. Let $a_{0}$ be the one giving the maximum $P_{w}$. 2. 2.

Test two more values of $a=(a_{0}-0.005,a_{0}+0.005)$. The 0.005 distance is chosen such that it is of the order of the minimum of the 68% width of $P_{w}(a)$ when scanning along $L_{m}$. Let $a_{1}$ be the one giving the maximum of $P_{w}$ among all the trials with $0.005\leq a\leq 0.025$; 3. 3.

Test $a=a_{g}$ the position of the maximum of the gaussian passing through the three points $a_{1}-0.005,a_{1},a_{1}+0.005$ (all tested in the previous steps); 4. 4.

Let $A_{m}(a_{m},\gamma_{m})$ be the point giving the maximum $P_{w}$ among the 6 trials.

$A_{m}$ is expected to lie on the major-axis and we need one more point to define it. Figure 7 shows $\gamma_{0}$, the intercept of the major-axis with the $\gamma$-axis, as a function of $a_{m}$. The correlation between the two can be modeled with $\gamma_{0}=0.4+100a_{m}$. We use this correlation to choose the point $A_{0}(a,\gamma)=(0,0.4+100a_{m})$ that, together with $A_{m}$, defines the major axis $L_{M}$.

An example of the $P_{w}$ variation along $L_{M}$ is shown in Figure 6. The $L_{M}$ scan allows us to reach a higher $P_{w}$ than the $L_{m}$ scan, but, because the latter goes through ${\cal E}_{90}$, the gain in $P_{w}$ is rather modest (less than 10% by definition of ${\cal E}_{90}$). To optimize the definition of the second step, we look at the maximum relative gain in $P_{w}$ along $L_{M}$ versus its relative position with respect to $a_{m}$. As can be seen in Figure 8, the gain in $P_{w}$ is on average about 1%, reaching at most 5% for few pulsars. Regarding the optimal position $a$, it can be either smaller or larger than $a_{m}$. This would imply at least two more trials, on top of the six ones already performed during the $L_{m}$ scan, inducing an additional trial correction of $\log_{10}(8/6)=0.125$.

For pulsars on the verge of pulsation detection, i.e. at the $4\sigma$ level, corresponding to $P_{w}\sim 4.2$, the additional trial correction corresponds to about 2% of $P_{w}$, larger than the average 1% gain of the $L_{M}$ scan. As a consequence, we choose not to perform the $L_{M}$ scan. The model weight pulsation significance corrected for the 6 trials is thus:

[TABLE]

As in the simple weight method, the 6 trials are partially correlated and, for the same reasons, we choose to ignore them and to use a conservative definition of $P_{m}$.

5 Results

To test the performance of the simple and model weight methods, we apply them to the sample of 144 LAT pulsars (117 pulsars from 2PC and 27 post-2PC detected pulsars). For each pulsar, we perform the binned maximum-likelihood fit presented in Section 4.2 to estimate the TS of the pulsar and compute the model weights. 12 pulsars are found to have $\mathrm{TS}<25$.

To perform the pulsation search, we select Pass 8 SOURCE class events within 5 degrees of the pulsar above 60 MeV. We use the Fermi plugin (Ray et al., 2011) to the pulsar timing software Tempo2 (Hobbs et al., 2006; Edwards et al., 2006) and the pulsar ephemeris provided by the PTC (Smith et al., 2008, 2018) to convert the event arrival time into a pulsar rotational phase. The pulsation probabilities are computed using only the data collected during the validity period of the ephemerides.

We also compare the simple and model weight results to the one obtained with the original Kerr (2011) method, which simply corresponds to the model weight method with the weights computed using the spectral parameters given by the maximum-likelihood fit. This is possible only when the source is significantly detected by the fit, i.e. when $\mathrm{TS}>25$. We name $P_{\mathrm{fit}}$ the corresponding pulsation significance. Compared to $P_{m}$, $P_{\mathrm{fit}}$ has the obvious advantage of being the result of only one trial.

To compare the performance of the new methods to a standard unweighted pulsation search, we also perform a grid search, with each grid point corresponding to a set of cuts in energy and distance to the pulsar. We use the grid parameters defined by Kerr (2011), where it is reported that the event weighting method improves the pulsation sensitivity by a factor 1.5–2.

5.1 Model weights

The pulsation significance $P_{m}$ obtained with the model weights is shown in Figure 9 as a function of TS. All 144 pulsars have $P_{m}>5.4$, corresponding to $4.6\sigma$. The general trend is that $P_{m}$ on average increases with TS. This is not surprising since, on average, the larger the TS, the stronger the pulsation signal. The scatter around this trend is mainly due to the pulse shape diversity. The most interesting result is the clear detection of pulsation from all pulsars with $\mathrm{TS}\lesssim 50$, which proves that the model weight method is able to detect pulsation even when the pulsar spectrum information is not available or not fully reliable. The $P_{m}$ and $P_{s}$ values are given in Table 1 for the 12 pulsars with $\mathrm{TS}<25$.

Ignoring the PSF event type information when computing the model weights leads to a loss of sensitivity that decreases with TS, from 20% on average for faint pulsars to 10% for the brightest.

The comparison between $P_{m}$ and $P_{\mathrm{fit}}$ is shown in Figure 10. On average the model weight method provides the same pulsation significance, within $\leavevmode\nobreak\ \sim$5% for most of the pulsars. Only one has $P_{m}$ lower than $P_{\mathrm{fit}}$ by more than 10%.

On the contrary, the model weight method often improves significantly over $P_{\mathrm{fit}}$, with a gain greater than 20% for 8 pulsars. These 8 pulsars all have $\mathrm{TS}>60$ so it is unlikely that the improvement is explained by a poor estimation of the spectral parameters by the fit. A possible explanation is the presence of a significant off-pulse emission: the spectral parameters derived by the fit correspond to the spectrum of the sum of the pulsed and unpulsed gamma-ray components and not to the pulsed component alone, while the model weight method is able to find spectral parameters closer to the pulsed component spectrum.

As reported in Abdo et al. (2013), 34 2PC pulsars have a significant off-peak emission, whose spectrum is compatible with a simple power-law spectrum for 13 of them. Moreover, the off-peak emission of these 13 pulsars is generally soft, with a spectral index $\sim 2$. A flat and soft off-peak emission that is not negligible compared to the pulsed emission could lead to a total emission spectrum significantly different from the pulsed emission spectrum. Investigating that it is also the case of the largest $P_{m}/P_{\mathrm{fit}}$ pulsars would require performing the analysis of the off-pulse emission, which is out of the scope of this paper and will be presented in the forthcoming third LAT Gamma-ray Pulsar Catalog (3PC).

We instead estimate the significance of the curvature of the spectrum. As in Abdo et al. (2013) and Acero et al. (2015) we perform a maximum-likelihood fit assigning the pulsar a power-law spectrum and compute $\sigma_{\mathrm{curv}}=(\mathrm{TS}-\mathrm{TS}_{\mathrm{PL}})^{-1/2}$, where TS and $\mathrm{TS}_{\mathrm{PL}}$ correspond to the maximum-likelihood fit of a power-law with and without exponential cutoff, respectively. We note however that $\sigma_{\mathrm{curv}}$ measures the significance of the curvature and not how flat or soft the spectrum is.

A very simple and crude estimator of the softness and flatness of the spectrum is provided by the sum $S=\gamma+\log_{10}E_{c}$, with $E_{c}=a^{-1/\beta}$. This estimator can be improved by taking into account the correlation between $\gamma$ and $\log_{10}E_{c}$. The analysis of the distribution of $\gamma$ vs $\log_{10}E_{c}$ for the 144 pulsar sample yields a positive correlation with a slope $\sim 1.418$ and the projection along the first principal axis is $S_{f}=0.576\gamma+0.817\log_{10}E_{c}$. In the case of a spectrum without significant curvature, the spectral parameters $\gamma$ and $E_{c}$ are not precisely measured. To take into account this case, we replace $\gamma$ and $E_{c}$ by their 68% upper limit and define the following simple “soft–and–flat” estimator:

[TABLE]

Figure 11 shows the ratio $P_{m}/P_{\mathrm{fit}}$ as a function of $S_{f}$. Although the correlation between these two quantities is not perfect, we note that all the 8 pulsars with $P_{m}/P_{\mathrm{fit}}>1.2$ lie in the right-hand tail of the $S_{f}$ distributions, meaning that their spectra are among the most soft and flat of the sample. Some of their properties are given in Table 2. Only two of them have a curvature significance above $4\sigma$: J0742$-$2822 and J1828$-$1101 with $\sigma_{\mathrm{curv}}=8$ and 4.9, respectively. All the others have $\gamma>2.4$. We note that only four of these 8 pulsars are in 2PC and only two are detected as a gamma-ray source (J0742$-$2822 and J1410$-$6132). Out of these two, only J1410$-$6132 is reported to have a significant off-peak emission.

The comparison of the model weights with the grid search is shown in Figure 12 as a function of TS for the 144 pulsar sample. We confirm the improvement of the weighting method over the unweighted approach, with a gain in sensitivity larger than 2 for the $\mathrm{TS}<25$ pulsars.

5.2 Simple weights

The comparison between the simple weights and the model weights is shown in Figure 13. As expected, the simple weights are less powerful than the model weights. The difference in performance decreases from about 30% for the brightest pulsars to an average of about 15% at $\mathrm{TS}\sim 300$. This difference never goes beyond 40%. This difference is relatively small given the simplicity of the simple weight implementation compared to the complex procedure of the model weights. The $P_{s}$ values are given in Table 1 for the 12 pulsars with $\mathrm{TS}<25$. Except for a few of the brightest pulsars, the simple model performs better than the standard grid search, especially for $\mathrm{TS}\leq 100$ where $P_{s}$ is greater than $P_{\mathrm{grid}}$ by at least 40%.

When defining the simple weights in Section 3.1, we set $\sigma_{w}=0.5$. To validate this choice, we perform a scan over $\sigma_{w}$ between 0.2 and 0.95, with a 0.05 step. For each pulsar and each value of $\sigma_{w}$, we compute $r(\sigma_{w})$, the ratio of $P_{s}$ divided by the maximum $P_{s}$ reached over the $\sigma_{w}$ scan. The average over the 144 pulsar sample of this ratio as a function of $\sigma_{w}$ is shown in Figure 14, as well as the lowest ratio, i.e. corresponding to the pulsar for which the choice of $\sigma_{w}$ is the worst. The two quantities have a maximum between 0.5 and 0.6 and do not vary much around the maximum, which validates the $\sigma_{w}=0.5$ choice.

To summarize the results, we show in Figure 15 the comparison of all methods that allows a clear ranking of the pulsation search methods, from the least sensitive unweighted grid scan to the most powerful Kerr (2011) and model weight methods, with the latter, contrary to the former, being able to detect pulsation from very faint pulsars.

We note that both the simple and model weight approaches can be used at other wavelengths (e.g. X-ray and TeV bands) but they need to be adapted to the specific context of each wavelength, taking into account the pulsar spectral parameter phase space, the PSF energy dependence and the typical background spectrum. In the case of the simple weight method, the derived general shape of $f(E)$ may be very different from the gaussian-like one obtained for the LAT energy band.

6 Conclusion

We have shown that it is possible to extend the event weighting technique to very faint gamma-ray sources when searching for pulsation and presented two approaches. The first one, the simple weight method, uses a very simple definition of the weights while the second one, the model weight method, fully takes into account the spatial and spectral information of the RoI around the pulsar. The key point for both methods is to explore efficiently the pulsar spectral parameter space, which is done with only 6 trials.

The model weight method reaches the same perfomance as the original event weighting method that was not applicable to the very faint gamma-ray sources. It can even be more powerful in the case of pulsars with a significant off-pulse emission.

The simple weight method is less sensitive than the model weight one but the loss of sensitivity is only $\sim$30% for faint sources. So this simple approach can be very useful, especially since it is straightforward to implement and is much less CPU intensive.

As a consequence the model weights method is on average the most sensitive method. It works for all pulsars (faint or bright, with or without off-pulse emission) and naturally benefits from any improvement of the instrument performance (e.g. Pass 8 PSF event types).

After ten years in orbit, Fermi-LAT continues to take data that are folded with the updated ephemerides provided by the Pulsar Timing Consortium. The two new methods presented in this paper are designed to help to detect new gamma-ray pulsars, including very faint ones, allowing us, for instance, to further test the existence of a spin-down power “deathline” (Smith et al., 2018) below which pulsars might cease to produce gamma-ray emission.

Acknowledgements.

We thank our Fermi-LAT collaborators Toby Burnett, Lucas Guillemot, David Smith and Matthew Kerr for fruitful discussions. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). The Parkes radio telescope is part of the Australia Telescope 533 which is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO. We thank our colleagues for their assistance with the radio timing observations. The Lovell Telescope is owned and operated by the University of Manchester as part of the Jodrell Bank Centre for Astrophysics with support from the Science and Technology Facilities Council of the United Kingdom. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

Appendix A Monte Carlo estimated probability distribution of $H_{20}$

The asymptotic $H_{20}$ probability distribution has been estimated using Monte Carlo (MC) by de Jager et al. (1989) and de Jager et al. (2010). It corresponds to a simple exponential $P(H_{20}>x)=e^{-\lambda x}$ with $\lambda\sim 0.4$. This result is very close to the analytic solution found by Kerr (2011) that yields a practical formula (valid for $m\geq 10$) with $\lambda=0.398405$.

Pulsation searches very often require multiple trials corresponding to different event selections. Some of these can lead to small data samples. As a consequence it is important to estimate the $H_{20}$ probability distribution for low $N$, the number of events used when computing $H_{20}$. To do so, we ran $10^{9}$ MC realizations for various $N$, assigning a random phase to each event. Figure 16 shows the logarithm of the cumulative distribution, $\log_{10}P(H_{20}>x)$, for $N=20,50,100$ and 1500. The asymptotic behavior is reached in $N=1500$ case but the distribution of the other cases shows a clear departure from the pure exponential above $x\sim 20$.

To characterize this departure, we fit $\log_{10}P(H_{20}>x)$ with a linear polynomial over the interval corresponding to $-7<\log_{10}P(H_{20}>x)<-4$ (i.e. $x\gtrsim 23$). The slope $\lambda_{1}$ of this linear polynomial is shown in Figure 17 as a function of $N$. We find the following parameterization for $\lambda_{1}$ (shown in Figure 17):

[TABLE]

where $\lambda_{0}=-0.398405/\log(10)=-0.173025$.

We first approximate $\log_{10}P(H_{20}>x)$ with a broken linear polynomial whose slopes at small and large $x$ are $\lambda_{0}$ and $\lambda_{1}(N)$, respectively. We find that the break position is compatible with $x=22$ for all $N$. For low values of $N$ the broken linear polynomial approximation does not work well around the break position, as shown in Figure 18. To have a better representation of $\log_{10}P(H_{20}>x)$, we use a double broken linear polynomial approximation, with the same slopes at small and large $x$ as for the single one, with break positions at $x=15$ and $x=29$:

[TABLE]

When considering the $x$ interval such that $P(H_{20}>x)>10^{-7}$, the maximum of the absolute difference between the parameterization and the MC result is less than 0.1 for $N\geq 20$. For $N=10$, the maximum difference reaches 0.25 at $x=29$. As a consequence, the approximation given by Equation 17 is valid for $N\geq 20$ with a 0.1 precision for $P(H_{20}>x)>10^{-7}$.

Appendix B Monte Carlo estimated probability distribution of $H_{20w}$

To estimate the probability distribution of $H_{20w}$, we ran $10^{8}$ MC realizations of a 5 degree region of interest with a uniform background following the Galactic diffuse emission spectrum defined in Section 3.1 (corresponding to a broken power law, with spectral indices of $\sim 1.6$ and $\sim 2.5$ below and above $\sim 3$ GeV, respectively), convoluted with the Fermi-LAT effective area. These realizations are performed for various numbers of events. For each realization, we compute $H_{20w}$ using the simple weights definition of Equation 11 with $\sigma_{\mathrm{psf}}$ set to 1 degree independently of energy. For CPU efficiency’s sake, each realization is used twice, with $\mu_{w}$ set to 2.5 and 3, each choice leading to a value of $H_{20w}$.

When parameterizing $\log_{10}P(H_{20}>x)$ in the previous Section, we use the number of events because it drives the level of fluctuations. In the case of the weighted version of $H_{\mathrm{test}}$, the number of events is not useful anymore. We use instead the sum of the weights, $W$, that we compute under the prescription that the maximum weight is 1. Figure 19 shows $\log_{10}P(H_{20w}>x)$ for some of the MC configurations corresponding to $W\sim 20,50,100,700$. As in the unweighted case, $\log_{10}P(H_{20w}>x)$ departs from a pure exponential at low $W$ and can be approximated to first order by a broken linear polynomial. We note that the various MC configurations (number of events, $\mu_{w}$ choice) were chosen to explore a range for $W$ between about 10 and 1000.

We estimate $\lambda_{1}$, the slope over the interval corresponding to $-7<\log_{10}P(H_{20w}>x)<-4$. Figure 20 shows how the $\lambda_{1}$ variation with $W$ compares to the unweighted-case $\lambda_{1}$ parameterization of Equation 16. This parameterization works well for large values of $W$ but not at low $W$. We find that using $W+5$ instead of $W$ gives a better match, as shown in Figure 21. The fact that increasing $W$ leads to a better match is not surprising since, in the weighted case, the number of events that play a significant role in the $H_{20w}$ computation (i.e. with a relatively large weight) is on average larger than $W$.

As a consequence, we use the following parameterization of $\log_{10}P(H_{20w}>x)$, that is obtained by simply replacing $N$ by $W+5$ in Equation 17:

[TABLE]

When considering the $x$ interval with $P(H_{20w}>x)>10^{-7}$, the maximum of the absolute difference between the parameterization and the MC result is less than 0.1 for $W\geq 10$. As a consequence, the approximation given by Equation 18 is valid for $W\geq 10$ with a 0.1 precision for $P(H_{20w}>x)>10^{-7}$.