The LOFAR Tied-Array All-Sky Survey: Timing of 21 pulsars including the first binary pulsar discovered with LOFAR

Here we describe the observational strategy and setup, as well as the timing analysis of the isolated pulsars discovered by LOTAAS. Whenever a pulsar is detected and confirmed by a follow-up observation with LOTAAS (see Sanidas et al. 2019 for a description of the follow-up observation strategy), it is added to the monthly LOTAAS timing campaign. In the campaign, the pulsars were observed with all 24 HBA stations of the LOFAR core. The lengths of individual observations of each pulsar depend on the brightness of the pulsar and are listed in Table 1. The dual-polarization complex voltage data were recorded 400 sub-bands of 195.3 kHz each, centred at an observing frequency of 149.8 MHz. The data were then processed with the automated Pulsar Pipeline (PulP; Stappers et al., 2011; Kondratiev et al., 2016), which coherently dedisperse and fold the data using dspsr (van Straten & Bailes, 2011) and psrchive (Hotan et al., 2004), respectively to produce an archive from each observation. The archive files are data cubes folded with 1024 phase bins across one pulse period, five seconds time sub-integrations and 400 sub-bands across the observing band of 78.1 MHz. The setup is more sensitive than the LOTAAS survey setup, where only six HBA stations of the LOFAR Superterp are used.

The pulsars were also all observed with the Lovell Telescope at Jodrell Bank, at a central observing frequency of 1532 MHz. Initially, each pulsar was observed 4$-$5 times over a span of about 10 days. Each observation was between 30 min and 1 hour long. If the pulsar was detected, a lower cadence observing campaign was continued, with each source typically observed once every two weeks, depending on the availability of the telescope. The observation length is shown in Table 1, and was chosen based on the brightness in the initial detection. The observations were conducted at a central observing frequency of 1532 MHz, with 768 channels across a bandiwidth of 384 MHz. They were processed with the digital filterbank backend (DFB), where the data were incoherently dedispersed to the DM of the pulsar and folded at the best known period. The folded archive files have sub-integrations that are 10 s long and 1024 phase bins. If a pulsar was not detected with the Lovell Telescope at 1532 MHz, a dense observing campaign with the full LOFAR core was conducted. Four observations of 15 minutes each were conducted over a period of 10 days, with at least two of the observations over consecutive days. The data recording and processing followed the same setup used for the monthly timing campaign.

The folded archive files obtained from the observations were analysed with psrchive. The pulse times-of-arrival (TOAs) of each pulsar were measured by cross-correlating the integrated pulse profiles of each observation with an analytical, noise free template. The templates were created by fitting between one and four von Mises functions¹¹1The von Mises function is a close approximation to a wrapped Gaussian function. See Evans et al. (2000) for a detailed description of the function., depending on the complexity of the pulse profiles (Kramer et al., 1994), with paas of psrchive. An initial template was created by fitting the pulse profile of the observation with the highest signal-to-noise ratio (S/N) across the full LOFAR bandwidth to provide an interim timing solution. The LOFAR observations were then averaged into two frequency channels, with central frequencies of 128 MHz and 167 MHz in order to measure the $\mathrm{DM}$. Different templates were made for each observing frequency, taking into account possible profile evolution across the bandwidth (Hassall et al., 2012), fitting for the pulse profiles averaged across the observing span. The templates were then aligned at the half-maximum of the most prominent component, in order to take into account possible effects of scattering that might shift the peak of the profile at low frequencies to slightly later phase. For the isolated pulsars, only the average DM across the observing span is measured, without modelling the DM variation between individual observations. A phase-connected timing solution was then determined by modeling the TOAs with tempo2 (Edwards et al., 2006; Hobbs et al., 2006). Any offsets between the TOAs obtained with LOFAR and the Lovell telescope were corrected by fitting a jump between the two sets of data. The Solar System ephemeris DE405 was used to model the motion of the bodies in the Solar System in order to convert the topocentric TOAs into barycentric TOAs.

As there are uncertainties in the initial spin periods of the pulsars from follow-up observations, there is an ambiguity in the number of pulsar rotations between successive monthly observations. In order to resolve the ambiguity and phase-connect the monthly TOAs, the TOAs from the dense campaign were used to model the rotation of the pulsars over short timescales. We were able to model the TOAs of all the pulsars with the TOAs from the dense observations with just a refinement of the spin period of the pulsars, producing coherent timing solutions that take into account all pulsar rotations within the short observing span. This allowed us to extrapolate the solutions to longer timescales without introducing ambiguities, which we then used to model the positions and spin period derivative of the pulsars.

The large apparent period derivative measured from the confirmation observation of PSR J1658$+$3630 strongly suggested that the pulsar is in a binary system. In order to model the binary orbit, more frequent observations of the pulsar are required, on top of the regular monthly observations. An initial dense observing campaign was first conducted to model the binary parameters of the system. The observing strategy is presented in Table 2 as Cam. 1.Only the Stokes I data were recorded for these observations, with a central frequency of 149.8 MHz and 6400 12.2 kHz channels. The data were incoherently dedispersed to the DM of the pulsar and folded with the best known ephemeris at the time of the observations, producing archive files with 10 s sub-integrations and 256 phase bins.

The binary orbit of PSR J1658+3630 was first determined by modelling the changes in the apparent spin period of the pulsar. The apparent spin periods at each epoch were measured by refolding the data at slightly different spin periods from the best known period and finding the spin period that produces the highest S/N using pdmp of psrchive. We then used the fitorbit tool of psrtime²²2http://www.jb.man.ac.uk/pulsar/observing/progs/psrtime.html to model the apparent spin period of the pulsar at different epochs. We found that the variation of the apparent spin period is roughly sinusoidal, suggesting the pulsar is in an almost circular orbit and can be tentatively modelled assuming zero eccentricity. An initial timing solution was produced, modelling the orbital period, projected semi-major axis and the epoch of ascending node of the system and the spin period of the pulsar. This incoherent timing solution was then used to further refine the rotational, astrometric and binary properties of the pulsar using the timing procedure in Section 2.1 to obtain a coherent timing solution.

The incoherent timing solution of PSR J1658$+$3630 suggests a possibility of measurement of Shapiro delay of the pulsar signal due to its companion. Hence, a second dense campaign for PSR J1658$+$3630 was conducted to provide a large coverage of orbital phases, specifically at orbital phase 0.25, where the binary companion is in front of the pulsar. The observing strategy is described in Table 2 as Cam. 2.These data were recorded, processed and analysed using the same method as the monthly timing observations, as described in Section 2.1.

PSR J1658$+$3630 was also observed with five of the international LOFAR stations in Germany, namely stations in Unterweilenbach (telescope identifier DE602), Tautenburg (DE603), Bornim (DE604), Jülich (DE605) and Norderstedt (DE609), operated by the German LOng Wavelength consortium. The observing strategy and the individual stations involved are presented in Table 2. The observations were conducted at a central frequency of 153.8 MHz and with a bandwidth of 71.5 MHz across 366 sub-bands. The data from the stations in Unterweilenbach, Tautenburg and Jülich were recorded on machines at the Max-Planck-Institut für Radioastronomie in Bonn, while the data from stations in Bornim and Norderstedt were record at machines in Jülich Supercomputing Centre. They were recorded using the LOFAR und MPIfR Pulsare (LuMP4) software³³3https://github.com/AHorneffer/lump-lofar-und-mpifr-pulsare as channelized complex voltages and then coherently dedispersed to the DM of the pulsar and folded using the best ephemeris of the pulsar available in 2017 July, producing archive files with ten-second sub-integrations and 1024 phase bins. A summary of the different observing strategies on PSR J1658$+$3630 is shown in Table 2.

The timing observations conducted with both LOFAR and the Lovell Telescope were used to study the pulse profile evolution and the spectral properties across different frequencies. Additionally, the pulsars were observed at a central observing frequency of 334 MHz with the Lovell Telescope to provide a larger frequency coverage in these studies. Each pulsar was observed once for 30 minutes. These data have a bandwidth of 64 MHz across 512 frequency channels and were processed with the DFB. The folded archive files have sub-integrations of 10 s and 512 phase bins.

The pulse profiles of the pulsars at 149 and 1532 MHz were obtained by adding the integrated pulse profiles obtained from each observation together with psradd of psrchive, after aligning the profiles from different observations using the timing solutions presented in Table 3. The pulse profiles at 334 MHz were obtained from the single observations made. We measured the pulse widths and the duty cycles of the pulsars at 10 and 50 per cent of the maximum, W₁₀, $\delta_{10}$, W₅₀, $\delta_{10}$, respectively, using the following method. First, noise free templates of the pulse profiles were generated by fitting a number of von Mises functions, each with a different height, width and position using paas. The templates formed are the summation of individual von Mises functions that best represents the overall shapes of the profiles. The observed off-pulse region was then used to generate a noise distribution. One thousand simulated pulses were generated by creating new noisy profiles combining the noise free template with the simulated noise. These were then fitted with the noise free templates, allowing height, width and position of each individual function to vary. The widths of the remodelled templates at 10 and 50 per cent of the maximum from the 1000 trials were then used to determine the mean width and error. The widths are defined as the outermost components of the templates that are at 10 per cent and above, and 50 per cent and above the maximum.

The LOFAR observations were flux-calibrated using the method described in Kondratiev et al. (2016) in order to measure the flux densities of the pulsars. The calibrated data of individual observations were then averaged into two frequency bands to measure the flux densities of the pulsars at 128 and 167 MHz respectively. The average flux density of each pulsar was then computed from these observations. We found that the offset between the initial pointing position of the observations of most of the pulsars and the measured position obtained through pulsar timing is of the order of one arcmin. The pointing offset was corrected for several of the pulsars studied here in later observations. We attempted to correct the measured flux densities due to loss in sensitivity as a result from the offsets from the early observations. However, we found that pulsars with initial pointing positions that are offset from the timing position by less than 1 arcmin did not show any improvement in S/N after updating the pointing position, while those with offset of more than 1 arcmin showed S/N improvement that is less than expected. We suspect that this could be due to two issues. First, the ionosphere is known to induce a jitter in the position of the pulsar in the sky up to 1 arcmin around the actual position of the source. We also suspect that the shape of the tied-array beam formed cannot be modelled as a 2D-Gaussian due to the complexity of the beam shape and the beam not being fully coherent. Hence, we decided that, for pulsars with an initial pointing offset less than one arcmin, flux density measurements for all observations were used without applying any correction, while for those with larger offsets, only the observations after the correction of pointing position were used to measure the flux density of the pulsars. The error in the measured flux densities of each frequency band of a single observation is conservatively estimated to be 50 per cent (Bilous et al., 2016; Kondratiev et al., 2016).

The flux densities of pulsars detected at 334 and 1532 MHz were estimated using the radiometer equation (Lorimer & Kramer, 2005). The receiver temperature T_rec of the 334 MHz receiver is 50 K and the sky temperature T_sky in the direction to the relevant pulsar was estimated by extrapolating the T_sky at 408 MHz (Haslam et al., 1982) with a spectral index of $-2.55$ (Lawson et al., 1987; Reich & Reich, 1988). The Gain, G, of the receiver is 1 K Jy^-1 and the bandwidth was 64 MHz. As for the 1532 MHz observations, T${}_{\mathrm{rec}}=25\,\mathrm{K}$, G$=1\,\mathrm{K}\,\mathrm{Jy}^{-1}$ and the bandwidth was 384 MHz. We estimated the average RFI fraction of each observation at 1532 MHz to be 20 per cent of the bandwidth and 5 per cent of the observing length, while at 334 MHz, the fraction is 20 per cent of the bandwidth and 20 per cent of the observing length. Upper limits for the non-detections were estimated using the radiometer equation with a threshold S/N of 10 and estimated pulse width based on the measured W₅₀ at 149 MHz. The uncertainty on the flux density of a single observation is estimated to be 20 per cent.

A notable exception to this sample of pulsars is PSR J0210$+$5845. It is located in a more sparsely populated region of the diagram, with a characteristic age of $\tau_{\text{c}}\approx 200$ kyr and surface dipole magnetic field of $\text{B}\approx 1.1\times 10^{13}$ G. The pulsar also shows large timing residuals, which could be timing noise related to its small characteristic age (Hobbs et al., 2010). Additionally, the large timing residual suggest that high order frequency derivatives are potentially measurable in the future with a larger observing span of a few years.

PSR J0107$+$1322 shows three distinct peaks in its profile at 149 MHz, which are well modelled with three components that increase in intensity from the leading to the trailing peaks. However, the profiles at 334 and 1532 MHz only show two distinct peaks; a weak leading peak and a strong trailing peak. The profiles are better modelled with two components rather than with three components. It is possible that the profiles at higher frequencies has 3 separate components which are not visible due to low S/N of the observations. While the increase in W₁₀ from 149 to 334 MHz suggest a profile evolution contrary to the RFM model, higher S/N observations of PSR J0107$+$1322 at higher frequencies are required to confirm the relationship, as the low S/N observation at 334 MHz could also introduce an extra uncertainty to the determination of the off-pulse baseline to determine W₁₀.

PSR J1643$+$1338 shows two distinct peaks in its profiles at all three observing frequencies, all of which are well modelled with just a single component per peak. Both W₁₀ and W₅₀ of the pulsar are shown to increase with increasing frequency. We also measured the separation between two components (Table 6). We found that the separation between the components increases at higher frequencies as well, suggesting that the pulsar does not conform to the RFM model, where the component separations are expected to decrease with increasing frequencies. While the behaviour of PSR J1643$+$1338 is unusual, the anti-RFM profile evolution has been seen in several other pulsars before (Hassall et al., 2012; Chen & Wang, 2014; Noutsos et al., 2015; Pilia et al., 2016).

The pulse profiles of PSR J1749$+$5952 show two distinct peaks at 149 and 334 MHz, while at 1532 MHz, the trailing peak seems to merge with the leading peak. The profiles at 334 and 1532 MHz can be modelled relatively well with just two components. However, the profile at 149 MHz requires three components; the two main components describing each peak and a bridge component that is fitted across the profile. We measured two distinct sets of W₅₀ values of the profiles at both 149 and 334 MHz. We found that this is due to the trailing peak of both profiles having roughly half the intensity of the leading peak. The smaller values in Table 5 correspond to the width of the main peak at half maximum while the larger values correspond to the width of the whole profile.

The measurements of the separation of the components corresponding to the two peaks seen in the profile of J1749$+$5952, shown in Table 6 suggest that the separation between components decreases with increasing frequency, in agreement with RFM. However, the values of W₅₀ that correspond to the width of the main peak at 149 and 334 MHz suggest that the width of the main peak increases with increasing frequency. The increase in width of the main peak at higher frequency could be affected by the exact beam pattern of the pulsar, and pulse longitude dependent spectral index effects.

The integrated pulse profiles of PSR J1810$+$0705 show significant frequency evolution. At 149 MHz, the profile shows two distinct peaks of about equal intensity and at 1532 MHz, the intensity of the leading peak relative to the trailing peak is reduced. The pulsar shows two distinct peaks at 334 MHz as well, however, the low S/N of the observation does not allow us to get a good estimate of the relative intensity between the peaks. The W₁₀ of the pulsar also shows a large decrease with increasing frequency, suggesting that the separation between the two components decreases with increasing frequency as expected. Due to the complex structure of the profile of PSR J1810$+$0705 at 149 MHz, we are unable to model the two peaks as distinct components to measure the component separation. The templates that were used for pulsar timing requires several components for each peak to produce an adequate model of the pulse profiles.

PSR J1916$+$3224 shows a strong main peak and a trailing, wider component at all three observing frequencies. The evolution of the pulse profile is found to follow the RFM model, where the W₅₀ of the profiles, which correspond to the width of the main peaks, decreases as observing frequency increases. While the measured W₁₀ suggests the profile width is constant between 334 and 1532 MHz, the separation between the two main components (Table 6) decreases with increasing observing frequency, in agreement with RFM model.

The emission of PSR J1657$+$3304 shows several interesting properties. Firstly, it shows long-term flux density variation over the observation span. We measured its flux density at 149 MHz across 15 observations over 347 days, shown in Fig. 6, and it varies by a factor of 10 over approximately 300 days before staying constant at low flux density (1.5 mJy) after MJD 58300. The magnitude of the variation is much larger than any other pulsars discussed in this work. We attempted to identify any potential period derivative change that correlates with the flux density variation, as the changes in flux density of the pulsar could be due to long-term mode changing that is often accompanied with changes in the period derivative of the pulsar (Kramer et al., 2006; Lyne et al., 2010). We did not detect any noticeable change in period derivative.

We also looked into the possibility that the flux density variation is due to interstellar scintillation in the strong regime. We found that the variation is unlikely to be due to diffractive interstellar scintillation as the expected scintillation bandwidth at 149 MHz is less than 1 kHz according to the NE2001 model (Cordes & Lazio, 2002), much smaller than the bandwidth of the LOFAR timing observations of 78 MHz. The long-term flux density variations of pulsars are often due to refractive interstellar scintillation (RISS; Rickett et al., 1984). However, we found that this is unlikely for PSR J1657$+$3304, as the expected timescale of RISS, $t_{\textrm{RISS}}$ has a frequency dependence of $t_{\textrm{RISS}}\propto\nu^{-2.2}$, which means that the timescale in the changes in flux density of the pulsar between the top (167 MHz) and bottom half (128 MHz) of the LOFAR band is expected to be different by a factor of 1.8. We found that the changes in flux density of the pulsar to be consistent across the LOFAR bandwidth. Hence, it is more likely that the flux density variation of PSR J1657$+$3304 is intrinsic to the pulsar.

Individual observations of the pulsar also showed nulling over duration between several pulses to a few minutes (Fig. 7). We estimated the nulling fraction of each observation of PSR J1657$+$3304 following the method of Wang et al. (2007), where the pulse energy distribution of the pulsar is compared with the off-pulse energy distribution. We also tested whether the obtained nulling fractions are the same for both high and low flux density observations, by increasing the noise level of the high flux density observations to have the same S/N as the low flux density observations and then calculating the nulling fraction. We find that the measured nulling fraction changes with S/N as it is more difficult to separate the pulses from the nulls. Hence we only estimated the nulling fraction from the four observations with flux densities above 10 mJy. The average nulling fraction is found to be 34 per cent. While the presence of nulling could account for some of the flux density variation seen between observations, we found that the flux density of individual sub-integrations also varied between observations.

In addition to nulling, we observed short-term mode changing in some observations, where the pulse profile switches between the commonly occurring double-peaked structure to only the leading peak being present (Fig. 8). The duration of the occurrence of each instance of the less common mode ranged between 1-5 minutes, similar to the duration of the nulls. For example, the less common mode occurred between third and sixth minute of the observation in MJD 58202, as shown in Fig. 7. Most of the mode changing occurred between nulls, with several exceptions such as at the beginning of the observation in MJD 58118 where the pulsar changes mode before a period of nulling.

PSR J1657$+$3304 was observed with the Lovell Telescope at both 334 and 1532 MHz. The observation at 334 MHz is strongly affected by RFI, and the pulsar is not detected in four 30-minute observations at 1532 MHz. However, the observations at 1532 MHz coincide with the period when PSR J1657$+$3304 showed low flux density at LOFAR observing frequencies. This gives an upper limit to the flux density of the pulsar at 1532 MHz of 0.06 mJy when it is in the low flux density state.

The low eccentricity of the orbit of PSR J1658$+$3630 led us to use the ELL1 binary timing model (Lange et al., 2001) to model the TOAs. The parameters of the timing model are shown in Table 8 and the timing residuals as a function of both time and binary phase are shown in Fig. 9. The pulsar has a spin period of 33 ms and spin period derivative of 1.16$\times 10^{-19}$ s s^-1, placing it in the location of the $P$-$\dot{P}$ diagram (Fig. 2) that is populated by the intermediate-mass binary pulsars (IMBP), suggesting that PSR J1658$+$3630 is part of the IMBP population. The average DM of the pulsar is 3.044 pc cm^-3, indicating DM-distances of either 0.22 (Yao et al., 2017) or 0.49 kpc (Cordes & Lazio, 2002), making PSR J1658$+$3630 one of the closest-known pulsar binaries. The orbital period is 3.016 days and the companion has a minimum mass of 0.87 $M_{\odot}$. The minimum mass of the companion is sufficiently high that the measurement of Shapiro delay is possible if the system has a high inclination angle. However, we found that the pulsar showed temporal variation in DM, which has a large effect on the TOAs at LOFAR frequencies, we are unable to measure any Shapiro delay, and that the timing residuals in terms of binary orbital phase (Fig. 9) does not show signs of Shapiro delay (See Demorest et al. 2010 for the effects of Shapiro delay on the timing residuals).

We modelled the DM variation of PSR J1658$+$3630 across the observing span by measuring the average DM value over spans of 30 days. The measured DM value in each epoch is then compared to the average across the observing span and is shown in Fig. 10. We found that for most of the observing span, the DM value fluctuates between $-0.0002$ and $+0.0001$ pc cm^-3 around the average DM value of 3.0439 pc cm^-3, with the observations since MJD 58280 showing larger DM increase of between $+0.0002$ and $+0.0005$ pc cm^-3 more than average. The timing solution took into consideration the measured DM variation. The other pulsars studied in this work are likely to show DM variation as well. However, due to their lower timing precision compared to PSR J1658$+$3630, the variations are too subtle to be measured.

PSR J1658+3630 is coincident (within $1.2\arcsec$) with SDSS J165826.50$+$363031.1, a $r=22.1$, $g-r=0.56$ star in the Sloan Digital Sky Survey (SDSS, York et al., 2000; Albareti et al., 2017). This optical counterpart is also seen in images from the Panoramic Sky Survey Telescope and Rapid Response System (Pan-STARRS, Chambers et al., 2016; Flewelling et al., 2016). Using white dwarf cooling models by Bergeron et al. (2011)⁴⁴4http://www.astro.umontreal.ca/~bergeron/CoolingModels/, as described in Holberg & Bergeron (2006); Kowalski & Saumon (2006); Tremblay et al. (2011) and Bergeron et al. (2011), we find that the SDSS $u-g$, $g-r$ and $r-i$ colors of the counterpart are consistent with those of a 0.9 M_⊙, $T_{\mathrm{eff}}\sim 5250$ K white dwarf. Given this, and the low stellar density in SDSS (of order 9 stars per square arcminute with $r<23$), we conclude that the counterpart is the white dwarf companion of PSR J1658+3630. Given it’s high mass of $M_{\mathrm{WD}}>0.87$ M_⊙, it is likely to be a Carbon-Oxygen or a Oxygen-Neon-Magnesium type white dwarf. We found that the position of the pulsar obtained from the timing solution is different from both the SDSS observation taken in 2000 and the Pan-STARRS observation taken in 2014, suggesting that the system has moved across the sky over the years. With this assumption, we derived a proper motion of $+0\aas@@fstack{\prime\prime}015\pm 0\aas@@fstack{\prime\prime}008$ yr^-1 in right ascension and $-0\aas@@fstack{\prime\prime}065\pm 0\aas@@fstack{\prime\prime}008$ yr^-1 in declination. The overall proper motion indicates the pulsar has a transverse velocity of between 69-154 km s^-1 based on the DM-distances obtained. The proper motion estimated here is large enough that a sinusoidal variation with increasing amplitude will be seen in the timing data if it is not accounted for. As we only have one and a half years of timing data, we were not able to refine the measured proper motion through pulsar timing yet. Hence, we fixed the proper motion of the pulsar at the values measured while refining the other properties.

While the timing model is able to broadly describe the rotational and orbital properties of the pulsar, we note that there remains a significant timing residual in the TOAs and a poor $\chi_{\mathrm{red}}^{2}$ value. There are several possible contributions to the large timing residuals. First, the proper motion measured has a significant uncertainty of $0\aas@@fstack{\prime\prime}008$ yr^-1 in both right ascension and declination. The proper motion can be measured to a higher precision with pulsar timing, but this will require a longer observing span on the order of several years. We also noticed that there is an unmodelled extra delay on the timing residuals at around MJD 58165. This is likely a short duration increase in the DM of the pulsar on a time-scale shorter than the 30 days used to model the DM variation of the pulsar.

The average pulse profile integrated over the observing span shows frequency-dependent profile evolution. Fig. 11 shows the integrated pulse profiles at four different LOFAR sub-bands and at 334 MHz, overlaid with the model that best describes the pulse profiles. The profiles in the LOFAR bandwidth are described by three von Mises components, corresponding to the leading bump, the main peak and a small trailing component that increases in relative intensity to the main peak as the observing frequency increases. The profile at 334 MHz has a much lower S/N and was fitted with just 2 components. We are unable to identify if the trailing component is present at 334 MHz. We also measured the width of the pulses at 50% and 10% of the maximum, shown in Table 9.

The pulsar is also observed to undergo diffractive scintillation in some of the observations, where the scintillation bandwidth appeared to be smaller than the bandwidth of a single sub-band of 195 kHz. This is much smaller than the expected bandwidth of 1.4 MHz predicted by the NE2001 model. While the detailed study of the scintillation properties of the pulsar will be presented in a separate paper, we note that the presence of diffractive scintillation might affect the overall pulse profile shape of individual observations due to the variation in flux densities in different parts of the bandwidth at different observing epochs. This will subsequently affect the timing precision of these observations due to profile shape changes relative to the template used. The DM variation can affect the overall shape of the profile as well, as different parts of the bandwidth may not be aligned at the fiducial point of the template.

We also calculated the average flux densities of PSR J1658$+$3630 at the same four LOFAR sub-bands and the flux density from the 334 MHz observation in which the pulsar is detected. The average flux densities at the LOFAR frequencies were calculated from 41 observations obtained from the second dense campaign and are shown in Table 9. The diffractive scintillation and possibly refractive scintillation will result in variation in the flux densities. We model the spectral index of the pulsar using a single power law, shown in Fig. 12, and measured $\alpha=-2.5\pm 0.7$. The pulsar is therefore a steep spectrum source, similar to the other MSPs discovered by LOFAR (Bassa et al., 2017; Pleunis et al., 2017). The spectral index predicts a radio luminosity at 1400 MHz, $L_{1400}$ of between 0.0015 and 0.007 mJy kpc², depending on the distance obtained by the electron density model used (Cordes & Lazio, 2002; Yao et al., 2017). This is similar to another recently discovered nearby, low-luminosity MSP J2322$-$2650 ($L_{1400}=0.008$ mJy kpc^-2; Spiewak et al., 2018), strengthening the argument that there could be a large population of low-luminosity Galactic MSPs that can only be detected at very low radio frequencies.