Follow-up of 27 radio-quiet gamma-ray pulsars at 110-190 MHz using the international LOFAR station FR606

The Fermi Large Area Telescope has detected over 260 gamma-ray pulsars. About one quarter of these are labeled as radio-quiet. In the population of nonrecycled gamma-ray pulsars, the fraction of radio-quiet pulsars is higher, about one half. Most radio observations of gamma-ray pulsars have been performed at frequencies between 300 MHz and 2 GHz. However, pulsar radio fluxes increase rapidly with decreasing frequency, and their radio beams often broaden at low frequencies. As a consequence, some of these pulsars might be detectable at low radio frequencies even when no radio flux is detected above 300 MHz. Our aim is to test this hypothesis with low-frequency radio observations. We have observed 27 Fermi-discovered gamma-ray pulsars with the international LOw Frequency ARray (LOFAR) station FR606 in single-station mode. We used the LOFAR high band antenna (HBA) band (110-190 MHz). On average, we use 9 h of observation per target after the removal of affected datasets, resulting in a sensitivity for pulse-averaged flux on the order of 1-10 mJy. We do not detect radio pulsations from any of the 27 sources, and we establish stringent upper limits on their low-frequency radio fluxes. These nondetections are compatible with the upper limits derived from radio observations at other frequencies. We also determine the pulsars' geometry from the gamma-ray profiles to see for which pulsars the low-frequency radio beam is expected to cross Earth. This set of observations provides the most constraining upper limits on the flux density at 150 MHz for 27 radio-quiet gamma-ray pulsars. In spite of the beam-widening expected at low radio frequencies, most of our nondetections can be explained by an unfavorable viewing geometry; for the remaining observations, especially those of pulsars detected at higher frequencies, the nondetection is compatible with insufficient sensitivity.


Introduction
The Large Area Telescope (LAT) on the Fermi satellite has strongly increased the number of known gamma-ray pulsars (see e.g., the "Fermi 2nd Pulsar Catalog" Abdo et al., 2013, hereafter 2PC). A major update, "3PC," which is in preparation, will characterize at least 260 gamma-ray pulsars 1 . For the Fermi LAT pulsar detections, three different approaches have been used (Pletsch et al., 2012a): (1) For some pulsars known from radio observations, gamma-ray pulsations have been detected using radio ephemerides, see for example Smith et al. (2019). (2)   majority of cases, however, no radio counterpart has been detected. Approximately half of the 142 nonrecycled gamma-ray pulsars are radio-loud (S 1400 ≥ 30 µJy), the rest are radio-quiet .
If detected, radio emission from a gamma-ray pulsar provides extra information, such as the pulsar's dispersion measure (DM). Combined with a geometrical model of the Galaxy (e.g., Cordes & Lazio, 2002;Yao et al., 2017), this provides an approximate source distance. This has motivated a number of radio follow-up studies for known gamma-ray pulsars 2 .
Most of these follow-up studies have been performed between 820 and 1500 MHz. Studies at lower frequencies are interesting for two reasons. First, the spectra of most pulsars are usually well described by a power law, but different pulsars have different spectral indices (Sieber, 1973;Maron et al., 2000;Jankowski et al., 2018). Some of the gamma-ray pulsars might have a steeper than average spectrum with a very low (possibly undetected) flux density at high frequencies, but a high (and potentially detectable) flux at lower frequencies. Second, the radio beams of many pulsars are wider at low frequencies (e.g., Sieber et al., 1975), such that a beam that narrowly misses Earth at higher frequencies could still intersect Earth at lower frequencies. Such beam widening is expected, for example, from the radius-to-frequency mapping model (i.e., higher frequency emission originating closer to the neutron star than lower frequency emission, Cordes, 1978).
With beam widening in mind, Maan & Aswathappa (2014) performed follow-up observations of Fermi LAT detected pulsars at 34 MHz using the Gauribidanur telescope. However, many pulsars show a spectral break or turnover between 50 and 100 MHz (e.g., Sieber, 1973;Kuzmin et al., 1978;Izvekova et al., 1981), so that observations at 34 MHz may miss pulsars which might be detectable at frequencies above this spectral turnover. In addition, scatter broadening may lead to a pulse width exceeding the pulse period, making detection difficult. Scatter broadening is usually assumed to scale as ∝ f −4 , where f is the observing frequency, so that low-frequency observations are particularly strongly impacted. Indeed, Maan & Aswathappa (2014) did not detect any periodic signal above their detection threshold (8 σ) at 34 MHz, but they were able to determine flux density upper limits.
Flux density upper limits at 150 MHz were derived for all observable radio-quiet gamma-ray pulsars using the GMRT allsky survey TGSS (Frail et al., 2016). However, the integration time per pointing was only 15 minutes, leading to relatively weak constraints on the flux density.
For these reasons, we systematically followed up on all 27 northern sky nonrecycled radio-quiet gamma-ray pulsars with targeted observations and searched for radio pulsations at low frequencies. For this, we used the international LOw Frequency ARray (LOFAR) station FR606 in Nançay in stand-alone mode, observing in the high band antenna (HBA) frequency range . This paper is organized as follows: Section 2 discusses the interest of low-frequency observations (higher flux at low frequencies in Subsection 2.1 and beam widening at low frequencies in Subsection 2.2). Section 3 describes the methods used: We describe how the targets were selected (Subsection 3.1), how the observations were performed (Subsection 3.2), how the data were processed (Subsection 3.3), and how flux density upper limits were calculated (Subsection 3.4). In Section 4, we present the results of our observations. In Section 5, we discuss the pulsars' geometry: We first discuss the consequences of beam widening for a uniform distribution of pulsar geometries (Subsection 5.1), then explain how we determined the pulsars' geometry from the observed gamma-ray profiles (Subsection 5.2), and finally discuss the consequences of beam widening for those specific geometries (Subsection 5.3). Section 6 closes with some concluding remarks.

Radio beam basics
2.1. Radio spectra At frequencies 200 MHz, the spectra of most pulsars are well described by a single spectral index (Sieber, 1973;Maron et al., 2000;Jankowski et al., 2018). However, different pulsars have different spectral indices, with a wide spread. Throughout this work, we adopt an average spectral index of −1.6 with a standard deviation of 0.54 (Jankowski et al., 2018). We assume that the spread of spectral indices is at least partially physical; for this reason, we use their standard deviation rather than their standard error. This range of values formally encompasses those given by other works, such as Maron et al. (2000) and Lorimer et al. (1995). At frequencies below 200 MHz, the spectral index flattens and most nonrecycled pulsars show a turnover between 50 and 100 MHz (e.g., Sieber, 1973;Kuzmin et al., 1978;Izvekova et al., 1981;Bilous et al., 2016;Jankowski et al., 2018;Bilous et al., 2020;Bondonneau et al., 2020a). This is not accounted for in our calculations, and the extrapolated flux density limits at frequencies < 200 MHz are probably slightly overestimated.
With a spectral index of −1.6 ± 0.54, we can expect flux densities to be higher at 150 MHz than at 1400 MHz by 1-2 orders of magnitude. This gain is partially offset by the increased sky temperature against which the pulsar has to be detected; even so, for a comparable telescope gain, or effective area, we can expect to detect fainter pulsars. This is particularly true for those that have a steeper than average spectral index.
No matter how steep the spectrum may be, a pulsar that is intrinsically faint or too distant remains undetectable. A consensus about pulsar radio luminosity L r does not exist, but it is frequently assumed that it depends on the rotation period (P) and period derivative (Ṗ) as L r = L 0 P aṖb , albeit with large uncertainties on the parameters a and b and with a large spread, that is to say P aṖb and L r are poorly correlated. As an example, Johnston & Karastergiou (2017) found a, b such that L r =Ė 1/4 whereĖ is the spindown power. Distance d affects detectability in two ways: First, S 1400 ∝ L r /d 2 . Second, a large d implies a large dispersion measure, especially at low Galactic latitudes. For pulsars with periods 500 ms and observations at 150 MHz, we are limited to DM 500 pc cm −3 ; for higher DM values, the expected scattering time exceeds the pulsar's rotation period (see Section 3.3). Thus some pulsars are not seen even if the spectrum is steep and/or the beam broadens into the line of sight, as discussed below.

Beam widening
A simple radio beam model suffices for this work. Following Cordes (1978) or, equivalently, the Pulsar Handbook (Lorimer & Kramer, 2005, Section 3.4), radio emission is centered around the neutron star magnetic dipole axis. The dipole is inclined with an angle α from the rotation axis, and the line of Article number, page 2 of 17 sight from Earth makes an angle ζ with the rotation axis. Highenergy electrons follow the magnetic field lines and radiate at radio frequencies that depend on the field line curvature. Hence, radio emission is a cone-shaped beam of radius where r em is the height of radio emission and r LC = cP/(2π) defines the "light cylinder," the radius at which an object corotating with the neutron star with a spin period P would reach the speed of light, c. Typical values of r em = 300 km and P = 0.1 s give ρ ≈ 20 • . The radio beam sweeps Earth only if (α+ρ) > ζ > (α−ρ), making the pulsar potentially radio-loud (RL). Following Gil et al. (2002), the observed pulse has a width W given by where β = (ζ − α) is the angle between the magnetic axis and the line of sight. For RQ pulsars, |β| > ρ and W is undefined. If the beam grazes Earth, W may be so small that the radio flux integrated over the narrow pulse is below a given radio telescope's sensitivity, and if detected it is likely very faint: We call these borderline cases radio-faint (RF). From the Pulsar Handbook, ρ = A 1 f −q + A 0 or r em = B 1 f −p +B 0 , where f is the radio frequency. For positive values of the index p, Eq. (1) leads to beam widening at low radio frequencies. This is called radius-to-frequency mapping (RFM). The parameters (B 0 , B 1 , p) have different values for different pulsars. For example, Thorsett (1991) studied seven multi-component pulsars that have been observed over a wide frequency range, finding the outer components to be separated by angles of 5-10 • at 1400 MHz, but up to 30 • at 150 MHz. Similarly, Xilouris et al. (1996) studied eight pulsars. Their Figure 2 shows the profile widths to approximately double between 1000 and 100 MHz.
While W increasing at low frequency f is common, this is not the case for all pulsars. Pilia et al. (2016) determined the index δ for 100 pulsars, such that the observed width at 10% of the pulse amplitude is w 10 ∝ f δ (see their Figures 3 and 7). For ∼ 20% of the pulsars, they find δ > 0, meaning that the pulse narrows with decreasing frequency. However, they did indeed obtain δ < 0 for 80% of the pulsars and found an average value of δ ∼ −0.1 (weighted mean over all values of δ). Thus, a significant fraction of their sample roughly agrees with the expectation from simple RFM.
Based on these observations, we assume RFM to hold. More specifically, we used the parameterization of Story et al. (2007, Equations (9) and (10)), which is based on the model of Kijak & Gil (1997, 1998. In terms of the above parameterization, this is equivalent to B 0 = 0 and p = 0.26. Using P = 100 ms andṖ = 1 × 10 −15 s s −1 , this model gives a beam width of 16.8 • at 1400 MHz and 22.5 • at 150 MHz. In Section 5.2, we combine these values with geometrical constraints (α, ζ) extracted from the gamma-ray profiles to predict radio detectability. The comparison of these model results to pulsed emission discovered with FR606 (or the upper limits thereof) provides useful data for emission models.

Target selection
To a good approximation, the effective area of the FR606 antenna array is where z is the zenith angle of the source and A max eff = 2048 m 2 (van Haarlem et al., 2013, Appendix B). For the telescope's latitude (λ = 47.35 • ), a declination limit of δ > −10 • ensures cos 2 z > 0.28 at pulsar culmination, so that the telescope can be used with a meaningful sensitivity for a few hours on a given target each day. We thus selected the 20 radio-quiet pulsars in 2PC with δ > −10 • , excluding Geminga (J0633+1746), which has already been extensively explored at low radio frequencies (Maan, 2015, and references therein). In addition, we included seven pulsars discovered after 2PC (Pletsch et al., 2013;Clark et al., 2015Clark et al., , 2017Wu et al., 2018). This leaves us with a total of 27 targets, which are listed in Table 1 (column 1).
All of our targets have already been observed by radio telescopes, with observing frequencies between 34 and 2000 MHz (see Table A.1 for the full list of observations). Table 1 (columns 6-8) shows the flux density upper limits for the most constraining of the previous observations. Detections rather than upper limits are denoted as such; in those cases, we indicate the equivalent flux density, assuming W/P = 0.1, rather than the measured flux density, with the measured value for the fractional pulse width W/P (see Appendix A). Column 9 gives the corresponding flux density limit at our observing frequency of 150 MHz for an assumed spectral index of −1.6.

Observations
The observations were carried out on the International LOFAR Station in Nançay, FR606. LOFAR is fully described in Stappers et al. (2011) andvan Haarlem et al. (2013). LOFAR stations have two different frequency bands: We used the HBA band (i.e., 110-190 MHz, with a center frequency of 149.9 and a total bandwidth of 78.125 MHz), for which the station consists of 96 antenna tiles, each of which is made up of 16 dualpolarization antenna elements. The signals from individual HBA tiles were coherently summed, creating a digital telescope.
Long observing sessions were split into individual observations typically lasting one hour. Each target was observed between 7×1h and 16×1h, amounting to a total of 346h of telescope time (on average approximately 13h per target). To ensure phase coherence, the time span of our radio observations was at most 15 days for each pulsar. The analysis of test data showed that daytime observations contained intense radio frequency interference (RFI) and they did not provide data of a sufficiently good quality. To determine which fraction of the day allowed for high-quality observations, we observed a well-known pulsar (B0105+68, period P ∼ 1.07 s, DM∼ 61.06 pc cm −3 ) with the same setup and pipeline as for sources of interest. Even though the pipeline includes RFI cleaning, the pulsar was not detected in daytime observations, while it was clearly detected in all nighttime observations. For this reason, all observations of the sources of interest were taken at night.
Despite this precaution, a number of observations showed high RFI levels, and they had to be discarded subsequently (observations were discarded if RFI led to false candidates at apparent S /N > 7 even after RFI cleaning, see next section). In the end, only 255h of observations with low RFI levels were re-Article number, page 3 of 17 Table 1. Northern radio-quiet gamma-ray pulsars. Notes. Column 1: name of the pulsar (pulsars with know radio emission are denoted with * ). Columns 2 and 3: Galactic coordinates for each pulsar. Column 4: pulsar period (P). Column 5: spindown luminosityĖ. Columns 6 and 7: observing frequency and measured mean flux density (or upper limit) for the most constraining observation (assuming a spectral index of −1.6). Column 8: Reference (for columns 6 and 7). Column 9: flux density (or upper limit) from columns 6 and 7 extrapolated to an observing frequency of 150 MHz (using a spectral index of −1.6). Column 10: the number of good observations (usually one hour) used in the final analysis (i.e., with a sufficiently low RFI level). Column 11: sky temperature. Column 12: the equivalent duration of the dataset under the assumption of nominal gain, as given by Eq. (5). Column 13: flux density upper limit at 150 MHz (in the band 110-190 MHz) as determined in this work. Column 14: limit on the spectral index s derived from the most constraining previous detection (columns 6 and 7) and our flux density limit (column 13).
References. (1)  tained (apparent S /N< 7), and the final analysis was based on between 2 and 13h per target (Table 1, column 10 gives the number of observing sessions lasting typically 1 hour; column 12, the equivalent duration of a zenith observations, t eff obs , is explained in Section 3.4).

Data analysis
For the data analysis, we used the standard pulsar tools tempo2 (Hobbs et al., 2006), DSPSR (van Straten & Bailes, 2011), PSRCHIVE (van Straten et al., 2012), and COASTGUARD 3 (Lazarus et al., 2016). We folded the data at the period determined by the LAT rotation ephemeris, thus leaving the pulsar's dispersion measure (DM) as the only free parameter during the data analysis. In detail, we proceeded as described in the following.
obtained an ephemeris for each of the 27 pulsars that is valid at the epoch of the FR606 observation.
During the observation, the full data stream was split on four different data acquisition machines (each receiving one quarter of the total bandwidth). The LuMP software 4 was used to recorded the data.
After the observation, the four datasets (one per data acquisition machine) were each folded at the initial Fermi LAT period using DSPSR (van Straten & Bailes, 2011) and channelized into 2440 channels of ∼10 kHz. The number of channels was chosen as a compromise between the available computing resources and the expected loss of sensitivity due to DM smearing (the associated factor β δDM is quantified in Section 3.4). Subsequently, the borders of the bands were removed (low sensitivity of the antennas near the edge of their bandpass), so that in total 8000 useful frequency channels were kept (the frequency range of . Then, the data were cleaned of RFI using COASTGUARD (Lazarus et al., 2016). Typically, ∼ 3% of the data retained in the final step were flagged (with a maximum flagged fraction of 8.4%).
Subsequently, the data were rebinned to subintegrations of 300 s in order to reduce the data volume and speed up the postprocessing. In the next step, the data from the four acquisition machines were combined.
A few months after the observation, an improved ephemeris based on Fermi LAT data was constructed (based on the extended Fermi LAT time series). The datasets (the 300 s subintegrations) were refolded with this improved ephemeris.
Using pdmp, we searched each individual observation session (∼1 hour) for a dispersed, periodic radio signal at the Fermi LAT period. This search was performed for all pulsars, including those already detected in radio: The DM precision of the high frequency observations might not be good enough (the effect of a small DM error is 100 times stronger at 150 MHz than at 1500 MHz); also, the precise DM value could have changed since the detection. We used 3980 trial DM values between 2 and 400 pc/cm 3 in steps of 0.1 pc/cm 3 (incoherent dedispersion). For 19 of our 27 lines of sight, the Galactic electron density model NE2001 (Cordes & Lazio, 2002) predicts a maximum DM of less than 400 pc/cm 3 . For the eight remaining lines of sight (all close to the Galactic plane), the corresponding estimated scattering time at 150 MHz (assuming a scatter broadening spectral index of −4) already exceeds the pulsar's period, rendering an extension to even higher DM values useless. The result is similar for the YMW16 model (Yao et al., 2017): For 17 lines of sight, the maximum DM is below 400 pc/cm 3 , and for most of the remaining ten lines of sight, the estimated scattering time at 150 MHz exceeds the pulsar's period. The only exceptions are J1958+2846, J2021+4026, and J2030+4415, for which a maximum DM of ∼500 would have been a slightly better choice.
Manual inspection showed that all observation sessions for which pdmp produced a high signal-to-noise ratio (S/N) were in fact corrupted by remaining RFI and showed high values for a wide range of DM values. For this reason, all observation sessions with a value of apparent S /N ≥ 7 were discarded. The remaining observation sessions of each target (between two and 13 sessions depending on the target) were added together using psradd.
For the combined files (one per target), we searched again for dispersed, periodic signals using pdmp. All candidate detec-tions with a value of S /N ≥ 5 (a total of 1657 candidates) were inspected manually.
This procedure was tested and validated on a well-known pulsar (B0105+68, period P ∼ 1.07 s, DM ∼ 61.06 pc cm −3 ). We then processed the observations of the 27 pulsars of Table  1 in this way; no pulsed radio signal was detected. To quantify our nondetections, we calculated the flux density upper limits, taking the elevation of the pulsar during each observing session into account (see following section).

Flux density upper limits
In none of our observations was a pulsar radio signal detected. Upper limits for the pulse-average flux were calculated following the Pulsar Handbook: where a value of S /N min = 5 was assumed as a threshold in the S/N required for a detection, T rec = 422 K is the receiver temperature (as derived from measurements by Wijnholds & van Cappellen, 2011). The sky temperature T sky depends on the sky position of the observed target; it is taken from the sky map at 408 MHz by Haslam et al. (1982) and was scaled to our frequencies using f −2.55 (Lawson et al., 1987). The values for T sky are given in Table 1 (column 11). Furthermore, β δDM is the loss of sensitivity due to imperfect DM gridding (see below), G nom is the nominal telescope gain for a pointing near zenith, n p is the number of polarizations (in all of our observations, n p = 2), t eff obs is the effective observing time for each pulsar (see below), and ∆F eff is the effective frequency bandwidth of each observation after RFI cleaning. Finally, P is the pulsar period and W is the width of the pulse; for all pulsars, we assumed W/P = 0.1.
With our grid of trial DMs, the DM error δDM of the best trial is at most 0.05 pm/cm 3 . This leads to a slight loss in the S/N and thus in sensitivity. Using Cordes & McLaughlin (2003, eq. (12-13)), a center frequency of 149.9 MHz, a total bandwith of 78.125 MHz, and W/P = 0.1, we find values in the range from 0.75 ≤ β δDM ≤ 0.99 for our sample, depending on the pulsar's period. These values were taken into account for the calculation of the flux density limits in Table 1. For the nominal gain, we used G nom = A max eff /(2k B ) = 0.74K/Jy, which is equivalent to the effective area of A max eff = 2048m 2 of the international LOFAR station FR606 for a pointing near zenith at an observing frequency of 150 MHz (van Haarlem et al., 2013, Appendix B). As the elevation of the pulsar varies during an observation, the projected area of the antenna array varies over time, cf. eq. (3). Assuming the effective area is approximately constant over an observation of approximately one hour, the effective observing time after N obs observations (with an effective area A eff,i for the i th observation) is given by It is important to note that this implies that observations at a low elevation (and thus with A eff,i ≪ A max eff ) contribute only marginally to the combined observation. The values of t eff obs are given in Table 1 (column 12). The values of A eff,i and T sky were calculated using the LOFAR calibration tool lofar_fluxcal.py (Kondratiev et al., 2015). The software makes use of the Hamaker beam model (Hamaker, 2006, and references therein) and the mscorpol 5 package to calculate the Jones matrices for a given frequency and sky coordinates.
With this, the upper limits on the flux density, S min , were calculated, and the values are given in Table 1 (column 13). The implications are discussed in Section 4.

Results
No radio emission was detected for any of our targets. At our observing frequency of 150 MHz, we obtained flux density upper limits between 0.8 mJy and 26 mJy (the values are given in column 13 of Table 1). The large differences in the upper limits for different targets result from four effects: (a) The number of 1h observations was different for different targets; (b) during some nights, the RFI conditions were worse than during others, forcing us to remove a larger number of observations from processing; (c) the low elevation of some sources led to a low effective observing time t eff obs despite a large number of observations N obs ; and (d) depending on the direction, the background sky temperature T sky varied by more than one order of magnitude. In particular, we note that our flux density limits are less constraining for J1838−0537 (26 mJy), J1846+0919 (7.6 mJy), J1906+0722 (8.0 mJy), and J1907+0602 (7.1 mJy) than for our other targets. This is related to their proximity to the Galactic plane (Galactic coordinates are given in Table 1, columns 2 and 3), which causes a high sky temperature (T sky , Table 1 column 11), and at the same time this reduces the effective area of the telescope because of their low declination (and thus leads to a small t eff obs , column 12). We obtained our most constraining flux density limit of 0.8 mJy for J0357+3205, where the sky temperature is particularly low (far off the Galactic plane) and for which RFI conditions were good, so that we have a large number of usable observations. Table 1 shows the most constraining previous observation for each pulsar, assuming a spectral index of −1.6. The observing frequency and the upper limit for the flux density are given in columns 6 and 7, and the equivalent flux density at 150 MHz is shown in column 9.
The same data are displayed in Figure 1 6 , in which we compare our flux density upper limits (empty blue squares) to upper limits (empty red triangles) and detections (filled green circles, with a typical error bar of 30%) obtained by previous observation campaigns (all values are summarized in Appendix A and Table A.1). As can be seen in Figure 1, our nondetections are compatible with the nondetection of those 21 of our targets that have been studied with imaging observations at 147.5 MHz using GMRT (Frail et al., 2016); our observations do, however, provide more stringent upper limits (by approximately one order of magnitude). The remaining six pulsars have never before been observed in this frequency range.
To test whether our upper limits are compatible with previous observations at different frequencies, we assume a spectral index of −1.6 ± 0.54 (see Section 2.1). To compare this to previous nondetections, the dark blue lines in Figure 1 represent flux density limits equivalent to our observation, assuming an aver-age spectral index of −1.6, and the shaded areas between the light blue lines correspond to a spectral index of −1.6 ± 0.54 (see Section 2.1). Figure 1 shows that our upper limits are compatible with all previous upper limits. For J0554+3107, our upper limit is more constraining than the previous observation by a factor ∼ 2.
Five of the pulsars we observed have been detected at higher radio frequencies, but they remain undetected in our observations. In Table 1, these pulsars are denoted with * in column 1, and limits for the spectral index s are given in column 14. For these pulsars, the dark green lines in Figure 1 represent flux density limits equivalent to the most constraining detection, assuming an average spectral index of −1.6, and the shaded areas between the light green lines correspond again to a spectral index of −1.6 ± 0.54 (see Section 2.1).
For J0002+6216, radio emission has been reported at frequencies of 1400 and 2000 MHz (Clark et al., 2017;Wu et al., 2018). When combined, these two observations hint at a spectral index shallower than average (−1.2). However, as both observing frequencies are not widely separated, the spectral index is not very well constrained. Our nondetection is compatible with both previous observations and constrains the spectral index to values > −2.1.
The pulsar J0106+4855 has been detected at 820 MHz (Pletsch et al., 2012b). With our setup, a nondetection is expected as long as the spectral index is > −2.2.
J0631+0636 has been detected at 327, 1400, and 1510 MHz. The two available flux measurements hint at a spectral index that is steeper than average (−1.9). This spectral index is still compatible with our nondetection, which constrains the spectral index to values > −2.3.
For J1907+0602, very faint radio emission has been reported at 1510 MHz (Abdo et al., 2010;Ray et al., 2011). This is consistent with our nondetection and constrains the spectral index to values > −3.1.
Finally, for J2032+4127, faint radio emission has been reported at 2000 MHz (Camilo et al., 2009;Ray et al., 2011). Extrapolating their flux measurement to our observing frequency with an assumed spectral index of −1.6 gives a flux value 20% higher than the upper limit derived from our observations. Considering typical errors on flux density measurements and the variability of pulsars, our nondetection is compatible with the detection at 2000 MHz. It can also be explained by either variations in the pulsar's emission (either intrinsic or due to scintillation on a timescale larger than the duration of the observation), or a spectral break or turnover in the 100-200 MHz range. Alternatively, it could be undetectable because of a flatter spectral index than the average. Indeed, for a spectral index shallower than −1.5, which is compatible with the range of spectral indices discussed above, the extrapolated flux is compatible with our nondetection.
In all cases, our observations provide the most constraining upper limits at 150 MHz for the set of 27 pulsars we observed. Based on the comparison above, J0631+0636 and J2032+4127 would be the most interesting targets for potential reobservations at 150 MHz.  Fig. 1. Flux density upper limits for 27 radio-quiet gamma-ray pulsars, based on observations with LOFAR/FR606. Empty red triangles (with arrows): Flux density limits from previous observations. Filled green circles: Flux density measurements from previous observations (a typical error bar of 30% is indicated). Empty blue squares (with arrows): Flux density limits from this study. Dark blue line: Equivalent flux density limits for our observation (assuming a spectral index of −1.6). Shaded area between the light blue lines: Equivalent flux density limits for our observation for a spectral index of −1.6 ± 0.54. Dark green line: Same as the dark blue line, but with respect to the most constraining detection. Shaded area between the light green lines: Same as the shaded area between the light blue lines, but with respect to the most constraining detection.

Beaming, revisited
Article number, page 7 of 17 A&A proofs: manuscript no. 40841final range (150 MHz). The angle ζ can have values in the range from 0 • ≤ ζ ≤ 90 • . For any given value of α, the pulsar is geometrically visible (potentially RL) if the line of sight is not too far from the magnetic axis, (α − ρ) < ζ < (α + ρ). Uncertainties in α and ζ (∆α and ∆ζ) are accounted for, such that the pulsar is only classified as RL if this condition is satisfied for all values of α ± ∆α and ζ ± ∆ζ. If ζ < α − ρ or α + ρ < ζ (again, accounting for the uncertainties in α and ζ), the pulsar is considered as RQ. In all other cases (the pulsar could be either RL or RQ depending on the precise values of α and ζ within the uncertainties), the pulsar is considered RF. Figure 2 shows the (α, ζ) plane. The values shown for the pulsars are explained in the next section. A few of these pulsars should be geometrically visible. This is expected, for example, for the pulsar with (α = 85 • , ζ = 84 • ) in Figure 2b. As the beam is wider at a low frequency (ρ 150 > ρ 1400 ), a pulsar can be RQ at 1400 MHz, but RL at 150 MHz. If we assume a fully uniform distribution of α and ζ and no uncertainties (∆α = ∆ζ = 0), approximately 20% of the pulsars that are RQ at 1400 MHz could be RL at 150 MHz. For our sample of 27 pulsars, this could mean on the order of five detections at 150 MHz with FR606.
Of course, a uniform distribution of α and ζ is not expected, and the uncertainties of α and ζ have to be taken into account. In the following, we estimate the values of α and ζ for the 27 pulsars of our sample, based on their gamma-ray profiles, and we attempt to improve our estimate of how many pulsars may have been detectable.

Pulsar geometry from gamma-ray profiles
Gamma-ray beams are very narrow in neutron star longitude due to concentration of the gamma-radiating electrons and positrons along "caustically" focused magnetic field lines; however, they are very broad in latitude, being brightest near the neutron star equator, and fading toward the poles. Beam shapes depend on several parameters: The open field line configuration varies with α and r LC , and electron acceleration gap sizes depend on magnetic field strength. Gamma-ray emission models (see below) differ in how they exploit these parameters, but all have in common that "in fine" an observed profile depends on α and ζ. In the following we use the observed gamma-ray profiles to constrain the pulsar geometry (in terms of α and ζ), and thus refine our prediction of how many, and which, of our pulsars may be radio-detectable at low frequency.
We first made weighted gamma-ray pulse profiles for our sample, integrated above 100 MeV. Phases were calculated using tempo2 and rotation ephemerides were derived from LAT data for use in 3PC using the methods of Kerr et al. (2015). Weights, that is the probability that a given photon comes from the pulsar and not from a background source, were calculated using the methods of Kerr (2011). In all cases, the lightcurves were compatible with those previously published, such as in 2PC.
The weight calculation requires a spectral and spatial analysis of the region surrounding the pulsar. For our analysis, we used the results of the FL8Y source list 7 applied to Pass 8 (P8R2) LAT data corresponding to events within 3 • of each pulsar spanning 2008-08-04 to 2016-08-04. Our timing solution for J1932+1916 folds poorly beyond 2013-07-04. and so we used just under 5 years of data for that pulsar. Our LAT data sets include events belonging to the SOURCE class as defined under the 7 See https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/. P8R2_SOURCE_V6 instrument response functions 8 , with energies between 0.1 and 100 GeV and zenith angles ≤ 90 • .
We then compared the gamma-ray profiles with the predictions from the outer gap (OG) model and the two-pole caustic (TPC) model, generated over a broad parameter space. Specifically, we used the same simulations and pulse-profile fitting techniques as Wu et al. (2018) (for more details on the simulation and fitting see Johnson et al., 2014). References to the OG and TPC models are given there. For J0002+6216, J0106+4855, J0631+0646, J1907+0602, and J2032+4127, we jointly fit the gamma-ray and radio profiles, using the 1400, 820, 1400, 1400, and 2000 MHz pulse profiles, respectively. As in Wu et al. (2018), we could not produce acceptable fits for J0631+0646. Compared to Wu et al. (2018), our analysis includes more gamma-ray photons, leading to updated results, especially for J0002+6216 and J2017+3625. Table 2 lists the resulting α and ζ angles, with columns 2 and 3 for the OG model, and columns 4 and 5 for the TPC model. These results are also shown in Figure 2. Most of the fits are limited by systematic uncertainties of ∼ 10 • (see Johnson et al., 2014;Pierbattista et al., 2015, for more details).
From the best-fit geometries and using the radio cone model of Story et al. (2007), we predicted the radio detectability for each pulsar at 1400 MHz (Table 3, columns 2 and 4 for models TPC and OG) and 150 MHz (Table 3, columns 3 and 5 for models TPC and OG, respectively). RL means that our line of sight should intersect a bright part of the cone. RF suggests that our line of sight skims the cone edge, and RQ means that the cone completely misses Earth. Aligned and nearly aligned rotators are considered as RF.
The characteristic width of the simulated radio beam at 0.1% of the peak intensity is given by Equations (9) and (10) of Story et al. (2007). For our simulations, we used P = 100 ms andṖ = 1 × 10 −15 s s −1 , which gives a width of 16.8 • at 1400 MHz and 22.5 • at 150 MHz.
Column 6 in Table 3 reports (− ln(likelihood) TPC − (− ln(likelihood) OG ). When this value is positive, the OG model describes the data better than the TPC model, and vice versa. Based on their experience fitting many pulsars, Johnson et al. (2014) found that an absolute value of ≥15 was needed in log likelihood difference to determine that one model was significantly favored over another. Based on this, the last column gives the preferred model for each pulsar (indicated in parentheses if the log likelihood difference is <15). Similar to the conclusions drawn by Johnson et al. (2014), Pierbattista et al. (2015), and others, we did not find that one model predominantly describes gamma-ray pulse profiles better than another. Table 2 shows our results for the geometric angles ζ and α for the models TPC and OG. The implications for (geometrical) detectability of radio emission at 150 MHz and 1400 MHz are shown in Table 3. The results are also shown in Figure 2, where pulsars detected in radio are shown with large filled circles and pulsars without detected radio emission are shown with small empty circles.

Geometry results and discussion
In Table 3, pulsars with detected radio emission are denoted with a * in column 1. It can be seen that all of these are labeled as either RF or RL. Equivalently, all detected pulsars fall within or close to the shaded region in Figure 2.  Pulsars on the dashed diagonal line (α = ζ) beam their radio emission to Earth regardless of their beam opening angle. Radio emission from pulsars between the two colored dashed lines is only detectable if their respective beam width ρ ≤ 16.8 • , whereas pulsars between the two dash-dotted lines are visible if ρ ≤ 22.5 • . Pulsars detected in radio are shown with large filled circles; pulsars without detected radio emission are shown with small empty circles.
Only one of the pulsars labeled as RL at 1400 MHz (J2238+5903) has not been detected in radio. However, a nondetection does not rule out the presence of low-level radio emission. 27 pulsars at 150 MHz is compatible with an average spectral index and with RFM-like beam-widening for most pulsars.

Conclusion
We have followed up on 27 radio-quiet gamma-ray pulsars at low radio frequencies . No pulsed radio emission was detected. We have established stringent upper limits on their low-frequency radio flux density, which are considerably more constraining than previous limits in comparable frequency ranges.
Despite the beam-widening expected at low radio frequencies, the comparison to simulations shows that most of our non-detections can be explained by an unfavorable viewing geometry; for the remaining observations (especially those of pulsars detected at higher frequencies), the sensitivity of our setup was not sufficient.
Based on our geometrical simulations and flux density upper limits, our best candidates for follow-up observations at 150 MHz are J0633+0632 and J2032+4127. It should be noted that follow-up observations of pulsars detected in radio, even if only at a higher frequency, are easier to analyze as no search for the correct DM value is required. With this in mind, all five pulsars already detected in radio should be reobserved at 150 MHz, potentially except for J1907+0602 for which a low-frequency Article number, page 10 of 17 detection is only possible if the spectrum turns out to be exceptionally steep.
Acknowledgements. LOFAR, the Low-Frequency Array designed and constructed by ASTRON, has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy. Nançay Radio Observatory is operated by Paris Observatory, associated with the French Centre National de la Recherche Scientifique and Université d'Orléans. We acknowledge the use of the Nançay Data Center computing facility (CDN -Centre de Données de Nançay). The CDN is hosted by the Station de Radioastronomie de Nançay in partnership with Observatoire de Paris, Université d'Orléans, OSUC and the CNRS. The CDN is supported by the Region Centre Val de Loire, département du Cher. 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. Work at the Naval Research Laboratory is supported by NASA.
A&A proofs: manuscript no. 40841final Appendix A: List of radio observations Table A.1 List of all previous radio observations for our 27 targets. The (upper limit) spectra of Figure 1 are based on the data in this table, as are columns 6-9 in Table 1. For detected pulsars, Figure 1 and Table 1 both use the nominal values (assuming W/P = 0.1). Column 1: pulsar name. Column 2: observing frequency. Column 3: pulsar detected in observation? ("y" if detected, blank otherwise). Column 4: mean flux density or upper limit. For uniformity, this is the equivalent flux density (assuming W/P = 0.1). When available, the measured flux density (with the measured values of P and W) is given in a footnote. Column 5: telescope name. Column 6: frequency bandwidth of the observation. Column 7: observing time t (for observations of this work, we give the effective observing time t eff obs instead, see eq. (5). Column 8: minimum S/N assumed by the authors. Column 9: duty cycle (W/P) assumed by the authors.