Pulsars with NenuFAR: backend and pipelines

NenuFAR (New extension in Nan\c{c}ay upgrading LoFAR) is a new radio telescope developed and built on the site of the Nan\c{c}ay Radio Observatory. It is designed to observe the largely unexplored frequency window from 10 to 85\,MHz, offering a high sensitivity across its full bandwidth. NenuFAR has started its"early science"operation in July 2019, with 58\% of its final collecting area being available. Pulsars are one of the major topics for the scientific exploitation of this frequency range and represent an important challenge in terms of instrumentation. Designing instrumentation at these frequencies is complicated by the need to compensate for the effects of both the interstellar medium and the ionosphere on the observed signal. Our real-time pipeline LUPPI (Low frequency Ultimate Pulsar Processing Instrumentation) is able to cope with a high data rate and to provide real-time coherent de-dispersion down to the lowest frequencies reached by NenuFAR (10\,MHz). The full backend functionality is described, as well as the main pulsar observing modes (folded, single-pulse, waveform, and dynamic spectrum). This instrumentation allowed us to detect 172 pulsars in our first targeted search below 85\,MHz, including 10 millisecond pulsars (6 of which detected for the first time below 100 MHz). We also present some of the"early science"results of NenuFAR on pulsars: a high frequency resolution mapping of PSR B1919$+$21's emission profile and a detailed observation of single-pulse sub-structures from PSR~B0809$+$74 down to 16\,MHz, the high rate of giant-pulse emission from the Crab pulsar detected at 68.7\,MHz (43 events/min), and the illustration of the very good timing performance of the instrumentation, allowing us to study dispersion measure variations in great detail.


Introduction
Pulsars are rapidly rotating, highly magnetised neutron stars, which emit collimated beams of radiation across the electromagnetic spectrum. Most of known pulsars have been discovered and studied in the radio band. According to the ATNF Pulsar Catalogue 1 (Manchester et al. 2005), about 2800 pulsars have been found to date. Most of them have been discovered in the L-band (1-2 GHz) with, e.g. the High Time Resolution Universe survey -HTRU (Keith et al. 2010;Cameron et al. 2020), the Parkes Multi beam Pulsar Survey -PMPS (Manchester et al. 2001;Lorimer et al. 2015), and the Pulsar Arecibo L-band Feed Array survey -PALFA (Lazarus et al. 2015;Parent et al. 2019)). A few hundred sources were unveiled by searches around 300-400 MHz (Hessels et al. 2008;Lynch et al. 2013), and recently the low frequency domain was explored intensively with LOFAR, discovering 73 new pulsars in the window 119-151 MHz with the LOTAAS survey (Sanidas et al. 2019).
In the radio band, pulsars are known to have relatively steep spectra, with power-law index in the range of −1 to −3 (Malofeev et al. 2000;Maron et al. 2000;Bilous et al. 2016). Thus, despite the observational challenges caused mostly by deleterious effects of propagation in the ionised interstellar medium, low-frequency pulsar surveys are promising. New low-frequency discoveries help sample the low luminosity end of the pulsar population and explore better the properties of the ionised interstellar medium, to which low frequencies are most sensitive. In addition, measuring pulsar flux densities at low radio frequencies helps constraining the so-called low-frequency turnover (Sieber 1973;Malofeev 1993;, which is often present in pulsar spectra. Finally, the properties of radio emission at low frequencies become very dynamic, providing unique insight into conditions in the pulsar magnetosphere. Several wideband pulsar studies have been performed recently at frequencies below 100 MHz using pioneering radio telescopes such as LOFAR-LBA (Hassall et al. 2012;Pilia et al. 2016;Kondratiev et al. 2016;, UTR-2 (Zakharenko et al. 2013) and the LWA (Stovall et al. 2015). LOFAR-LBA has a large fractional bandwidth, but highly non-uniform sensitivity, with antenna response peaking in the 45-75 MHz range. UTR-2 operates at very low frequencies, 10-30 MHz, but records a single linear polarisation. The LWA has proven to be very powerful, sharing the same antenna radiator design as NenuFAR (Hicks et al. 2012) with a homogeneous frequency coverage across the band, but with only 256 antennas (to be compared with the present 1064 optimised NenuFAR antennas, that will eventually increase to 1824).
NenuFAR 2 is a compact phased-array and interferometer, formed of hexagonal groups of 19 antennas called Mini-Arrays (MA). Its core counts today 56 MA = 1064 antennas. It will be extended to 80 MA = 1520 antennas by the end of 2020, and eventually reach 96 MA = 1824 antennas (Zarka et al. 2020). Among current low-frequency telescopes, NenuFAR is well equipped to undertake an exhaustive census of the pulsar population in the low radio frequency band.

Observing pulsars with NenuFAR
Similarly to LOFAR's high frequency antennas (HBA, see van Haarlem et al. 2013), NenuFAR is designed to deal with the received radio signal on three different levels. While single antennas see the entire sky above them, analogically phased MAs of 19 antennas have a field of view of 8 to 69 degrees depending on the frequency (85-10 MHz). This can be computed from the full width at half maximum (FWHM) λ D , with λ the wavelength and D the MA diameter ≈25 m. The full NenuFAR core with its 400 m diameter, allows to synthesise digital beams which are 0.5 to 4.2 degrees wide. The geocentric coordinates of the phase centre are the same as for LOFAR international station FR606: (x,y,z) = (4324017.054, 165545.160, 4670271.072 m). The time reference comes from a rubidium clock source in Nançay. The NenuFAR beamformer produces 768 simultaneous beamlets with 195 kHz bandwidth at a time resolution of 5.12 µs at any desired frequency in the band 10-85 MHz and in any direction, within the analog beam of the phased MAs. This setup allows us in particular to cover the full bandwidth (75 MHz) with two digital beams (768×195 kHz = 2×75 MHz), or to observe using four digital beams with a bandwidth of 37.5 MHz each (768×195 kHz = 4×37.5 MHz). By default, these four digital beams use all MAs; alternatively, different sub-arrays (with different sets of MAs) can be used for each of these four beams, allowing for a fly's eye mode with up to four totally independent directions on the sky, corresponding to four analogical beams. For all cases, the full polarisation information of the signal is recorded in the form of four Stokes parameters. More detail on NenuFAR are given in Zarka et al. (2020). Flux and polarisation calibration are in progress and will be described in the full instrument paper to come (Zarka et al., in prep.).

De-dispersion
The main effect of the ionised interstellar medium is the dispersion of the radio waves, where the dispersive delay is proportional to the line-of-sight integrated electron content, the so-called dispersion measure (DM). Being also proportional to ν −2 , where ν is the observing frequency, the dispersive delay has a strong impact on low frequency radio observations and requires dedicated instrumentation to process the signal and detect pulsed emission. The lowest frequency channel that can be coherently de-dispersed is directly related to the block size of the instrumentation (≈ 10 seconds) and the intra-channel dispersive delay corresponding to the DM of the observed pulsar, i.e. 12 MHz for DM∼10 pc cm −3 , 25 MHz for DM∼100 pc cm −3 , 53 MHz for DM∼1000 pc cm −3 .
There are two ways to correct for dispersive delays (see e.g. Lorimer & Kramer 2004;Hankins & Rickett 1975), as schematically shown in Fig. 1. Incoherent de-dispersion (in red) corrects for inter-channel delays only and, thus, requires a large number of narrow channels to minimise intrachannel smearing. It is generally applied off-line, once the signal has been recorded. Coherent de-dispersion (in blue) is designed to correct for intra-channels delay, managing the waveform signal by changing its phase in the Fourier domain. This process makes it possible to have wide channels while keeping a high time resolution. This method does not introduce dispersion smearing into the data in contrast with incoherent de-dispersion. It is clear that coherent de- shows a strongly dispersed pulsation from a pulsar in black, and the incoherently and coherently dedispersed signal in red and blue, respectively. The channels are separated by grey dotted lines (195 kHz wide for both LOFAR and NenuFAR). Panels (B) and (C) represent the de-dispersed and frequency integrated profiles using both methods. Amplitude scales are arbitrary but the areas below the red and blue curves should be similar.
dispersion is necessary to preserve the integrity of the pulsar pulse profile, especially at low frequencies. To correct for the intra-channel dispersion at low frequency (10-85 MHz) this method requires long Fourier transforms and a coherent dispersion code dimensioned for the corresponding buffer. As example the intra-channel dispersion delay of a signal with a DM value of 10 pc cm −3 is about 0.03 s at 85 MHz, 0.13 sec at 50 MHz, and 16.21 s at 10 MHz.

The UnDySPuTeD LUPPI
The new real-time code LUPPI (Low frequency Ultimate Pulsar Processing Instrumentation, Bondonneau et al. 2018) operating on NenuFAR was adapted from NUPPI (Nançay Ultimate Pulsar Processing Instrument, Desvignes 2014), the Graphical Processing Unit (GPU) cluster for the Nançay Radio Telescope (NRT) observations since August 2011. NUPPI is designed to handle dispersive delays in 4 MHz channels in L-or S-bands for a total instantaneous bandwidth of 512 MHz. The NUPPI software was derived from GUPPI (Green Bank Ultimate Pulsar Processing Instrument, Demorest 2014 andDuPlain et al. 2008).
LUPPI runs on two composite machines named UnDyS-PuTeD (Unified Dynamic Spectrum Pulsar and Time Domain receiver), each hosting two Intel Xeon E5-2620 CPUs (Central Processing Units) and two Nvidia GeForce GTX 1080 GPUs.
By using together the two UnDySPuTeD machines, LUPPI has at its disposal a total of 32 CPU cores, 512 GB of error-correcting code RAM memory and 32 BG of GPU memory. Because of the observing frequencies below 100 MHz, LUPPI has to deal with additional difficulties compared to GUPPI and NUPPI.  Fig. 2. Operational scheme of the LUPPI pipeline including its main elements: Initialisation (A) of the observation parameters, the data processing cycle organized in threads around a shared memory and communicating through semaphores and post-processing of the observation with production of "quicklooks" plots. The network thread (B) catches the packets send by the beamformer, the computing thread (C) processes the data in the Fourier domain on GPUs and the write thread (D) exports data products to disks.
The huge intra-channel dispersion imposes extremely long Fourier transforms (typically more than 20 seconds, i.e. 200 times longer than for similar computations at 1.4 GHz) and leads to memory management difficulties (both RAM and GPU) in the NenuFAR framework. These difficulties have been mitigated by modifying the size of the buffers as well as the variables that contain it in the source code. These modifications enable the correction of intra-channel dispersive delays up to 10.737 s.
LUPPI captures the 768 beamlets produced by the NenuFAR beamformer (described in Zarka et al. 2018;Dumez-Viou et al. 2019) as a series of network packets. The two complex-voltage values (one per polarisation) at each time step are used to compute the four Stokes parameters of the signal.
The operational scheme of LUPPI is shown in Fig. 2. The pipeline is divided into four real-time elements and a post-processing analysis: A) The initialisation phase prepares the parameters of the observation (DM, channels, source name, etc.) and allocates a total of 128 GB 3 of shared memory. We use shared 3 buffers size: 4 (4×37.5 MHz) × 2 (input and output buffer) × 4 (number of blocks per buffer) × 2 (dual linear polarisation) × 2 (2 byte complex) × 2 22 (number of points per Fast Fourier Transform (FFT) × 192 (number of channels) 8×12.88 GB. The buffer size is rounded up to 8×32 GB to be compatible with a future 16-bit mode. memory to implement ring buffers 4 in order to allow threads to exchange data and process in parallel. Furthermore it is during this phase that the coherent dispersion chirps 5 are calculated.
B) The network thread catches UDP (User Datagram Protocol) packets sent by the beamformer and reads the packet index to look for any missing packets and takes it into consideration to preserve the linearity of the time series. The packets are then stored in the shared memory and a semaphore is used to unleash the computation thread.
C) The computation thread is dedicated to feed GPUs with data and get the results back. At this stage, the data blocks are sent to four GPUs (one processing unit per band of 37.5 MHz).
First, a fixed overlap 6 is added to absorb the intrachannel dispersive delay. For a given GPU hardware, the value of the overlap results as a compromise between computing efficiency (which is better when the overlap is small) and the possibility to reach low frequencies (for a given DM value, a larger overlap allows to process data at a lower frequency). With our GPUs being limited to 8 BG of memory, we chose an overlap of 50% of the total block size, i.e. 10.737 seconds. With this, we can de-disperse frequency channels of frequency >12 MHz for a DM of ∼10 pc cm −3 , >25 MHz for a DM of ∼100 pc cm −3 , >53 MHz for a DM of ∼1000 pc cm −3 .
Then, the complex voltages are transformed to the Fourier domain using the NVIDIA CUDA Fast Fourier Transform library (cuFFT https://docs.nvidia.com/cuda/cufft) in order to multiply the signal by the transfer function corresponding to the coherent dispersion chirp. One of the unique features of this instrumentation is the availability of a Faraday coherent de-rotation allowing to correct for the rotation of the linear polarisation within each channel (Bondonneau et al., in prep).
Furthermore, the data are transformed back to the time domain with a inverse Fast Fourier Transform (iFFT) and folded at the period of the observed pulsar or integrated in time depending on the observing mode. Finally, the result is pushed into the output buffer and a semaphore is activated to launch the write thread.
D) The write thread writes the data product out from the output buffer into a RAID 0 (Redundant Array of Independent Disks using data striping) storage unit.
The last step is the post-processing to verify that the observation output products are present on the server and to create diagnostic plots ("quicklooks", see §2.4.4). The quicklook plots are sent by email to the observer within minutes. If data files are missing on the server a warning email is sent to the observer. The source code of LUPPI is available on the Github 7 .

Currently available observing modes
NenuFAR and LUPPI currently offer several observing modes for pulsars. We describe the main ones below.

Folded mode
Time series are coherently de-dispersed in every 195 kHz channel and folded at the apparent period of a pulsar using an up-to-date ephemeris. The default sub-integration time is 10.737 s. Data are stored in PSRFITS 8 format (Hotan et al. 2004) and a quicklook plotting program displays the main features of the observation after removal of radio frequency interference (RFI). This is the main observing mode. It is used to characterise the pulse profile and polarisation as a function of frequency, analyse the spectral energy distribution, and study long-term variations, e.g. related to DM fluctuations or intrinsic behaviour in the pulsar magnetosphere.

Single-pulse mode
In this mode a time series is de-dispersed and integrated in time by a chosen factor (128 by default) to enable the analysis of pulse-to-pulse variations. Data are stored in PSRFITS format and a quicklook program displays the main features of the observation after RFI removal (using the rfifind software from PRESTO 9 , Ransom 2001). This mode is used to observe bright pulsars that show interesting features in their emission beam, such as drifting subpulses, mode changing, or anomalously intense pulses. It will also be used to search for pulsed emission when the DM is higher than 100 pc cm −3 and the scattering is important. Moreover, single pulse data can easily be converted into folded data.

Waveform mode
This mode consists in writing the raw complex-voltage data to disk, bypassing the computing thread on Fig. 2. It keeps the full time, frequency and polarisation information from the original raw signal (thus a very large data rate). This mode can be used to perform non-standard analyses, such as scattered giant pulses searches, a census of the Globular Clusters pulsar population, or interstellar and interbinary scintillation studies. Additionally, waveform data can be converted into folded or single pulse data by coherent dedispersion and integration. The two UnDySPuTeD machines also provide LUPPI with additional, independent modes described below:

Dynamic spectrum mode
The 195 kHz wide beamlets sent by the beamformer are GPU-processed to provide finer frequency resolution and time integration. The frequency resolution can be chosen between 0.10 and 12.20 kHz (corresponding to FFT lengths from 2048 to 16 samples). The time resolution (corresponding to a number of integrated FFTs restricted to be a power of 2, and being at least 4) can be chosen between 0.30 and 84.0 ms (depending on the frequency resolution). An apodisation function (windowing) can be applied among the few ones proposed (such as Hann, Hamming, or others), the default being Hamming. The resulting dynamic spectra are recorded in a very simple binary format. This dynamic spectrum mode is used for solar or Jupiter observations, exoplanet searches, and in the context of a pulsar blind survey or any radio transient search. This mode can be used to observe pulsars without de-dispersion, for instance the bright pulsar B1919+21 (Fig.3). In addition, it supports very fine channels (e.g. 1.5 kHz) enabling searches over wider DM ranges without coherent de-dispersion (see a detection example in Fig. 4).

Compressed UDP waveform mode
For a small number of projects we use a second waveform mode, in which UDP packets are directly dumped to disk when they arrive on the network socket. In this mode, we use real-time lossless compression with zstd 10 , which allows us to reduce the data volume by a factor of ∼2-3. Files written in this mode do not contain any header information. This mode is used, for example, for low-frequency Very Long Baseline Interferometry (VLBI) with international LOFAR stations on baselines of ∼1000 km.

The post-processing pipeline
In order to analyse observations we have developed two post-processing pipelines for folded and single-pulse observations. The role of these pipelines is to ensure that the transfer between the real-time processing machine and the data storage is successful and to redirect the data into the different branches of processing according to their mode of observation: waveforms, folded or single-pulses.

Waveform mode
The waveform mode corresponding to the raw data requires large processing and storage resources (540 BG per hour for a bandwidth of 37.5 MHz), consequently it is not handled by the quicklook pipeline even if the transfer remains automatic.

Folded mode
The post-processing pipeline for the folded observations is based on Nenuplot.py 11 . Based on the PSRCHIVE 12 library. This program is fed with the folded FITS files to perform the RFI mitigation and get out a quicklook in PNG format as described in the next section. Finally, the quicklook is automatically sent by e-mail to the observer.

Single-pulse mode
Single pulse observations are converted from 32 bits to 8 bits and normalised in time and frequency domains with a fits file converter 13 . Furthermore, the blocksize of a PSR-FITS file is reduced with RescaleTime.py 14 in order to 10 https://github.com/facebook/zstd 11 Nenuplot.py: a completely automatised pipeline for folded observations. Based on the PSRCHIVE library (Hotan et al. 2004)   Top-right panels: reduced χ 2 relative toṖ for the best P0 and relative to P0 for the bestṖ .
allow more accurate RFI mitigation with rfifind. Then, the PSRFITS file is analysed with PRESTO tools using a fine DM grid adapted to low frequencies to search for intense pulses in time and DM. PRESTO produces a quicklook analog to that shown in Fig. 6 which is sent by e-mail to the observer. Moreover, we use dspsr 15 to transform single-pulse 8-bit PSRFITS files into folded and de-dispersed PSRFITS files allowing us to use Nenuplot.py and produce a quicklook (see in §2.4.4) and send it to the observer with the quicklook of PRESTO.

The quicklook
The creation of a quicklook is a crucial part of the observation because it makes it possible to monitor the operation of the pipeline as well as the quality of the observation with respect to, e.g., RFI. Its aim is to allow quick validation of an observation and to reduce the response time in case of problem. This python software is based on the PSRCHIVE library and it is designed to be a complete pipeline compatible with all radio telescopes generating folded PSRFITS files. A quicklook of an observation of PSR B0329+54 is shown in Fig. 5. The quicklook plot presents the observation data in all its dimensions: time, phase, frequency, and polarisation in various panels. Panel 1 shows the pulse profile with the total intensity (in black), the linear polarisation (in red), and the circular polarisation (in blue). The profile allows us to determine if the observed source is detected or not. Panel 2 shows the metadata of the observation including the source name, DM, rotation measure (RM), P0, observation duration, bandwidth, signal-to-noise ratio (SNR), percentage of RFI, elevation of the source, etc. Panel 3 shows the temporal variation of the frequency-integrated profile. This plot is useful to highlight changes in pulsar period with respect to the reference ephemeris used for folding the observation or a time variation of the pulse profile. Panel 4 is the time-integrated bandpass (i.e. the spectral response of the telescope) for the polarisation parameters XX and YY (in blue and red, respectively). Channels flagged because of the presence of RFI are highlighted in grey in this panel. Panel 5 shows the time-integrated profile as a function of frequency for Stokes I (total intensity). This plot is important because at low frequency we are particularly sensitive to changes in DM. Panels 6, 7, and 8 show the frequency resolved profiles for the Stokes parameters Q, U, and V. The colour corresponds to the sign of the parameter (red for positive values and blue for negative). These plots are helpful to characterise the RM used to correct for the Faraday rotation between channels. In the example of Fig. 5 there is a frequency dependent structure in the parameters Q and U. This is the signature of an imperfectly corrected Faraday effect caused by ionospheric variations (i.e. the total RM is different from the value in the reference ephemeris used during the observation). Panel 9 shows the dynamic spectrum of the amplitude of the pulsation (recorded in the "ON-pulse" region), which can be used to quantify mitigation and data quality. Panel 10 is the dynamic spectrum of the bandpass (spectral response from the "OFF-pulse" regions) that can be used to assess the quality of the observation. The scintillation of the source is visible on panel 9, while panel 10 shows a 6-min periodic drop of the bandpass (gain) that result from the analog pointing of the instrument and can be corrected in post-processing.

Data archive
Currently, NenuFAR pulsar data are transferred to the Nancay Data Center (CDN) for temporary storage and distribution within the community. At the CDN, the data are stored on a RAID 6 storage cluster. Pulsar data recorded in folded mode ( §2.3.1), single-pulse mode ( §2.3.2), and dynamic spectrum ( §2.3.4) mode will be kept at the CDN, whereas waveform data ( §2.3.3) and data recorded in UDP waveform mode ( §2.3.5) are usually only kept temporarily and are removed after processing. A dedicated NenuFAR Data Centre is currently being commissioned, it will have open access user interfaces to data products.

Commissioning and first early science results
Early Science phase of the instrument will last until the end of 2021. During this time we will be conducting a set of observational programs with NenuFAR, under the umbrella of the "Pulsars Key Science Project" (ES03). Up to now, all our observations were made with a limited version of the full instrument, i.e. the 56-MA first step configuration that started operating with a stable setup in February 2019 (Zarka et al. 2020). A 80-MA configuration will be operational by the end of 2020.

First NenuFAR pulsar catalogue
A first targeted survey has already allowed us to observe 650 pulsars with DM < 100 pc cm −3 above −20 • of declination, of which about 26% were detected (Bondonneau et al., in prep). Where the non-detections are mainly due to the amplitude of the scattering at theses frequencies. These observations are an important complement to recent studies of the pulsar population observable below 85 MHz such as  and .
This first NenuFAR catalogue, combined with the first measurements of flux density and polarisation, is now used as a basis for a set of follow-up programs, e.g. the study of pulse-to-pulse variations with time and frequency, the effects of scintillation due to the interstellar medium, the characterisation of the pulsar population and the eclipses of certain systems. A set of 41 bright pulsars is tracked monthly to study the long-term variations of their average profiles, to characterise the DM evolution, and to study their spectral energy distribution in detail. Finally, a blind survey of the northern sky visible from Nançay started in Summer 2020 and will provide a more complete census of the local pulsar population at frequencies below 100 MHz.

Details of the emission profile of PSR B1919+21
The high sensitivity of NenuFAR and its stable gain across the band make it ideally suited to study average profiles in detail. Figure 7 shows the evolution of the observed profile for the historical pulsar B1919+21, discovered by Jocelyn Bell in 1967 (Hewish et al. 1968). Using NenuFAR, this pulsar has been detected with a high SNR (> 2000) down to 20 MHz. Obtaining good-quality data at such low frequencies is important for modelling profile evolution and constraining emission heights. About three dozen sources have been observed with such frequency resolution as in Fig. 7 in integrations of only one hour (Bondonneau 2019;).

Crab pulsar giant pulses
It has long been known that the Crab pulsar emits so-called giant pulses (GPs), very brief bursts of radio emission with intensity many times exceeding the intensity of an average pulse (e.g. Howard et al. 1968;Cordes et al. 2004;Hankins et al. 2016;Rudnitskii et al. 2016). However, at such low frequencies these pulses are highly scattered (at 60 MHz, the scattering timecale is 5 times the pulsar's spin period) and consequently difficult to observe in folded mode. There have been a few single-pulses campaigns to observe giant pulses at low frequencies: for example, Popov et al. (2006) reported on the detection of 45 pulses in 12 hours with the Ukrainian UTR-2 radiotelescope below 30 MHz; Eftekhari et al. (2016) announced the detection of 1400 pulses in 76 hours of observation with the LWA1. With NenuFAR we detected 600 pulses in only 13.8 minutes at a frequency of 68.7 MHz (Fig. 6). This detection rate and the telescope's sensitivity are very promising for future NenuFAR blind searches for new pulsars in the Northern hemisphere. A strong individual pulse of this observation as a function of the observed frequency is shown in Fig. 8. At 52 MHz the scattering tail of this pulse is roughly 20 times longer than the pulsar spin period.

Drifting sub-pulses in PSR B0809+74
Drifting patterns are a common phenomenon in pulsars (e.g. Weltevrede et al. 2006Weltevrede et al. , 2008Basu et al. 2019  of them show systematic variations of their sub-structures in an organised manner. Among the ∼70 known pulsars that exhibit this behaviour, an interesting prototype of this phenomenon is the bright and extensively studied pulsar B0809+74. It was observed down to 16 MHz with LOFAR LBA (Stappers et al. 2011;Hassall et al. 2013), and detected in the band 18-27 MHz with UTR-2 (Zakharenko et al. 2013). The drifting sub-structures and the anoma- lously intensive pulses of PSR B0809+74 (AIPs, Ulyanov et al. 2006) are resolved in great details by NenuFAR ( Fig. 9) and are visible all along its frequency domain, from 85 down to 10 MHz. Because of the large fractional bandwidth and good sensitivity, NenuFAR data are very well suited to study frequency-dependent phase delay of drifting sub-pulses. This delay can provide novel constraints on the location of emission regions in the magnetosphere (Bilous et al. 2019;Maan 2019).

Detection of 10 MSPs
At frequencies below 100 MHz, only five millisecond pulsars (MSPs) have been detected prior to this work. The first detection of PSR J2145−0750 was published in Dowell et al. (2013), PSRs J0030+0451 and J0034−0534 in Stovall et al. (2015), PSR J1400−1431 in Swiggum et al. (2017) and finally, PSR J0437−4715 in Bhat et al. (2018). The first three MSPs were also reported to be detected with the LOFAR LBA core by Kondratiev et al. (2016). Thanks to NenuFAR, we already doubled this list (Bondonneau et al., in prep). As an example, in Fig. 10 we show the observed average pulse profile of PSR J1022+1001, obtained during the commissioning phase of the instrument. It was detected with a DM of 10.25356(3) pc cm −3 , which is compatible with observations at higher frequencies and the SNR of 46 corresponding to a mean flux density of 28 ± 14 mJy (following the same method as the one used to calculate the mean flux density limit in section 3.9).

DM variations from pulsar timing
The study of time variations in DM for different lines of sight enables studies of the ionised interstellar medium and its turbulence (see, e.g. Phillips & Wolszczan 1991;Backer et al. 1993;Donner et al. 2019). This is particularly interesting at low frequencies, where the dispersive time delay is amplified by its dependence on ν −2 . Such observations demand a large instantaneous fractional bandwidth, high sensitivity, and good clock stability, all of which are provided by NenuFAR. DM variations are typically studied via pulsar timing and the calculation of pulse times-of-arrival (TOAs). First tests have shown that for NenuFAR, the typical uncertainties of TOAs are ∼5-8 times smaller than for the LOFAR station FR606 (used in LBA mode, i.e. in the same frequency range). This translates to an increased precision for the study of DM variations. Fig. 11 shows DM vari- ations for PSR B1919+21 as measured by NenuFAR and LOFAR FR606 obtained with a 2D template fit (a phasefrequency template extending over 118 channels, built from the longest observation of PSR B1919+21 by NenuFAR, showing a SNR > 6000). The data from both telescopes are comparable, but the NenuFAR data provide a much higher precision (of the order of 10 −5 pc cm −3 ). As the template used for pulsar timing is derived from the data itself, the procedure only gives relative DM values (sufficient for the analysis of DM variations) rather than absolute DM values (Donner et al. 2019). Series of DMs are well aligned especially if we look at the simultaneous observation in the left panel in Fig. 11.

Interstellar Scintillation Studies
By recording data in waveform mode, an additional filterbank step can be applied to increase frequency resolution (similar to the dynamic spectrum recording mode) before phase-averaging the observations (as in the folded mode). This hybrid combination of observing modes allows dynamic spectra to be constructed of the pulsed SNR with much higher frequency resolution than the standard folded mode allows; and with more sensitivity than the standard dynamic spectrum mode allows. This is particularly useful for the study of interstellar and interbinary scintillation (Rickett 1990), since the size of scintles is strongly dependent on the observing frequency, causing typical scintillation bandwidths to be much smaller than 195 kHz at NenuFAR frequencies. Figure 12 shows an example of a dynamic spectrum that was obtained with this method, on PSR B0809+74. The frequency resolution was set to 1.5 kHz and the data were folded into 10-second integrations. The strong frequency scaling of the scintle size is apparent from the two 3-MHzwide segments of data shown here: in the 50-53 MHz range the scintillation bandwidth (i.e. the size of scintles in frequency) is 3.9 ± 0.7 kHz; in the 80-83 MHz band it is 23 ± 1 kHz. The scintillation timescale (the size of the scintles in time) scales less strongly: from 3.79 ± 0.07 min in the 50-53 MHz subband to 6.1 ± 0.3 min in the 80-83 MHz band. This high sensitivity and the fact that even a narrow frequency range can achieve numerous scintles, allows thorough and self-consistent statistical analyses of scintillation scaling laws and enables precise monitoring of scintillation variations and their possible relation to other interstellar effects.
3.8. Telescope characterisation 3.9. Mean flux density limit The detection limit of NenuFAR can be estimated using the radiometer equation (Lorimer & Kramer 2004) with the same set of parameters as in  for the LOFAR single station FR606, but using the collecting area and the instrument temperature of NenuFAR for 56 MA. The result is shown in Fig. 13 which compares the minimal flux density detectable by NenuFAR for 56 MA (blue curve) with the corresponding limit of FR606 (red curve) and the pulsars detected by FR606. In this Figure the minimum detectable flux density is calculated for a pointing outside the Galactic plane (gl =0 • , gb =90 • , T sky = 2350 K), assuming 4 hours of integration time, on a 75 MHz bandwidth, a SNR of 5 and a pulse duty cycle of 7%. The pulsars displayed are detected by the FR606 station with duty cycles greater than 7% ). There is a large jump in sensitivity between the in- Fig. 12. Dynamic spectra of PSR B0809+74, in two 3 MHz wide bands, of a one-hour NenuFAR observation. The colour scale indicates the pulse SNR ranging from blue (no signal) to yellow (high-SNR), which is heavily modulated due to diffraction in the interstellar medium. The high-SNR "islands" are commonly referred to as scintles and provide information on the turbulent interstellar plasma.  struments, which is mainly due to a factor of 11 in the collecting area. This jump will allow for a large number of new detections during the NenuFAR targeted survey.

RFI environment at NenuFAR site
Like all radio telescopes NenuFAR is subject to RFI. The post-processing pipeline for the folded mode (Section 2.3.1) is equipped with an automatic RFI cleaner (Section 2.4.2). Based on the observations of the last two years, Fig. 14   72 MHz. A few frequencies are permanently occupied by narrow-band RFI (e.g. 35.5 MHz, 72.5 MHz, and 78 MHz). Above 87.5 MHz, the band is saturated by frequency modulation broadcasting. Below 30 MHz Fig. 14 shows broadband interference (10-25 MHz) which is much brighter and extends to higher frequencies in the daytime than during nighttime. At daytime, RFI heavily affects the observations of pulsars below 20 MHz; at nighttime, observations down to 10 MHz remain possible.

Summary
We have provided a full description of the NenuFAR pulsar instrumentation, including the four main modes used for pulsar observations. We have shown a number of first results from our observations, part of which were obtained already during the commissioning phase. These early science results are very promising and demonstrate the potential of NenuFAR for pulsar science; more detailed follow-up studies are currently being performed within the NenuFAR "Pulsars Key Science Project". We expect significant contributions from NenuFAR, in particular in the fields of pulsar emission mechanisms and the study of the interstellar medium (ISM) and ionosphere.
With its high sensitivity at low frequencies, NenuFAR is also a milestone on the route to the Square Kilometer Array (SKA), hence its pathfinder label 16 . It will help refining the SKA science goals, particularly concerning emission mechanism and propagation effects in the context of highly dispersed and scattered signals. As such, it will contribute to maintain a useful bridge of pulsar studies across the next decade and support the community of pulsar astronomers focusing on emission processes, populations, and ISM characterisation. These science goals are complementary to the ones carried out at L-and S-bands for high precision pulsar timing (gravitational wave detection and tests of General Relativity in the context of the International Pulsar Timing Array collaboration, IPTA). When the SKA will be operational, joint studies between the SKA and NenuFAR will be performed in the two hemispheres in complementary frequency ranges, providing the potential for additional discoveries.