1. Introduction
The radio spectra of pulsars offer unique insights into the nature of the mysterious pulsar emission mechanism, and for this reason have been the subject of extensive study for many decades (e.g. Sieber Reference Sieber1973; Malofeev & Malov Reference Malofeev and Malov1980; Izvekova et al. Reference Izvekova, Kuzmin, Malofeev and Shitov1981; Lorimer et al. Reference Lorimer, Yates, Lyne and Gould1995; Maron et al. Reference Maron, Kijak, Kramer and Wielebinski2000). Furthermore, studies of pulsar spectra are useful resources for planning surveys of the Galactic pulsar population with the Square Kilometre Array (SKA). In particular, spectral modelling can be used to inform pulsar population studies, which play a key role in estimating the yields of future surveys (e.g. Xue et al. Reference Xue2017; Keane et al. Reference Keane2015). The SKA is expected to discover thousands of new pulsars and millisecond pulsars (MSPs), which will be critically important in developing pulsar timing arrays (PTAs) for the detection of ultra-low-frequency gravitational waves (e.g. Manchester et al. Reference Manchester2013). PTAs are one of a broad range of high-profile applications that have led to pulsars being ranked as a headline science for the SKA (e.g. Janssen et al. Reference Janssen2015).
It is well known that the majority of pulsars exhibit steep power-law spectra ( $S_\nu\propto\nu^{\alpha}$ ; where $\alpha$ is the spectral index and $S_\nu$ is the flux density at frequency $\nu$ ) with a mean spectral index of $\langle{\alpha}\rangle={-1.60\pm 0.03}$ (Jankowski et al. Reference Jankowski2018). However, the radio spectra of some pulsars have been observed to reach a peak in flux density at low radio frequencies (a spectral ‘turnover’), which is most commonly attributed to either synchrotron self-absorption (Sieber Reference Sieber1973) or thermal free-free absorption (Malov Reference Malov1979). Interestingly, studies of the MSP population have found that they do not turn over like long-period pulsars (e.g. Kramer et al. Reference Kramer1999; Kuzmin & Losovsky Reference Kuzmin and Losovsky2001; Toscano et al. Reference Toscano, Bailes, Manchester and Sandhu1998), but rather continue as power-laws to the lowest observable frequencies (e.g. Erickson & Mahoney Reference Erickson and Mahoney1985). Only a handful of MSPs have shown hints of turning over (e.g. Kuniyoshi et al. Reference Kuniyoshi2015; Dowell et al. Reference Dowell2013), all of which are predicted to peak well below ${100}\,\mathrm{MHz}$ , which suggests that turnovers in MSPs may occur at much lower frequencies than long-period pulsars. Deviations from a simple power-law (such as spectral flattening and low-frequency turnover) are referred to as spectral features. Jankowski et al. (Reference Jankowski2018, hereafter JvSK+18) performed a study of the spectral properties of 441 pulsars and found that 21% of pulsars show spectral features, increasing to 44% for pulsars with good low-frequency coverage (i.e. more than two data points below ${600}\,\mathrm{MHz}$ ). Building upon this work, Johnston et al. (Reference Johnston2021) classified an additional 44 pulsars, 20% of which show spectral features.
Unfortunately, the vast majority of catalogued pulsars lack reliable flux density measurements below ${400}\,\mathrm{MHz}$ , which puts poor constraints on spectral models at these frequencies. A recent version of the Australia Telescope National Facility (ATNF) pulsar catalogueFootnote a (version 1.67, released in 2022 March; Manchester et al. Reference Manchester, Hobbs, Teoh and Hobbs2005) shows that ${\sim}70\%$ of the over 3300 known pulsars do not have flux density measurements available below ${400}\,\mathrm{MHz}$ . This is despite the emergence of the latest generation of low-frequency aperture array telescopes such as the Murchison Widefield Array (MWA; Tingay et al. Reference Tingay2013; Wayth et al. Reference Wayth2018), the LOw-Frequency ARray (LOFAR; van Haarlem et al. Reference van Haarlem2013; Stappers et al. Reference Stappers2011), and the Long Wavelength Array (LWA; Taylor et al. Reference Taylor2012; Ellingson et al. Reference Ellingson2013), which have led to the publication of flux densities for a considerable number of known pulsars at low frequencies (e.g. Murphy et al. Reference Murphy2017; Xue et al. Reference Xue2017; Bondonneau et al. Reference Bondonneau2020; Bilous et al. Reference Bilous2020; Bilous et al. Reference Bilous2016; Stovall et al. Reference Stovall2015).
In general, pulsar flux density measurements can be obtained in two ways: (1) measurements of pulsed emission from beamformed detections, which often have inaccurate absolute flux calibrations and therefore limit constraints on spectral models; and (2) measurements from interferometric continuum images, which are less sensitive, but provided that data are recorded for a large field-of-view (e.g. the MWA Voltage Capture System (VCS); Tremblay et al. Reference Tremblay2015), can be robustly calibrated by virtue of the ability to perform in-field calibration against hundreds of other calibrator sources. Continuum measurements may also suffer from source confusion (i.e. blending with other sources in the field), depending on the spatial resolution. The vast majority of flux density measurements are made from pulsed emission, however the advent of widefield low-frequency aperture arrays has enabled continuum flux density measurements to be made for a small number of pulsars (e.g. Murphy et al. Reference Murphy2017).
The low-frequency component of the SKA (SKA-Low), to be built at the Murchison Radio-astronomy Observatory (MRO) in Western Australia, will operate between $50{-}350\,\mathrm{MHz}$ , and will be capable of providing substantially improved constraints on pulsar emission properties at low frequencies, including the nature of the spectra. Currently operational at the MRO are two second-generation prototype SKA-Low stations: the Aperture Array Verification System 2 (AAVS2; van Es et al. Reference van Es, Marshall, Spyromilio and Usuda2020; Macario et al. Reference Macario2022) and the Engineering Development Array 2 (EDA2; Wayth et al. Reference Wayth2022). Aside from providing early insights into the capabilities of SKA-Low technologies, these stations are being used to prepare for science to be conducted with the SKA-Low. The stations are particularly useful tools for performing pulsar monitoring due to the real-time beamforming capability, which yields more manageable data volumes and enables faster data processing turn-around times in comparison to the MWA-VCS.
In this work, we demonstrate the early pulsar detection capabilities of the AAVS2 and EDA2 by performing low-frequency observations of a selection of bright southern-sky pulsars, measuring flux densities for the detected pulsars, and modelling their radio spectra. In Section 2, we describe the facilities, target selection, and observations. In Section 3, we outline the flux density measurement and calibration methodology, as well as the methods used to model and classify the pulsar spectra in a robust and unbiased way. In Section 4, we present the results of a shallow all-sky census and multi-frequency follow-up observations and report flux density measurements and spectral classifications for the detected pulsars. Finally, we summarise and give conclusions in Section 5.
2. Observations
2.1 SKA-low prototype stations
2.1.1 The AAVS2
The AAVS2 is a low-frequency aperture array that comprises 256 log-periodic dual-polarisation SKALA-4.1 antennas distributed pseudo-randomly over a circular ground plane with a station diameter of ${\sim}{38}\,\mathrm{m}$ . A description of the AAVS2 is provided in van Es et al. (Reference van Es, Marshall, Spyromilio and Usuda2020), and a characterisation of its performance can be found in Macario et al. (Reference Macario2022). The antenna arrangement replicates the EDA2 and the previous generation stations—the AAVS1 (Benthem et al. Reference Benthem2021) and the EDA1 (Wayth et al. Reference Wayth2017)—to enable direct comparisons to be made, however the diameter of the AAVS2 was upscaled to accommodate the larger antennas. The 50– ${350}\,\mathrm{MHz}$ observing band is split into coarse channels which are separated by ${\approx}{0.781}\,\mathrm{MHz}$ and oversampled up to a bandwidth of ${\approx}{0.926}\,\mathrm{MHz}$ . For pulsar observations, the station forms a phased-array beam by coherently summing the complex voltages of each array element. The station beam can be electronically steered to an arbitrary pointing direction by applying phase delays to each of the antenna signals using digital beamforming implemented in the firmware of Tile Processing Units (TPMs; Naldi et al. Reference Naldi2017). The beam response varies as a function of both frequency and pointing direction, with the sensitivity generally being highest at the zenith. The initial system, which was employed for the work presented in the paper, allowed data capturing from a single coarse channel at a maximum time resolution of ${\approx}{1.08}\,\unicode{x03BC} \mathrm{s}$ . Since then, the system has been upgraded for a much larger bandwidth ( ${>}32$ coarse channels), the testing and commissioning of which is in progress.
2.1.2 The EDA2
Developed as a comparator to the AAVS2, the EDA2 uses the same analogue and digital signal chain. However, the EDA2 makes use of MWA-style dual-polarisation bowtie dipoles with modified low-noise amplifiers (LNAs) which extend the sensitivity down to ${\sim}{50}\,\mathrm{MHz}$ to replicate the frequency range of the AAVS2. Details of the design, calibration, and performance of the EDA2 are given in Wayth et al. (Reference Wayth2022). Its predecessor, the EDA1, was used in Bhat et al. (Reference Bhat2018) to observe two MSPs at low frequencies between ${\sim}50{-}300\,\mathrm{MHz}$ , giving an early demonstration of pulsar science with an SKA-Low prototype station.
2.2 Beamforming: SKA-Low stations vs the MWA
The MWA is a low-frequency aperture array telescope which serves as a precursor to the SKA-Low and has been operational at the MRO since 2013 (Tingay et al. Reference Tingay2013). The system is capable of producing high-time- and high-frequency-resolution observations, which has enabled a wide range of meaningful low-frequency pulsar science to be performed, from studies of emission properties, multipath propagation, and pulsar timing, to the first pulsar discoveries with the MWA (e.g. McSweeney et al. Reference McSweeney, Bhat, Tremblay, Deshpande and Ord2017; Meyers et al. Reference Meyers2018; Bhat et al. Reference Bhat2018; Kaur et al. Reference Kaur2019; Kirsten et al. Reference Kirsten2019; Swainston et al. Reference Swainston2021; McSweeney et al. Reference McSweeney2022). The SKA-Low stations have ${\sim}20\%$ of the MWA’s effective collecting area (as indicated by simulations at ${160}\,\mathrm{MHz}$ ), and ${\sim}3\%$ of the bandwidth, giving them ${<}4\%$ of the sensitivity. However, the MWA performs offline incoherent and tied-array beamforming on tile voltages recorded at $28\,\mathrm{TB\,h}^{-1}$ using the high-time-resolution VCS backend (Tremblay et al. Reference Tremblay2015). In comparison, the SKA-Low stations record beamformed data at only $12\,\mathrm{GB\,h}^{-1}$ (for a single coarse-channel). This makes the data easier to process due to the beamforming having already been performed in real-time, effectively enabling the stations to operate like single-dish telescopes. The more manageable data volumes and faster data processing turn-around times in comparison to the MWA-VCS make the stations useful tools for performing data-intensive science such as high-cadence pulsar monitoring.
2.3 Selection of target pulsars
The AAVS2 and EDA2 are relatively new instruments that can be employed for conducting a shallow all-sky census (i.e. sensitive primarily to bright pulsars, due to the limited sensitivity of the initial system). For this work, we compiled a modest sample of bright southern-sky pulsars with version 1.64 of the ATNF pulsar catalogue using the following constraints: (1) a declination limit of $\delta<+30^\circ$ (J2000), the same as the MWA; (2) a mean flux density at ${400}\,\mathrm{MHz}$ of $S_{400}>40\,\mathrm{mJy}$ or a mean flux density at 1400 MHz of $S_{1400}>5\,\mathrm{mJy}$ ; and (3) a dispersion measure (DM) cut-off of $450\,\mathrm{cm}^{-3}\,\mathrm{pc}$ . Since the vast majority of catalogued pulsars do not have flux densities reported below ${400}\,\mathrm{MHz}$ , the flux density constraints were calculated by extrapolating higher-frequency flux densities down to ${150}\,\mathrm{MHz}$ using a nominal spectral index of $\alpha=-1.6$ . Based on the expected sensitivities of the SKA-Low stations, we chose $200\,\mathrm{mJy}$ to be the minimum predicted flux density at ${150}\,\mathrm{MHz}$ . Finally, the DM cut-off was chosen to eliminate pulsars with very high DMs and to limit the number of pulsars to 100. The pulsars observed in this work are shown on a Galactic sky map in Figure 1.
The SKA-Low stations are capable of observing nearly three octaves of frequency range, which makes them useful instruments for studying pulsar emission properties over a large and uncommonly observed part of the radio spectrum. To this end, nine pulsars were selected as targets for follow-up observations at multiple frequencies. To select these targets, the successful detections from the initial all-sky census were ranked in order of signal-to-noise (S/N) ratio to give preference to the brightest pulsars, and daytime transiting pulsars were excluded. One of the selected pulsars (J1820–0427) did not yield any further detections from follow-up observations.
2.4 Observations
The pulsars selected for follow-up observations were observed at 18 different centre frequencies from 70.3–351.6 MHz. The observations for each pulsar were conducted over the course of a single night, one at a time, starting at the lowest frequency channel and incrementing to the next highest selected channel with each consecutive observation. The follow-up observations were 15 min in duration (half that of the 30 min census observations), separated by 2 min intervals.Footnote b PSRs J0835–4510, J0953+0755, and J1645–0317 were simultaneously observed with the EDA2 and the AAVS2, and the remaining pulsars were observed solely with the EDA2 while the AAVS2 was taken offline for a hardware upgrade. Further observations were taken of PSR J0835–4510 to assess the accuracy of the beam sensitivity simulations across a range of elevations. Data collected with the stations were saved onto dedicated on-site servers and later transferred to the Pawsey Supercomputing CentreFootnote c where they were stored and processed.
Despite the limited recording bandwidths ( ${\sim}{1}\,\mathrm{MHz}$ ), the SKA-Low stations have demonstrated the capabilities to make high-S/N detections of several bright pulsars. Figure 2 shows an EDA2 observation of PSR J0953+0755: a short-period ( $P\approx{0.25}\,\mathrm{s}$ ), low-DM ( ${\approx}2.96\,\mathrm{cm}^{-3}\,\mathrm{pc}$ ) pulsar which is known to exhibit strong variability in its observed intensity due to interstellar scintillation (Bell et al. Reference Bell2016; Sokolowski et al. Reference Sokolowski2021). The observation was taken while the pulsar was fortuitously in a brightened state, which resulted in a peak flux density of over 200 Jy. Single pulses have been detected in this observation, and their further investigation will be left for a future publication. The frequency versus phase plot shows negligible dispersion sweep ( $\Delta t\approx{2}\,\mathrm{ms}$ ) due to the small observing bandwidth and low DM, but it can still be discerned in the figure.
Pulsars are very highly polarised radio sources, and polarimetric analysis of their emission can provide useful information about their beam geometry and the intervening interstellar medium (ISM). In principle, the SKA-Low stations are capable of producing full-Stokes pulse profiles which can be used to assess the polarimetric calibration of the stations (e.g. Xue et al. Reference Xue2019). The pulse profile in Figure 2 shows the linear and circular polarisations and the corresponding linear polarisation position angle curve of the bright observation of PSR J0953+0755. Although the observation has not been derotated, this pulsar has a small rotation measure ( $\mathrm{RM}\approx{-0.66}\,\mathrm{rad\,m^{-2}}$ ), and over the small recording bandwidth the effect of Faraday rotation is minimal. McSweeney et al. (Reference McSweeney2020) published a polarisation pulse profile of this pulsar using MWA data at ${152}\,\mathrm{MHz}$ which shows first-order agreement with the EDA2 observation.
3. Data processing and analysis
3.1 Integrated pulse profiles
Raw voltages were saved into binary files containing 5 min of data each, which were subsequently merged and prepended with a header containing the pulsar ephemeris, period, DM, and other observing metadata. The data were then processed into archive files using the dspsr software packageFootnote d (van Straten & Bailes Reference van Straten and Bailes2011). The observations were coherently dedispersed, and a 256-point convolving filterbank was run to fine-channelise the coarse channel. The time series was folded into 64 phase bins for initial detections and 256 phase bins for flux measurements to increase the number of off-pulse bins for calibration. The psrchive data analysis libraryFootnote e (Hotan, van Straten, & Manchester Reference Hotan, van Straten and Manchester2004) was used to average the data in time and frequency and construct the pulse profiles shown in Figure 3.
3.2 Calculation of mean flux densities
To calculate a mean flux density S, the $N_\mathrm{on}$ on-pulse phase bins were added together and the result was normalised by the pulsar period P, i.e.
where $I_i$ is the Stokes I value of the i-th on-pulse phase bin, $\Delta t$ is the bin width, and $\Lambda$ is a calibration constant. The on-pulse was estimated using the formula
where $n_\mathrm{r,l}$ are the bin indexes of the right and left on-pulse boundaries, $n_\mathrm{max}$ is the index of the bin with the maximum Stokes I value across the profile, $N_{3\sigma}$ is the number of phase bins with a Stokes I value more than three standard deviations from the mean, k is a scaling constant, and N is the total number of phase bins. The $N_{3\sigma}$ value gives an estimate for the number of on-pulse phase bins, which is most accurate for profiles with short duty cycles, but typically underestimates the on-pulse for profiles with larger duty cycles. The value of k (typically 2 or 3) was selected to ensure that all on-pulse estimates were generous enough to pick up any residual emission. Each on-pulse estimate was checked visually to ensure sensibility.
The pulse profile of PSR J0835–4510 is heavily scattered at low frequencies, with the scattering tail dominating the off-pulse region at frequencies below around ${300}\,\mathrm{MHz}$ . This makes it difficult to calibrate reliably using measurements of the off-pulse noise. For the lowest-frequency detections, we processed the integrated pulse profile with a DM of zero to determine an upper-limit of the background noise level. We then used this estimate to calibrate the profile and calculate the flux density. Since these measurements are likely underestimates, they were excluded from the spectral fits. This method works under the condition that the dispersion sweep exceeds the pulse period. For a single coarse-channel, this condition is only met for J0835–4510 below ${185}\,\mathrm{MHz}$ .
3.3 Calibration procedure
Calibration of the flux densities required estimations of the system-equivalent flux density (SEFD; i.e. a measure of the system sensitivity) in the pointing direction and at the frequency of each observation. The SEFD is defined as
where $T_\mathrm{sys}$ is the system temperature and G is the gain. At low frequencies, $T_\mathrm{sys}$ is dominated by the sky background, which is a strong function of frequency ( $\propto\nu^{-2.55}$ ; Lawson et al. Reference Lawson, Mayer, Osborne and Parkinson1987). The gain is also a function of frequency due to its dependence on the effective collecting area of the array elements. To calculate the SEFD, the beam response was simulated using python codeFootnote f detailed in Sokolowski et al. (Reference Sokolowski2022). The simulation code finds $T_\mathrm{sky}$ and G, calculates the SEFD in the X and Y polarisations, and combines them in quadrature according to
where $\mathrm{SEFD}_I$ is the SEFD in the Stokes I polarisation. Equation (4) assumes an unpolarised target source and is only strictly valid in the cardinal planes ( $\unicode{x03D5} =0^\circ,90^\circ$ ; where $\unicode{x03D5}$ is the azimuth angle), which leads to errors at low elevations (Sutinjo et al. Reference Sutinjo2021). The radiometer equation was used to calculate the standard deviation of the background noise expected in images of the sky,
where $n_p=2$ is the number of polarisations, $t_\mathrm{int}$ is the integration time, and $\Delta\nu$ is the bandwidth. We have assumed $\sigma_\mathrm{sim}$ to be equal to the off-pulse noise of beamformed observations. By comparing $\sigma_\mathrm{sim}$ to the standard deviation of the off-pulse noise in the uncalibrated pulse profile ( $\sigma_\mathrm{off}$ ), a flux density calibration constant $\Lambda$ was calculated:
The results from the sensitivity simulation code were compared against a full-electromagnetic simulation of the EDA2 station beam in the commercial simulation package feko.Footnote g It was verified that the simulation code consistently calculated SEFDs within ${\sim}30\%$ of the feko simulation. Due to limitations in time and resources, the flux densities could not all be calibrated against feko simulations. Instead, the relative error in the SEFD estimates from the simulation code was accounted for in the uncertainty calculations.
A set of observations were performed at a range of elevations from ${\sim}20{-}70^\circ$ for a bright pulsar (J0835–4510) to test the flux calibration far from the zenith. The error was negligible at elevations greater than ${\sim}40^\circ$ , so no elevation-dependent corrections were implemented.
3.4 Compiling spectral data
Flux density measurements are often reported in pulsar discovery papers and census papers, however compiling these resources can often prove time consuming. In the interest of streamlining this process for future work, we have developed an open-source database of published spectral data called pulsar spectra,Footnote h which will be described in full in a future publication (Swainston et al. submitted). A complete list of publications included in the database can be found in the documentation.Footnote i The current version of the database (v1.4) comprises the majority of pulsar flux density measurements available at radio frequencies, with particular emphasis on compiling low-frequency measurements, where spectra often exhibit flattening or curvature. For the pulsars analysed in this work, our database is nearly up-to-date.
3.5 Modelling the pulsar spectra
Robust modelling of data combined from the literature is difficult due to the measurements being made with different telescopes and calibration procedures that have varying levels of reliability. In many cases, this introduces mismatches in systematic errors, which can make modelling the data complicated. Our parameter fitting method is based on the work of JvSK+18, who addressed these issues using methods from robust regression and information theory. The spectral modelling code used in this work is included in version 1.4 of pulsar spectra. Further details regarding the modelling method can be found in JvSK+18 and the references within, and a full description of our implementation of the method will be given in Swainston et al. (submitted). A high-level summary of the modelling method is given below.
We used a Gaussian likelihood function, which quantifies the probability of a data set given a set of model parameters. The likelihood was modified using the Huber loss function (Huber Reference Huber1964), which penalises outlier data points by reducing their contribution to the model fit. A robust cost function was then derived by taking the negative logarithm of the likelihood function. The optimal parameters were found by minimising the cost function, which is equivalent to maximising the likelihood function.
The cost function was minimised using the migrad algorithm from the minuit2 minimisation library, which was integrated into pulsar spectra via the python bindings in iminuit Footnote j (Dembinski et al. Reference Dembinski2020). migrad converges towards a local minimum using a combination of Newton steps and gradient descents and is called a set number of times before returning the best-fitting model parameters (James & Roos Reference James and Roos1975). The quality of the minimisation was verified by ensuring that the best-fitting parameters satisfied a convergence criterion, defined in terms of a specified tolerance value. The uncertainties were computed using hesse, an error calculator which computes the covariance matrix for the best-fitting parameters and determines the 1 $\sigma$ uncertainties as the square-root of the diagonal elements.
The five analytical models identified in JvSK+18 were used as a representative sample of the various morphologies reported in the literature. These models are simple power-law, broken power-law, log-parabolic spectrum, power-law with a low-frequency turnover, and power-law with a high-frequency hard cut-off. Following JvSK+18, the parameter defining the smoothness of the low-frequency turnover was restricted to a value of $0<\beta\leq 2.1$ . In each spectral model, the frequency was scaled by a reference frequency which was calculated as the geometric mean of the minimum and maximum frequencies in the spectrum.
In contrast with JvSK+18, we allowed the spectral index of the high-frequency cut-off model to be a fit parameter. This modification to the application of the model is in accordance with the original work (Kontorovich & Flanchik Reference Kontorovich and Flanchik2013) and enables the model to be fitted to a larger sample of spectra.
The best-fitting model was determined using the Akaike information criterion (AIC), which measures how much information the model retains about the data without overfitting. The model which resulted in the lowest AIC was selected as the best-fitting model. Furthermore, the probability that the selected best-fitting model is the true best-fitting model amongst those tested, $p_\mathrm{best}$ , was calculated as the inverse sum of the likelihoods of each model relative to the best-fitting model.
4. Results and discussion
4.1 Pulsar detections
Detections were made of 22 pulsars using data collected with the SKA-Low stations. The non-detection of the remaining 78 pulsars is due to a range of factors, including overestimated flux densities and interstellar scattering. Our initial estimates for the flux densities at ${150}\,\mathrm{MHz}$ (i.e. $S_{150}$ ) were made using a rudimentary simple power-law assumption (see Section 2.3) with single flux density measurements from the ATNF pulsar catalogue (i.e. $S_{400}$ or $S_{1400}$ ). In light of the development of pulsar spectra, we performed spectral fits for all of the non-detected pulsars to estimate how many of the non-detections can be attributed to overestimated flux densities due to spectral flattening or turnover. Only 29 ( ${\sim}40\%$ ) of non-detected pulsars are above the sensitivity limit of the stations (i.e. $S_{150}>200\,\mathrm{mJy}$ , conservatively). Out of these, 14 have high DMs ( ${\geq}100\,\mathrm{cm}^{-3}\,\mathrm{pc}$ ). Although DM smearing cannot be a factor (the observations were coherently de-dispersed), high-DM pulsars are much more likely to experience scattering by the ISM (e.g. Bhat et al. Reference Bhat, Cordes, Camilo, Nice and Lorimer2004), which can reduce the S/N and hence the detection rate of these pulsars. The remaining 15 bright ( $S_{150}>200\,\mathrm{mJy}$ ) low-DM ( ${<}100\,\mathrm{cm}^{-3}\,\mathrm{pc}$ ) non-detected pulsars include the Crab pulsar (PSR J0534+2200), which was likely not detected due to being heavily scattered at low frequencies (e.g. Meyers et al. Reference Meyers2017); five pulsars whose $S_{150}$ values were likely overestimated due to lack of spectral data below ${\sim}{600}\,\mathrm{MHz}$ to reveal the presence of any spectral flattening or turnover; and nine that cannot presently be accounted for.
Detections were made across the frequency range 70.3– ${351.6}\,\mathrm{MHz}$ with both stations, although most pulsars were only detected near the middle of the frequency band where the sensitivity is highest. The bright pulsars J0953+0755 and J1645–0317 were detected across the majority of the frequency band, with only a couple of non-detections near the low and high ends of the frequency range. Radio-frequency interference (RFI) affected a small number of observations, which were either excised or discarded depending on the extent of the contamination. Five pulsars were detected by only one of the stations. The non-detections were eventually traced to either a malfunctioning data capturing system or RFI. Therefore, 17 pulsars were simultaneously detected by both stations and direct comparisons can be made between these detections.
The highest-DM pulsar in the detected sample was PSR J1731–4744 ( $\mathrm{DM}\approx{123.06}\,\mathrm{cm}^{-3}\mathrm{pc}$ ), which was detected down to ${132.8}\,\mathrm{MHz}$ . It is also particularly noteworthy that PSR J0437–4715 was detected by both stations—this MSP’s exquisite timing stability has placed it as a high-profile target for PTAs (e.g. Manchester et al. Reference Manchester2013), and its high flux density enables detailed studies into the remarkable evolution of its integrated pulse profile (Bhat et al. Reference Bhat2018).
A gallery of integrated pulse profiles showing the best detection of each pulsar is displayed in Figure 3. The stations have demonstrated the capability to make relatively high-S/N detections for several bright pulsars. Particularly notable are the detections of PSR J0953+0755, the strongest of which exceeded a S/N of 250 (with 128 phase bins). Pulse profiles with a high S/N were re-processed at a higher time resolution.
The S/N ratio measured by each station during a simultaneous observation can be used to directly compare the station performances. For each observation, the S/N ratio was estimated using psrchive. On average, the AAVS2 produced higher-S/N detections than the EDA2. This was most apparent near the upper end of the observing band where the log-parabolic SKALA-4.1 antennas used by the AAVS2 are noticeably more sensitive than the MWA-style bowtie dipoles used by the EDA2.
Whilst the low sensitivity of the stations is often insufficient to study the features of pulse profiles, the detections of PSRs J0437–4715 and J2048–1616 each show two distinct pulse components that are consistent with their known structures in the published literature (Bhat et al. Reference Bhat2018; Xue et al. Reference Xue2017). Also notable is PSR J1900–2600, which is known to exhibit a dramatic evolution in the shape of its pulse profile and a complex, multi-component structure at higher frequencies (Mitra & Rankin Reference Mitra and Rankin2008). The AAVS2 and EDA2 detections of this pulsar show good agreement with the profile published by Xue et al. (Reference Xue2017) from an MWA observation at ${185}\,\mathrm{MHz}$ .
A closer examination of the European Pulsar Network (EPN) DatabaseFootnote k and studies of the pulsar population at low frequencies (e.g. Bondonneau et al. Reference Bondonneau2020; Xue et al. Reference Xue2017; Stovall et al. Reference Stovall2015) reveals that only a small fraction of pulsars have published detections below ${400}\,\mathrm{MHz}$ , and therefore detection of even a modest sample can make a useful contribution to the growing database. In fact, six detections in this paper were made at the lowest frequencies published to date: PSRs J0835–4510Footnote l ( ${159.4}\,\mathrm{MHz}$ ), J1116–4122 ( ${150}\,\mathrm{MHz}$ ), J1453–6413 ( ${101.6}\,\mathrm{MHz}$ ), J1456–6843 ( ${101.6}\,\mathrm{MHz}$ ), J1731–4744 ( ${132.8}\,\mathrm{MHz}$ ), and J1751–4657 ( ${150}\,\mathrm{MHz}$ ). Each of these pulsars were also detected by Xue et al. (Reference Xue2017) at ${185}\,\mathrm{MHz}$ using incoherently summed MWA data, however the SKA-Low station detections were made down to even lower frequencies, in addition to having coherently summed station beams. Moreover, the detections of the Vela pulsar (PSR J0835–4510) are marginally lower in frequency than the MWA detections published in Kirsten et al. (Reference Kirsten2019) at ${164.5}\,\mathrm{MHz}$ .
4.2 Mean flux densities
Mean flux densities were measured for 22 pulsars using the methods described in Sections 3.2 and 3.3. The results are given in Table 1. For eight pulsars, detections were made at multiple frequencies and mean flux densities were measured for each detection. Multiple simultaneous detections were made of PSRs J0835–4510, J0953+0755, and J1645–0317. Additionally, most of the pulsar detections from the initial census were made with both stations. There is therefore a substantial number of simultaneous observations which can be used to independently verify the accuracy of the mean flux density measurements. Figure 4 shows a comparison of the mean flux densities at frequencies for which observations were performed with both stations simultaneously. Nearly all measurements agree within a factor of 2, despite the stations using different antenna types and relying on beam models with different accuracies. The plot also indicates the sensitivities of the stations, with detections being made down to ${\sim}100\,\mathrm{mJy}$ (for a ${\sim}10\sigma$ detection of the pulsar in 30 min, at ${150}\,\mathrm{MHz}$ ). The agreement between the stations gives us confidence that the mean flux density measurements are reasonably accurate and that their calculated uncertainties are not underestimated.
a Pulsar period listed in the ATNF pulsar catalogue.
b DM listed in the ATNF pulsar catalogue.
c Centre frequency of the observation.
d Observation duration.
e S/N of the integrated pulse profile obtained from EDA2 and AAVS2 data.
f Mean flux density calculated from the integrated pulse profile obtained from EDA2 and AAVS2 data.
* Millisecond pulsar.
4.3 Flux density spectra
The spectra of all 22 pulsars which were detected with the SKA-Low stations were modelled using the method described in Section 3.5. The data sets for each spectrum had a sufficient number of points to fit all model types and spanned a sufficiently wide frequency range to make a meaningful fit. The results for the pulsars which were observed at multiple frequencies are presented in Table 2 and Figure 5, and each of these fits are analysed in detail in Section 4.3.2. For the pulsars which were only observed at single frequencies, the results are given in Table 3 and Figure 6.
The most common spectral types were power-law with low-frequency turnover and broken power-law, with seven pulsars being best fit by each of these models, followed by five pulsars exhibiting a log-parabolic spectrum. The power-law with high-frequency hard cut-off was seen to be exhibited by two pulsars, and the simple power-law model by one pulsar. Of the 22 analysed pulsars, nine were not previously classified in JvSK+18, and for eight we have updated the best-fitting models. These pulsars are indicated in Tables 2 and 3. As a result of the current sensitivity limitation of SKA-Low stations, our detections are primarily for brighter pulsars, most of which have some low-frequency measurements available. Furthermore, all but one pulsar (J1731–4744) show deviations from a simple power-law model. This suggests that pulsars with well-determined spectra are more likely to show spectral features, which is in line with the findings from JvSK+18.
4.3.1 Modelling without continuum flux densities
Measurements of flux density from continuum images can have much lower measurement uncertainties than traditional flux density measurements from beamformed detections if the field-of-view is large enough to enable robust in-field calibration to be performed. However, continuum measurements can be vulnerable to systematic errors caused by blending with emission from other radio sources within the beam. It can therefore be difficult to determine the intrinsic flux of the pulsar (e.g. Bell et al. Reference Bell2016). Observations of targets near the Galactic plane ( $|b|\lesssim5^\circ$ ) are more likely to be affected by source confusion, and the problem may be amplified for larger beam sizes, such as the ${\sim}3\,$ arcmin MWA beam used for the GaLactic and Extragalactic All-sky MWA (GLEAM) survey (Wayth et al. Reference Wayth2015).
a Probability that the best-fitting model is the true best-fitting model.
b Flux density calculated from EDA2 and AAVS2 detections. For frequencies where detections were made with both stations, a weighted average is given.
c The spectral index (PL), the break frequency in MHz (BPL), the turnover frequency in MHz (PL-T), the peak frequency in MHz (LPS), or the cut-off frequency in MHz (PL-C).
d The spectral index below the break (BPL), the spectral index (PL-T and PL-C), or the curvature parameter (LPS).
e The spectral index above the break (BPL), the smoothness parameter (PL-T), or the spectral index for the case of zero curvature (LPS).
f The spectral classification. PL: simple power-law. BPL: broken power-law. LPS: log-parabolic spectrum. PL-T: power-law with low-frequency turnover. PL-C: power-law with high-frequency cut-off.
g The frequency range of the fitted data set.
LDue to the significant scattering tail of PSR J0835–4510, the flux densities below $\text{200}\,\text{MHz}$ were calibrated using upper noise limits, and may be underestimates (see Section 3.2 for details).
EPulsars whose flux densities were calibrated against full-electromagnetic simulations.
UBest-fitting model updated from JvSK+18.
NNot modelled by JvSK+18.
Only five of the detected pulsars lie on the Galactic plane: J0835–4510, J1453–6413, J1820–0427, J1932+1059, J2018+2839. Of these pulsars, none show notable discrepancies between continuum measurements and our measurements. Nevertheless, for the spectral modelling presented in Figures 5 and 6, we have considered flux densities measured from pulsed emission and continuum images as two separate classes of measurements. In addition to the models fitted to the full data sets, secondary models were fitted without the contributions from continuum flux density measurements, which includes all measurements from Murphy et al. (Reference Murphy2017) and Bell et al. (Reference Bell2016) (i.e. the secondary fits excluded data from these two publications). For the most part, the two data sets suggested similar models, with a few exceptions which are discussed in Section 4.3.2.
a See notes of Table 2 for details.
b Centre frequency of the observation. These frequencies were chosen based on where the S/N was expected to be highest.
c Flux density calculated from EDA2 and AAVS2 detections at centre frequency $\nu$ . For frequencies where detections were made with both stations, a weighted average is given.
U Best-fitting model updated from JvSK+18.
N Not modelled by JvSK+18.
4.3.2 Analysis of pulsar spectra
PSR J0437–4715:
Being the only MSP in the sample, the spectral behaviour of this pulsar is of particular interest. Our modelling shows that the best-fitting model is a broken power-law with a spectral break at ${1900} \pm {100}\,\mathrm{MHz}$ and a steep spectral index of ${-2.4}\,\pm\,{0.2}$ above the break frequency. This is in agreement with JvSK+18 and lowers the uncertainties of these two parameters significantly. Interestingly, our analysis suggests that the spectrum of this pulsar shows a low-frequency turnover at ${285}\,\pm\,{5}\,\mathrm{MHz}$ when continuum flux densities are excluded from the fit. It is important to note that MSPs are generally not seen to turn over, and those that hint at a turnover are expected to peak well below ${100}\,\mathrm{MHz}$ . Furthermore, this pulsar shows high variability at low frequencies due to interstellar scintillation. Bhat et al. (Reference Bhat2018) measured a characteristic scintillation bandwidth of ${1.4}\,\pm\,{0.4}\,\mathrm{MHz}$ , which makes the SKA-Low stations susceptible to intensity variations between observations. This may explain the discrepancy between the EDA2 measurements and the Murphy et al. (Reference Murphy2017) measurements (around a factor of 2).
PSR J0630–2834:
This pulsar shows perhaps the most striking example of continuum measurements dominating the spectral fit. The primary model fit indicates that the pulsar turns over at ${98}\,\pm\,{4}\,\mathrm{MHz}$ , however the model departs significantly from the majority of data points, with the exception of the Murphy et al. (Reference Murphy2017) data. Excluding continuum measurements, we find that the best-fitting model is a broken power-law with a break frequency of ${1000} \pm {200}\,\mathrm{MHz}$ , which tracks more closely with the majority of the available data. Evidently, this pulsar would benefit from more low-frequency measurements to clarify these inconsistencies.
PSR J0820–1350:
This pulsar has a well-defined spectrum, showing a power-law with a break frequency of ${440}\,\pm\,{60}\,\mathrm{MHz}$ , which is in agreement with JvSK+18 and reduces the uncertainty by 40%. Our modelling also shows shallower power-laws on either side of the break frequency than those found in JvSK+18, although the values agree within their respective uncertainties. Interestingly, the exclusion of imaging measurements leaves the spectrum poorly constrained below ${400}\,\mathrm{MHz}$ , leading to the broken power-law peaking at ${300}\,\pm\,{40}\,\mathrm{MHz}$ , with a positive spectral index below the break.
PSR J0835–4510:
Although the Vela pulsar is one of the brightest and most well-studied pulsars, if we limit ourselves to traditional (beamformed) pulsar observations, its radio spectrum remains relatively sparsely populated and poorly constrained at low frequencies. Low-frequency measurements are limited by interstellar scattering, which leads to substantial temporal broadening at frequencies below ${300}\,\mathrm{MHz}$ and renders pulses undetectable below ${\sim}{150}\,\mathrm{MHz}$ . At these frequencies, the only measurements available are from continuum imaging. At higher frequencies (between ${600}\,\mathrm{MHz}$ and 10 GHz), the only strong constraints are provided by JvSK+18, Johnston & Kerr (Reference Johnston and Kerr2018), Zhao et al. (Reference Zhao2019), and van Ommen et al. (Reference van Ommen, D’Alessandro, Hamilton and McCulloch1997). Above 10 GHz, the model departs significantly from the data. Notably, this fit includes measurements from Mignani et al. (Reference Mignani2017) at millimetre and even sub-millimetre wavelengths made using the Atacama Large Millimeter/submillimeter Array (ALMA). Very few pulsars have been observed at these extremely high frequencies, and it is likely that the models considered in this work do not adequately describe the behaviour of the full spectrum. A more complex model such as a double broken power-law (see, for example, Bilous et al. Reference Bilous2016) may be more suitable. Notwithstanding, our modelling shows a broken power-law with a spectral break at ${970}\,\pm\,{50}\,\mathrm{MHz}$ and a spectral index below the break of ${-0.50}\,\pm\,{0.03}$ , which are in agreement with JvSK+18. Above the break, our modelling suggests a steeper spectrum than reported in JvSK+18, with a spectral index of ${-2.9}\,\pm\,{0.2}$ .
PSR J0837+0610:
This pulsar has a substantial amount of data available over a moderately wide frequency range. Our modelling shows that the best-fitting model is a log-parabolic spectrum that peaks at ${65}\,\pm\,{6}\,\mathrm{MHz}$ . It is clear that the measurements from Murphy et al. (Reference Murphy2017) dominate the spectral fit, and that the reliability of the model fit is heavily dependent upon these continuum measurements. Our flux density measurements fall in line with these measurements, barring some variability resulting from interstellar scintillation, an effect typically seen in low-DM pulsars such as this ( $\mathrm{DM}\approx 12.86\,\mathrm{cm}^{-3}\,\mathrm{pc}$ ). Nevertheless, excluding continuum measurements, we find that the spectrum exhibits a broken power-law with a break at ${185.00}\,\pm\,{0.04}\,\mathrm{MHz}$ , constrained strongly by a low-uncertainty measurement from Xue et al. (Reference Xue2017).
PSR J0953+0755:
Being one of the nearest and brightest pulsars, J0953+0755 has one of the best-determined radio spectra, with data available from 17 publications in our database and spanning three orders of magnitude in frequency. Our modelling shows that this pulsar exhibits a low-frequency turnover at ${104}\,\pm\,{7}\,\mathrm{MHz}$ , which is marginally higher than that which was found in JvSK+18. This pulsar is also known to exhibit significant variability (see Bell et al. Reference Bell2016), due to which measurements based on single observations may ideally need to account for additional sources of errors (typically a factor of ${\sim}2$ –3 at these frequencies). It is likely that scintillation is also the cause of the variability seen in measurements from Sieber (Reference Sieber1973) and Izvekova et al. (Reference Izvekova, Kuzmin, Malofeev and Shitov1981) at frequencies below ${\sim}{1}\mathrm{GHz}$ . For this spectrum, the contribution from continuum measurements is minimal.
PSR J1453–6413:
The best-fitting model for this pulsar is a steep power-law with a spectral index of ${-2.5}\,\pm\,{0.3}$ and a low-frequency turnover at ${140}\,\pm\,{40}\,\mathrm{MHz}$ . For this pulsar, our measurements are the lowest frequency available, and thus important in constraining the turnover. This pulsar was previously found to have a broken power-law with a break at ${320}\,\pm\,{30}\,\mathrm{MHz}$ (JvSK+18), however our data suggest that there may be significant curvature below ${300}\,\mathrm{MHz}$ .
PSR J1456–6843:
For this pulsar, we found that the best-fitting model is a log-parabolic spectrum with a peak at ${\lesssim}{30}\,\mathrm{MHz}$ ; however, it is unlikely that the spectrum truly peaks at such a low frequency. The far southern declination of the pulsar is perhaps the most likely reason for the relatively small number of measurements. More data would be beneficial to better constrain the spectrum.
PSR J1645–0317:
Our modelling shows that this pulsar exhibits a broken power-law with a peak at ${276}\,\pm\,{2}\,\mathrm{MHz}$ and a shallow positive spectral index of ${0.22}\,\pm\,{0.01}$ below the break. Above the break, the power-law has a spectral index of ${-2.15}\,\pm\,{0.01}$ , which is shallower than the spectral index of ${-2.6}\,\pm\,{0.2}$ reported in Johnston et al. (Reference Johnston2021). Excluding the continuum measurement from Murphy et al. (Reference Murphy2017), the peak of the broken power-law shifts to ${220}\,\pm\,{2}\,\mathrm{MHz}$ , with a steeper spectral index of ${0.89}\,\pm\,{0.01}$ below the break. Our data fills a large gap in the spectrum of this pulsar and plays an important role in constraining its spectrum.
PSR J1731–4744:
Our observations of this pulsar have yielded the lowest-frequency measurements in the spectrum, and our analysis hints at flattening or turnover at ${\lesssim}{200}\,\mathrm{MHz}$ . However, the large uncertainties of the measurements mean limited constraints on the model, which we find to be a simple power-law with a relatively shallow spectral index of ${-1.38}\,\pm\,{0.05}$ .
PSRs J0452–1759 and J1116-4122:
These are the only pulsars in our sample to show a high-frequency hard cut-off, a model proposed by Kontorovich & Flanchik (Reference Kontorovich and Flanchik2013) to describe the coherent emission of electrons when accelerated by the pulsar’s electric field. Following JvSK+18, we calculate the magnetic field strength in the centre of the polar cap as
where $m_e$ and e are the mass and charge of the electron, c is the speed of light, P is the pulsar period, and $\nu_c$ is the cut-off frequency obtained from the spectral fit. If we assume a dipolar magnetic field, the field strength drops off as $z^{-3}$ , where z is the distance from the centre of the neutron star. We can thus derive the emission height (i.e. the height of the centre of the polar cap) $z_e$ as
where $z_{\mathrm{surf}}$ is the radius of a canonical 1.4 $\,\mathrm{M}_\odot$ neutron star ( $12\,\pm\,2\,\mathrm{km}$ ; Steiner et al. Reference Steiner2018) and $B_{\mathrm{surf}}$ is the magnetic field strength at the surface of the star. Using the cut-off frequencies given in Table 3, we calculated $B_{\mathrm{pc}}$ and $z_e$ for both pulsars. The results are presented in Table 4, along with the magnetic field strength at the light-cylinder radius $B_\mathrm{LC}$ and the emission height normalised by the light-cylinder radius $z_e/R_\mathrm{LC}$ . Our fits predict $z_e\sim30$ – $50\,\mathrm{km}$ , which is in line with JvSK+18, who noted that the estimated emission heights are unreasonably low and that further constraints will be required at high frequencies to test the validity of the model and its underlying assumptions. Furthermore, for PSR J0452–1759, our fitted spectral index of ${-0.1}\,\pm\,{0.1}$ departs significantly from the assumed $\alpha=-2$ in Kontorovich & Flanchik (Reference Kontorovich and Flanchik2013). More data would be beneficial to better constrain this spectrum and clarify this inconsistency.
a Values listed in the ATNF pulsar catalogue.
5. Summary and conclusions
We have presented the first pulsar detections with the second-generation SKA-Low precursor stations, the AAVS2 and the EDA2. By performing a shallow all-sky census of the southern sky, we have detected 22 bright known pulsars at low frequencies. With the current modest station sensitivities ( ${\sim}{1}\,\mathrm{MHz}$ bandwidth and 30 min integration time), the stations are sensitive down to ${\sim}100\,\mathrm{mJy}$ for pulsar detections. Furthermore, follow-up observations were made across a wide frequency range ( ${\sim}70$ –350MHz), yielding detections at multiple frequencies for eight pulsars. Six pulsars observed in this work were detected at the lowest frequencies ever published (J0835–4510, J1116–4122, J1453–6413, J1456–6843, J1731–4744, J1751–4657). In general, the AAVS2 produced higher-S/N pulsar detections than the EDA2, which is likely a result of the differences in antenna types between the stations. We expect that the upcoming system bandwidth upgrade ( ${\gtrsim}{25}\,\mathrm{MHz}$ ) will enable a fivefold increase in sensitivity, and yield high-quality pulsar detections with only a fraction of the processing time of the MWA, making the stations strategically useful tools for regular pulsar monitoring.
Flux densities were measured for all 22 detected pulsars, which show close agreement between the two stations and with other measurements in literature. This gives us confidence that the sensitivity simulation code used to calibrate the observations produces consistently reliable estimates. Comparisons with full-electromagnetic simulations confirm this, with a typical discrepancy of no more than 30%.
We have revisited the spectral modelling of 13 pulsars previously analysed in JvSK+18 and further analysed an additional 9 by employing the robust modelling and classification method described in the original publication. By augmenting our results with a compilation of flux density measurements from the published literature, we have classified the analysed spectra into 5 morphological classes. For 17 pulsars, we have presented new or updated best-fitting models. In our sample of bright pulsars, all but one spectrum showed deviations from simple power-law behaviour, which suggests that populating the low-frequency spectra of pulsars will often reveal hidden spectral features. Furthermore, we have found that even very well-studied pulsars such as the Vela pulsar (PSR J0835–4510) often have poorly constrained spectra and could benefit from further measurements to fill gaps in the data and clarify inconsistencies between published values. Thus, it is important for more pulsar flux density measurements to be made at low frequencies, ideally with instruments that cover a wide frequency range to limit the mismatches in systematic errors introduced with combined data sets.
This work also provides an excellent demonstration that the capabilities of the SKA-Low precursor stations and pathfinder instruments (e.g. the MWA, LOFAR) can be employed to gain an improved understanding of pulsar radio spectra. A deeper analysis of the modelled spectra can be used to draw connections to the pulsar emission mechanism and to further constrain the location of the emission region. Moreover, the spectral modelling presented in this work can be expanded to a larger sample of pulsars to further our understanding of the low-frequency behaviour of pulsar spectra. These results will also help to inform pulsar population studies and make improved predictions of the expected yields of large pulsar surveys, which is particularly important in preparation for science with the SKA-Low (expected to be operational by the end of the decade).
Acknowledgement
We thank the referee for several useful comments that helped to improve this paper. We also thank Sam McSweeney for useful discussions. CPL was supported by an ICRAR Studentship and an EECMS Summer Scholarship through Curtin University. This scientific work makes use of the Murchison Radio-astronomy Observatory, operated by CSIRO. We acknowledge the Wajarri Yamatji people as the traditional owners of the Observatory site. We acknowledge the Pawsey Supercomputing Centre, which is supported by the Western Australian and Australian Governments. The AAVS2 and EDA2 are hosted by the MWA under an agreement via the MWA External Instruments Policy. The acquisition system was designed and purchased by INAF/Oxford University and the RX chain was design by INAF, as part of the SKA design and prototyping program. We acknowledge the support of the Curtin operations team, INAF group, and SKA-Low team in the development and ongoing maintenance of the facilities used in this work. We acknowledge the use of the following software/packages for this work: psrchive (Hotan et al. Reference Hotan, van Straten and Manchester2004) and dspsr (van Straten & Bailes Reference van Straten and Bailes2011); and the python packages iminuit (Dembinski et al. Reference Dembinski2020), psrqpy (Pitkin Reference Pitkin2018), astropy (Astropy Collaboration et al. Reference Astropy Collaboration2018), matplotlib (Hunter Reference Hunter2007), and numpy (Harris et al. Reference Harris2020). This work made use of NASA’s Astrophysics Data System and arXiv.