Simulating high-time resolution radio-telescope observations

Radio telescope observations are affected by radiometer noise. We model the frequency-dependent system noise as being drawn from a Gaussian distribution with amplitude defined by the radiometer equation:

where $G_{\rm dpfu}$ is the telescope gain in degrees per flux unit ($\rm K/Jy$), $T_{\mathrm{sys}}$ is the system temperature, $\tau$ the digital sampling time and $\Delta\nu$ the receiver bandwidth. Throughout the simulation software we assume that the signals and noise are represented in units of janskys.

High-time resolution data sets can be affected by low-frequency noise processes. Such noise can arise from gain variations, telescope-pointing jitter, spill-over variations, etc. Such noise can impact on the detection of radio transients like FRBs (e.g., Zhang et al. 2021) and requires de-reddening procedures in pulsar searches (e.g., Lazarus et al. 2015; Suresh et al. 2022). We provide the ability to model red noise with power spectral density

Here $\alpha$ is the spectral index and $S_{0}$ the power spectral density at 1 Hz. Our code also allows for a cut-off frequency in the simulated red noise (the modelled power spectral density below the cut-off is zero) or a frequency at which the spectrum flattens (the power spectral density below the flattening frequency equals the power at that frequency). An example is shown in Figure 2 where red and white noise has been simulated, the observing-frequency channels summed to form a time series and then Fourier transformed.

Radio frequency interference (RFI) is defined here as unwanted radio signals in the astronomical data sets (which may be self-generated at the observatory or from external sources). RFI is ubiquitous and dealing with RFI is a major challenge in current radio astronomy. We can divide RFI into three categories: 1) impulsive and wideband RFI, 2) narrowband and persistent RFI and 3) birdies (which produce periodic signals across the band from, for instance, the mains electricity supply). All these types of RFI can be simulated using our software.

Impulsive, broadband RFI is simulated based on a distribution of amplitudes and pulse widths detected at a specific observatory. The events are randomly distributed in time through the simulation. We used an analysis of the zero dispersion measure, impulsive RFI detected in representative observations of the High Time Resolution Survey (HTRU; Keith et al. 2010) to determine the typical amplitudes and widths of impulsive RFI detected at the Parkes observatory. Such events are shown as the narrow (only one or two time-samples wide), vertical stripes in Figure 3.

We can provide a simple simulation of the signals from transmission towers by increasing the system noise significantly in frequency bands where there is a strong transmitter. Examples are shown as the horizontal, relatively narrow-band features in Figure 3. The detectability of any astronomical source (such as the flare star simulated in this data set) in such bands will be significantly reduced (and often undetectable).

Most RFI is time dependent; either through intrinsic time variability or because the telescope pointing direction is not constant with respect to the location of the interfering source. We do not attempt a detailed simulation of all known types of time variable RFI. For instance, mobile transmission towers and handsets, WiFi, BlueTooth, aircraft navigational devices all produce strong signals that are time variable. Our simulation code allows for relatively simple representations of such signals. One example is shown near the centre of Figure 3 and represents a simplified representation of the signals observed from a mobile transmission tower.

The Parkes ultra-wide-bandwidth receiver system (Hobbs et al., 2020) is affected by RFI. Some of the most problematic are mobile transmission towers that produce impulsive RFI that can be extremely bright for short durations. Analogue-to-digital converters are close to being limited by the jitter noise of the sample clock (see Tuthill et al., in prep.). This leads to power in one part of the digitised band being “smeared” across the entire band. The Parkes data sets are affected by this issue and we are trialling various methods to mitigate this effect. In Figure 4 we show a simple simulation of 4G mobile transmission, which is highly impulsive in small sub-bands. We smear a tiny fraction ($10^{-8}$) of the signal strength across the entire band, which produces a signal that is similar to our actual observations in some sky directions.

Signals from satellites are a major RFI source at all observatories. The positions of individual satellites can be predicted using two line element (TLE) sets. However, for many of the satellites (such as the global positioning satellite systems), there are a large number of satellites above the horizon at a given time and the signal from each satellite is detected in the far side lobes of the telescope beam. To simulate such sources we assume that each satellite emits multiple signals each of which follows a sinc²-frequency-dependent function:

where $f$ is the observing frequency and $f_{0}$ is one of the emission frequencies from the satellite. Instead of attempting to model the structure of the far side-lobes of any particular telescope, we assume that the motion of the source through the far side-lobes of the telescope will introduce a sinusoidal variation in the signal amplitude:

where $A^{\prime}$ is the signal amplitude, $P$ is the variability timescale and $\phi$ is the signal phase. We show an example (close to 1500 MHz) in Figure 3 of the RFI caused by the global positioning satellites.

To simulate tones, we inject a sinusoidal signal (or, if requested, a rectified sinusoidal signal) with a given amplitude, phase and frequency. Such signals are typically only apparent after a periodicity search has been carried out on the final data product and examples of common signals are shown in Figure 3. In contrast some recent (or planned) surveys include strong periodic signals that are purposely injected into the data stream. This is because commensal high-time resolution and high frequency resolution surveys are now being carried out in order to maximise telescope efficiency by searching for pulsars, FRBs and spectral lines simultaneously (for instance, the CRAFTS survey being carried out with FAST; Li et al. 2018). Spectral line surveys traditionally make use of a switched noise source throughout their observations to track and account for telescope gain variations and for amplitude calibration. The simulation code can therefore include such a broad-band, switching noise source and we provide an example in Figure 3.

We know that high-time resolution data sets will include individual bright pulses from pulsars, fast radio bursts (FRBs) and related sources such as rotating radio transients (RRATs). Such transient events are parameterised using an event time, pulse amplitude, pulse width and dispersion measure (DM). The pulse shape can be simple (modelled using a single Gaussian or rectangular profile) or defined using multiple Gaussian components. We can also build band-limited (Kumar et al., 2021) and downward-drifting structure (Hessels et al., 2019) in the burst spectrum. Examples are shown in Figure 5 for bursts with different profile and frequency structure.

The aim of many high-time resolution surveys is to find fast transients (such as pulsars and FRBs). The use of low-bit quantisation in the output data products has required level setting procedures to account for longer term system noise variations (more details are provided in Section 4.2). This has ruled out the chance of finding slower transient signals. However, several classes of stars across the Hertzprung-Russell Diagram, including the Sun, produce intense bursts of non-thermal radio emission, powered by various forms of magnetic activity (e.g., Güdel, 2002). These include auroral activity (Zarka, 1998; Hallinan et al., 2015; Trigilio et al., 2011), coronal activity, driven by flaring, space weather, and other dynamic processes within stellar atmospheres and astropheres (e.g., Bastian, 1990; Pick & Vilmer, 2008; Benz, 2017); and the interaction of magnetic fields between two components in close binary systems such as RS Canum Venaticorum systems (Drake et al., 1989; Toet et al., 2021). Regardless of the driving mechanism, active stars can produce radio emission variable on timescales from milliseconds (e.g., Osten & Bastian, 2008) to days (Slee et al., 2003), with complex time-frequency structure. The ability to capture this variability across such a broad range of timescales remains under-explored. As such signals are likely to exist in our current archival data sets, yet the pulsar-based algorithms developed so far are unlikely to detect them, we have included the ability to simulate likely flare-star signatures. We parameterise such events as relatively broad-band sources covering a specified bandwidth, where individual components drift in time and frequency following a quadratic polynomial (allowing both linear and quadratic drifts). An example is shown in Figure 3.

Pulsars produce a sequence of periodic pulses. Apart from the brightest known pulsars, the majority of these pulses are so weak that they cannot individually be detected. Instead the search-mode data streams are de-dispersed, averaged across frequency channels and then Fourier transformed to search for periodic signatures within the data.

We simulate pulsars using a tempo2-style predictor (Hobbs et al., 2006) that can be used to determine the arrival times of pulses at a specific observing frequency and specific observatory. The use of predictors allows highly relativistic binary systems to be modelled. For instance, in Figure 6 we compare an actual observation with the simulation for PSR J0737$-$3039A, which was discovered by Parkes high-latitude multibeam pulsar survey (Burgay et al., 2003) and is a highly relativistic binary system (Kramer et al., 2021). The real data set was obtained from the CSIRO Data Access Portal⁴⁴4https://data.csiro.au/collections (DAP, Hobbs et al. 2011) and we folded the data at a nominal rotational period (upper panel in Figure 6). Our simulated result (lower panel) is based on a tempo2 predictor from an up-to-date timing ephemeris for this pulsar.

The pulse profile can be modelled using multiple Gaussian components. The resulting pulse shape is subsequently rescaled to ensure the area represents the mean flux density of the pulsar. Each individual pulse from the pulsar is simulated. We therefore can vary the intensity of each individual pulse. As the intensities of individual pulses are drawn from various distributions (as described for terrestrial use by Dawson et al. 2022), we allow the user to provide their own distribution that is then used within the simulation.

The observed astronomical signals have propagated through the interstellar medium (ISM). The time delay caused by dispersion for a pulse measured at two different frequencies of radio waves ( $\nu_{1}$ and $\nu_{2}$) is given by:

The dispersion measure (DM) is defined as the integral of electron density along the light of sight,

Typically $\alpha=-2$ for cold, diffuse, ionised gas. Dispersion is simulated in our software, where $\alpha$ and DM can be defined by the user. Note that negative dispersion measures are permissible, as they may be produced by astronomical signals with intrinsically negative drift.

The signals are also affected by scintillation. Various models for scintillation are possible and our simulation software allows the user to develop their own models of scintillation as required. Here we consider only diffractive scintillation in the case of strong scattering. For given scintillation time-scale and bandwidth, we simulate a dynamic spectrum following procedures described in Dai et al. (2016). In the dynamic spectrum the flux density as a function of time and frequency is $S_{\rm dyn}(t,\nu)$. We determine this from a 2-dimensional auto-covariance function defined using a user-provided scintillation timescale ($\tau$) and bandwidth ($\Delta\nu$) (Dai et al., 2016).

We note that FRB signals often have frequency-dependent amplitude fluctuations that cannot be modelled through scintillation. We therefore provide options for the user to specify any frequency evolution of the pulse events. An example of a FRB that is constrained to a small frequency range is shown as the event of example (4) in Figure 5.

Pulse events also undergo scattering in the interstellar medium. We therefore provide the ability to convolve a pulse signal with an exponential function with a time scale defined by a specified DM. We use the pulse-broadening function in Bhat et al. (2003) to model the scattering effect on integrated pulse profile,

where $\tau_{\mathrm{d}}$ is the pulse broadening time and $U(t)$ is the unit step function. To demonstrate the scintillation and scattering simulations we model a pulsar with a dispersion measure of 600 $\mathrm{cm}^{-3}\,\mathrm{pc}$, a diffractive bandwidth of 50 MHz and diffractive timescale of 1 minute (we note that these are not necessarily independent parameters for actual pulsars, but choose these parameters here to allow us to demonstrate multiple ISM-related effects in a single figure). The resulting profile is shown in Figure 7.

One aim of our simulation software is to provide tools to simulate unexpected, or unknown, sources. We simulate such “unknown unknown” signals in the following three ways: (1) generalising the burst event simulations, (2) using spline curves and (3) generic images ⁵⁵5In 1974, the active SETI used the Arecibo radio telescope to transmit a message towards the globular cluster M13. This message consisted of multiple generic images digitised in 1679 bits (Atri et al., 2011)..

We generalise the burst events by allowing positive or negative dispersion measures and allowing any choice of $\alpha$ in Equation (5). We also allow any user-defined frequency evolution of the burst intensity.

The high-time resolution data streams simulated here typically have relatively low frequency resolution. Therefore any unexpected source signals that will be detected in such data sets are likely to be relatively short in duration, but cover a wide band. We therefore allow for arbitrary broadband signatures to be simulated using cubic splines, which are defined by specific time and frequency control points.

We can embed information into the data set by producing an image (which could be obtained from the large number of available online datasets used for training ML algorithms) and using the pixels in that image to represent the time-frequency information in the simulated data set.

We provide an example of these three types of “unknown unknown” signals in Figure 8. This contains, from left to right, an arbitrary cubic spline event, a cartoon image of a cat and an FRB-like event with a negative dispersion measure.

The simulated source signals are either only present in specific sky directions (different pulsars are at different sky coordinates) or are always present (for instance, radiometer noise). In a multi-beam system, each receiver beam can be modelled independently and the sky position of that beam can be defined as a function of time (allowing for simulations of scanning surveys or tracking observations of a specific sky direction). The beam pattern on the sky is given by the telescope diameter and observing frequency and assumed to follow:

where $\theta=1.22\lambda/D$, $\lambda$ is the observing wavelength and $D$ the telescope diameter. As an example, we show in Figure 9, a simulation of the detection of the first FRB (the “Lorimer burst”). The left panel shows the actual data from the Parkes telescope (Lorimer et al., 2007) and the right panel shows our simulation. The event is seen primarily in Beam 6, but is also present in other beams (more detectable at lower frequencies where the beam is wider). Note that we model narrow-band interference as being detectable in all of the beams.

Our simulation software can also model a drift-scan survey in which the telescope is held fixed and astronomical sources will drift through a beam. In Figure 10 we show how a pulsar will be detected in such a survey. We have assumed here a sensitive, large telescope and show the pulsar pulses being detected in the side lobes as well as in the primary telescope beam.

The astronomy signal processor processes the incoming data streams and produces the final data products. Typically this involves channelising the data streams. Dispersion smearing within a frequency channel is modelled by first simulating more channels than required and then averaging those channels to the requested output channelisation. If a known pulsar is being observed then the data may first be coherently de-dispersed at the known DM of the pulsar and hence, in this case, there will be no channel-dependent dispersion smearing.

The data volumes can be enormous and so the output data samples are typically written using only 1-, 2- or 8-bit quantisation. Low-bit quantisation requires knowledge of the typical digital sample levels, which may change through an observation. The levels can be pre-defined (e.g., any time sample above zero is set to 1 and any signal below to 0 for 1-bit data), or user-defined. It is also possible to model automatic level-setting procedures. For instance, a running mean for each frequency channel (of specified number of samples) can be used to define the levels. This produces the observed change in the noise properties of the Lorimer burst detection after the event (in Figure 9) and was used in the Parkes multibeam surveys that used an analogue filterbank system (e.g,. Manchester et al. 2001). The levels can also be set from the first samples and then held fixed for the remainder of the observation. This is similar to the level setting used in the HTRU (Keith et al., 2010) and Survey for Pulsars and Extragalactic Radio Bursts (SUPERB; Keane et al. 2018) at the Parkes radio telescope.

Note that with low-bit digitisation it is common for a bright signal to saturate the system. This is shown in Figure 10 where the pulses being detected by the simulated primary beam are so strong that the time series saturates.