#  setigen: Simulating Radio Technosignatures for SETI

Bryan Brzycki Department of Astronomy, University of California Berkeley,
Berkeley CA 94720  Andrew P. V. Siemion Breakthrough Listen, University of
California Berkeley, Berkeley CA 94720  SETI Institute, Mountain View, CA
94043  Department of Physics and Astronomy, University of Manchester, UK
University of Malta, Institute of Space Sciences and Astronomy, Malta  Imke de
Pater Department of Astronomy, University of California Berkeley, Berkeley CA
94720  Steve Croft Department of Astronomy, University of California Berkeley,
Berkeley CA 94720  SETI Institute, Mountain View, CA 94043  John Hoang
Department of Astronomy, University of California Berkeley, Berkeley CA 94720
Cherry Ng Dunlap Institute for Astronomy & Astrophysics, University of
Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada  Department of
Astronomy, University of California Berkeley, Berkeley CA 94720  SETI
Institute, Mountain View, CA 94043  Danny C. Price International Centre for
Radio Astronomy Research, Curtin University, Bentley WA 6102, Australia
Department of Astronomy, University of California Berkeley, Berkeley CA 94720
Sofia Sheikh Department of Astronomy, University of California Berkeley,
Berkeley CA 94720  Zihe Zheng Goergen Institute for Data Science, University
of Rochester, Rochester, NY 14627

###### Abstract

The goal of the search for extraterrestrial intelligence (SETI) is the
detection of non-human technosignatures, such as technology-produced emission
in radio observations. While many have speculated about the character of such
technosignatures, radio SETI fundamentally involves searching for signals that
not only have never been detected, but also have a vast range of potential
morphologies. Given that we have not yet detected a radio SETI signal, we must
make assumptions about their form to develop search algorithms. The lack of
positive detections also makes it difficult to test these algorithms’ inherent
efficacy. To address these challenges, we present setigen, a Python-based,
open-source library for heuristic-based signal synthesis and injection for
both spectrograms (dynamic spectra) and raw voltage data. setigen facilitates
the production of synthetic radio observations, interfaces with standard data
products used extensively by the Breakthrough Listen project (BL), and focuses
on providing a physically-motivated synthesis framework compatible with real
observational data and associated search methods. We discuss the core routines
of setigen and present existing and future use cases in the development and
evaluation of SETI search algorithms.

astrobiology — technosignature — SETI — extraterrestrial intelligence

††software: NumPy (Oliphant, 2006), CuPy (Okuta et al., 2017), SciPy (Virtanen
et al., 2020), Astropy (Robitaille et al., 2013), Blimpy (Price et al., 2019),
H5py (Collette et al., 2017), Matplotlib (Hunter, 2007)

##  1 Introduction

Since the inception of radio SETI in the 1960s, technosignature searches have
greatly expanded to cover more sky area, wider frequency ranges, and a larger
variety of signal morphologies (Drake, 1961; Werthimer et al., 1985; Tarter,
2001; Siemion et al., 2013; Wright et al., 2014; MacMahon et al., 2018; Price
et al., 2018; Gajjar et al., 2021). Arguably the most developed branch of
radio SETI is the search for narrow-band technosignatures, with signal
bandwidths under 1 kHz, for which search algorithms are constantly being
produced and improved (Siemion et al., 2013; Enriquez et al., 2017; Pinchuk et
al., 2019; Margot et al., 2021). These algorithms operate on either voltage
time series data or time-frequency spectrogram data (i.e., dynamic spectra,
waterfall plots).

The incoherent tree deDoppler method is the primary search strategy for
Doppler-accelerated narrow-band signals in radio spectrograms (Taylor, 1974;
Siemion et al., 2013; Enriquez et al., 2017; Margot et al., 2021). An ideal
sinusoidal emitter will appear to exhibit a frequency drift over time due to
relative acceleration between the emitter and receiving telescope (Sheikh et
al., 2019). Under a constant relative acceleration, such a signal will have a
linear drift or slope in a spectrogram of Stokes I intensities. The tree
deDoppler algorithm efficiently integrates spectra over potential drift rates
and identifies signals above a threshold signal-to-noise ratio (SNR).
Breakthrough Listen, the most comprehensive SETI search program to date
(Worden et al., 2017), developed
turboSETI111https://github.com/UCBerkeleySETI/turbo_seti, an open-source
implementation of the deDoppler algorithm that serves as the backbone of many
technosignature searches (Enriquez et al., 2017; Enriquez & Price, 2019; Price
et al., 2020; Sheikh et al., 2020; Gajjar et al., 2021).

This method works well for signals with high duty cycles and linear drift
rates, but it can struggle to properly detect more complex signals (Pinchuk et
al., 2019). This is particularly problematic given the increasingly complex
radio frequency interference (RFI) environment within which these searches are
conducted. Moreover, the lack of robust, labeled, narrow-band signal datasets
can make it difficult to quantify a given implementation’s detection accuracy,
especially in light of RFI and variable bandpass responses.

For more complex signal morphologies, machine learning (ML) algorithms have
been proposed that use computer vision techniques to classify image-like
spectrograms. However, the same lack of labeled, narrow-band signal data makes
creating supervised ML models difficult. Zhang et al. (2019) used a self-
supervised approach in which spectrogram data was divided in time into two
halves, for which the ML task was to predict the second half given the first.
For an ML-based direction-of-origin filter, Pinchuk & Margot (2021) used a
separate non-ML method to detect signals and create an algorithmically-labeled
spectrogram dataset. In most cases, however, supervised approaches have relied
on generating synthetic signals of various classes in order to guarantee
correct labels (Harp et al., 2019; Brzycki et al., 2020; Margot et al., 2021).

To address these issues, we present setigen, an open-source Python library
that facilitates the creation of synthetic narrow-band signals and supports
injection into observational data. setigen is meant to provide a general-use
heuristic framework for creating mock radio SETI data. A primary design aspect
is ensuring that the synthesis process is grounded as much as possible in
physical quantities to better interface with real observations and search
algorithms. setigen makes heavy use of NumPy222https://numpy.org/ for
efficient matrix operations (Oliphant, 2006; Harris et al., 2020) and
blimpy333https://github.com/UCBerkeleySETI/blimpy for interfacing with data
products routinely used by BL (Price et al., 2019).

There are two main modules in setigen, “spectrogram” and “voltage,” dedicated
to the most common data formats used in radio SETI. The spectrogram module
works with Stokes I (intensity) data stored as time-frequency arrays and is
designed to be flexible and heuristic-based. It can be used to generate many
small snippets of data containing synthetic signals for quick algorithm test
cases or for full labeled datasets. The voltage module creates synthetic
antenna voltages, follows these voltages through a software-based signal
processing chain that models a standard single dish signal pipeline, including
quantization and a polyphase filterbank, and saves the final complex voltages.
This requires a lot more computational power, so voltage setigen routines can
be optionally GPU-accelerated via CuPy444https://cupy.dev/ (Okuta et al.,
2017). Since the voltage module models the signal processing chain, it can be
used to produce more “realistic” signals, test complex voltage processing
software, and evaluate how each signal processing element affects the final
signal sensitivity.

Radio SETI searches typically operate on data in spectrogram format, since it
compresses data and enables visualization and analysis of broader signal
morphology in time-frequency space (Enriquez et al., 2017; Margot et al.,
2018; Pinchuk et al., 2019; Price et al., 2020; Sheikh et al., 2020). As such,
setigen was initially written to create large datasets of radio spectrograms
for use in supervised ML search experiments. The library was later expanded to
support synthesizing raw voltage-level data to complement existing use cases.

setigen has already been used in a variety of applications, such as the
development and testing of search algorithms. It has been used to create
synthetic datasets with position labels for ML localization tasks in single
observations (Brzycki et al., 2020). setigen has also been used to inject
synthetic signals within ON-OFF cadences, each comprised of 6 consecutive
observations and used as a direction-of-origin filter for SETI. Ma et al.
(submitted) injected signals into ON-OFF cadences taken with the Robert C.
Byrd Green Bank Telescope (GBT; MacMahon et al., 2018) to train a
sophisticated variational autoencoder model that can classify cadences as
potential SETI candidates. Similarly, setigen was used extensively to produce
training and test data in BL’s first Kaggle ML
competition555https://www.kaggle.com/c/seti-breakthrough-listen, in which
contestants were tasked with classifying synthetic technosignature candidates
in ON-OFF cadences.

Outside of ML, synthetic setigen data is used in injection-recovery testing
for turboSETI as well as for a new search code,
hyperseti666https://github.com/UCBerkeleySETI/hyperseti. The voltage module
has been used to test and upgrade parts of the Allen Telescope Array’s (Welch
et al., 2009) software signal processing pipeline. Furthermore, setigen has
been used to test RFI rejection and detection techniques for the Parkes
Multibeam Galactic Plane Survey SETI search, helping to discriminate
terrestrial signals from different regions in the sky as SETI surveys with
multiple antennas or beams become more popular (Perez et al., in prep).

This paper is organized as follows. Section 2 outlines the standard signal
chain and processing pipeline used in single dish radio SETI observations to
motivate details behind setigen’s synthesis methods. Section 3 presents the
code methodology: Section 3.1 describes the spectrogram module for producing
and working with synthetic Stokes I time-frequency data, while Section 3.2
describes the voltage synthesis module in detail, connecting components of
typical radio signal chains to software analogues used in setigen. In Section
4, we discuss current limitations of the library and future directions for
signal synthesis for SETI.

##  2 Overview of Single Dish Signal Chains

To motivate the capabilities of setigen, we first give a broad overview of the
standard single dish data recording pipeline, as well as some details
pertinent to the Breakthrough Listen digital recorder (BL DR) system at the
GBT (MacMahon et al., 2018).

In a single-dish radio telescope, incoming radiation is reflected off the dish
surface toward a feed horn at the focus. The feed couples incident free-space
electromagnetic radiation to voltages within the telescope’s receiver system.

These voltages are passed to an analog down-conversion system containing a
heterodyne mixer, which shifts the signal from the target RF range into an
intermediate frequency (IF) range near baseband more suitable for receiver
hardware. The resulting voltages are then digitized by analog-digital
converters (ADC) to a specified number of bits $N_{\text{bits,d}}$ at a given
sampling rate $f_{s}$. The BL DR system digitizes voltages to 8-bit at a
sampling rate of $f_{s}=3$ GHz for each linear polarization (MacMahon et al.,
2018).

Radio telescope pipelines commonly use polyphase filterbanks (PFB; Bellanger
et al., 1976; Harris & Haines, 2011; Price, 2021) to help partition the usable
band and improve the spectral channel response of the system. For example, the
BL DR system uses an 8-tap PFB to divide the 1.5 GHz Nyquist range into
$N_{\text{coarse}}=512$ “coarse” spectral channels, which in turn are divided
among 8 compute nodes (MacMahon et al., 2018). This procedure performs a Fast
Fourier Transform (FFT) with a length of $P=2N_{\text{coarse}}=1024$. For
receivers with wide bandwidths, such as C-band at 3.95–8.00 GHz, multiple
copies of these elements, starting from the analog mixer, are employed to
cover the full band (NRAO, 2019).

The digital processing components of the BL DR system are done on custom
signal processing boards using field-programmable gate arrays (FPGAs),
provided by the Collaboration for Astronomy Signal Processing and Electronics
Research (CASPER; Hickish et al., 2016). These boards use fixed point
arithmetic and increase numerical bit size when doing computations (MacMahon
et al., 2018). Accordingly, both real and imaginary components of the
resulting complex voltages must be requantized (e.g. to $N_{\text{bits,r}}$)
before they are written to disk. The BL DR system records these as 8-bit
signed integers in GUPPI (Green Bank Ultimate Pulsar Processing Instrument;
DuPlain et al., 2008) raw format, based on FITS (Pence et al., 2010) and
stored as .raw files (Lebofsky et al., 2019).

Since raw voltage data comes at the highest resolution possible given the ADC
sampling rate, data volumes are large, especially during standard BL observing
campaigns. Therefore, we finely channelize or “reduce” raw data into
spectrograms (also known as dynamic spectra or “waterfall plots”), 2D arrays
of intensity (Stokes I) as a function of time and frequency (Lebofsky et al.,
2019). Multiple versions with different resolutions can be created from the
same set of raw data by varying the FFT length $N_{\text{fine}}$ and
integration factor $N_{\text{int}}$.

During fine channelization, an FFT of length $N_{\text{fine}}$ is performed on
complex raw voltages within individual coarse channels, resulting in
$N_{\text{fine}}$ fine channels each. So, we can express the full Nyquist
bandwidth as

$f_{N}=\frac{f_{s}}{2}=N_{\text{coarse}}N_{\text{fine}}\Delta f.$ (1)

This gives us an expression for the spectrogram’s frequency resolution:

$\Delta f=\frac{f_{s}/2}{N_{\text{coarse}}N_{\text{fine}}}.$ (2)

If the total observation length is $\tau$ and the number of time channels
(pixels) in the final spectrogram is $N_{t}$, then

$N_{t}=\frac{\tau}{\Delta t},$ (3)

assuming that $\tau$ is a multiple of the spectrogram’s time resolution
$\Delta t$. In practice, extraneous samples are truncated when necessary to
satisfy this requirement.

The integration factor $N_{\text{int}}$ is the number of spectra integrated in
the time direction. To get an expression for $\Delta t$, we can think in terms
of the total number of samples collected (for a single linear polarization):

$N_{s}=\tau f_{s}.$ (4)

The pipeline takes in $N_{s}$ real samples in time and, via a $P$-point FFT,
transforms the data into a complex 2D array in time-frequency space, with non-
integrated dimensions $N_{t}N_{\text{int}}\times PN_{\text{fine}}$.

$N_{s}=N_{t}N_{\text{int}}PN_{\text{fine}}=2N_{t}N_{\text{int}}N_{\text{coarse}}N_{\text{fine}}.$
(5)

Note that since the FFT is performed on real voltages, the unique frequency
extent is ultimately halved per the Nyquist range.

Combining Eqs. 2–5, we get

$\displaystyle N_{s}=\tau f_{s}$
$\displaystyle=2N_{t}N_{\text{int}}N_{\text{coarse}}N_{\text{fine}}$ (6)
$\displaystyle=2\frac{\tau}{\Delta
t}N_{\text{int}}N_{\text{coarse}}N_{\text{fine}}$ (7) $\displaystyle\Delta t$
$\displaystyle=\frac{2N_{\text{int}}N_{\text{coarse}}N_{\text{fine}}}{f_{s}}=\frac{N_{\text{int}}}{\Delta
f}.$ (8)

Although $N_{\text{fine}}$ and $N_{\text{int}}$ must both be integers, we
otherwise have fine control over $\Delta f$ and $\Delta t$ through Eqs. 2 and
8.

##  3 Code Methodology

As object-oriented software, setigen has a set of important classes and
routines that are described below. For more technical details and examples of
the API, see the full documentation777https://setigen.readthedocs.io/.

###  3.1 Spectrogram Module

The spectrogram module provides an interface for synthesizing Stokes I
(waterfall) data in a format common to radio SETI and is oriented around the
Frame class. A Frame object contains a 2D data array of intensities as a
function of time and frequency, as well as accompanying metadata, such as
starting frequency and time-frequency resolutions.

Data frames can be initialized from either saved observational data or frame
parameters. Frames can extract Stokes I data and observational metadata from
filterbank (.fil) or HDF5 files (.h5). The most important metadata for setigen
are the physical parameters of the underlying intensity data: resolutions and
ranges in both time and frequency. Empty frames can therefore be created
simply by specifying these parameters along with desired data array
dimensions.

####  3.1.1 Noise Synthesis

In most SETI applications, we search for statistically-significant signals
embedded in noise. Since voltage noise in the absence of RFI approximately
follows a zero-mean normal distribution (Thompson et al., 2017), the
radiometer noise in spectrogram data follows a chi-squared distribution
(McDonough & Whalen, 1995; Nita et al., 2007). When the time and frequency
resolutions are coarse enough, the spectrogram noise approaches a normal
distribution by the central limit theorem.

Specifically, suppose we have a sequence of input voltages $\\{x_{n}\\}$
following a Gaussian distribution with zero mean. During the coarse
channelization process, the polyphase filterbank applies, at its core, an FFT
to bring the voltages into frequency space:

$X_{k}=\sum_{n=0}^{N-1}{w_{n}x_{n}e^{-2\pi ikn/N}},\ k=0,\dots,N-1,$ (9)

where $N$ is the number of frequency bins and $\\{w_{n}\\}$ are coefficients
of a windowing function applied to improve the spectral response (Price,
2021).

More specifically, the filterbank sums over $M$ rows of $P$ samples before a
$P$-point FFT, so that the response of the $r$th row of $P$ samples is:

$X_{k,r}=\sum_{p=0}^{P-1}\left[\sum_{m=0}^{M-1}{w_{n^{\prime}}x_{n^{\prime\prime}}}\right]e^{-2\pi
ikp/P},$ (10)

where $n^{\prime}=mP+p$ and $n^{\prime\prime}=(r-M+m)P+p$ are indices of the
windowing coefficients and voltages samples in terms of $m$ and $p$. Here, we
assume that the $MP$ windowing coefficients are symmetric about the midpoint,
so that $w_{n}=w_{MP-n-1}$.

Ignoring quantization for the moment, we store the complex components of the
resulting FFT voltages, $\operatorname{Re}(X_{k})$ and
$\operatorname{Im}(X_{k})$, as raw voltage data. Since these are linear
combinations of independent zero-mean Gaussian variables (i.e. $x_{n}$), they
both follow zero-mean Gaussian distributions.

In the absence of a windowing function ($w_{n}=1$), for each channel besides
the real-valued DC and Nyquist bins, the variances of the real and imaginary
components are equal ($\sigma^{2}$; McDonough & Whalen, 1995). When a
windowing function is used, the underlying statistics can change such that the
variances of the complex components differ as a function of spectral bin (Nita
et al., 2007). However, for commonly chosen symmetrical windows (e.g.
Hamming), this effect is negligible in most spectral bins.

For a single linear polarization, the power is given by

$I_{x,k}=|X_{k}|^{2}=\operatorname{Re}(X_{k})^{2}+\operatorname{Im}(X_{k})^{2}$
(11)

Assuming both complex components have the same variance $\sigma^{2}$, the
power follows a chi-squared distribution with two degrees of freedom:

$I_{x,k}\sim\sigma^{2}\chi^{2}(2)$ (12)

During the fine channelization step, we integrate $N_{\text{int}}$ spectra in
the time direction and combine power from $N_{\text{pol}}$ polarizations.
Therefore, in the final Stokes I spectrogram, the total number of chi-squared
degrees of freedom is given by:

DOF $\displaystyle=2N_{\text{pol}}N_{\text{int}}=2N_{\text{pol}}\Delta f\Delta
t$ (13) $\displaystyle I_{k}$
$\displaystyle\sim\sigma^{2}\chi^{2}(2N_{\text{pol}}\Delta f\Delta t),$ (14)

using Eq. 8. For dual-polarization Stokes I data, $\text{DOF}=4\Delta f\Delta
t$. This allows us to generate synthetic chi-squared noise with the correct
number of degrees of freedom just from frame resolutions, which are either
directly specified or inferred from observations. Since non-calibrated
intensity values are arbitrarily scaled, we can simply scale the magnitudes of
synthetic chi-squared noise to match empirical observational noise
distributions.

The main function for noise synthesis across a frame is add_noise, which adds
random noise to every pixel in the data array. By default, it generates chi-
squared noise with a user-specified mean intensity $\mu$. Since the mean of a
chi-squared distribution equals the number of degrees of freedom, for dual-
polarization data, we have

$\displaystyle I_{k}$ $\displaystyle\sim\left(\frac{\mu}{4\Delta f\Delta
t}\right)\chi^{2}(4\Delta f\Delta t)$ (15) $\displaystyle\langle I_{k}\rangle$
$\displaystyle=\left(\frac{\mu}{4\Delta f\Delta t}\right)\cdot 4\Delta f\Delta
t=\mu$ (16) $\displaystyle\mathrm{Var}(I_{k})$
$\displaystyle=\left(\frac{\mu}{4\Delta f\Delta t}\right)^{2}\cdot 2\cdot
4\Delta f\Delta t=\frac{\mu^{2}}{2\Delta f\Delta t}.$ (17)

In addition to chi-squared noise, add_noise can also generate Gaussian noise.
By the central limit theorem, as the degrees of freedom increase, a chi-
squared distribution approaches a normal distribution. For example,
$N_{\text{int}}=51$ for BL’s standard high spectral resolution data product,
so $\text{DOF}=204$ and the resulting background noise is close to Gaussian.
Directly synthesizing Gaussian-distributed noise can save normalization steps
in data processing, but should be used carefully when comparing with real
observational data.

A useful extension of the noise synthesis function is add_noise_from_obs,
which draws from archived observational statistics to set realistic intensity
values. The observations were taken using the GBT at C-band and reduced to
(1.4 s, 1.4 Hz) resolution. For example, for chi-squared noise, the function
randomly selects an archived mean intensity, scales it to the appropriate
frame resolution, and populates noise per Eq. 15. An implementation detail of
BL’s fine channelization software,
rawspec888https://github.com/UCBerkeleySETI/rawspec, is that as part of the
FFT, intensity values are scaled up by a factor of the FFT length
$N_{\text{fine}}$. So, for observations going through the BL data pipeline
(i.e. the same digitization and coarse channelization hardware):

$\displaystyle\mu$ $\displaystyle\propto N_{\text{fine}}N_{\text{int}}$ (18)
$\displaystyle\propto N_{\text{fine}}\Delta f\Delta t$ (19)
$\displaystyle\propto\Delta t.$ (20)

Alternatively, the function also accepts user-provided arrays of background
noise intensity statistics from which to sample instead. This can be used for
synthesizing data with intensity ranges from other telescopes (e.g. Parkes) or
even GBT data at different frequency bands or sensitivities.

After noise synthesis, the frame will update class attributes storing the
estimated mean $\mu_{b}$ and standard deviation $\sigma_{b}$ of the background
noise. For an empty frame, the first noise synthesis function will set these
properties directly. For pre-loaded observational data and further noise
injection, the frame estimates the background noise through iterative sigma
clipping at the $3\sigma$ level to exclude outliers. For frames small enough
that noise statistics do not change over the frequency bandwidth, this enables
signal injection at desired SNR levels.

####  3.1.2 Signal Synthesis

For narrow-band signal synthesis, the add_signal function creates heuristic,
user-defined signals in spectrogram data. Our convention is that the
spectrogram data has time on the $y$-axis and frequency on the $x$-axis.

In spectrogram setigen, narrow-band signals have a “central” frequency at each
timestep and a unique spectral profile centered at that frequency. As such,
there are four main heuristic descriptors for a narrow-band signal in setigen:

  1. 1.

path – $I_{p}(t)$: Central signal frequencies as a function of time, e.g.
linear (constant) drift rate, quadratic drift rate

  2. 2.

t_profile – $I_{t}(t)$: Signal intensity as a function of time, e.g. constant
intensity, Gaussian pulses

  3. 3.

f_profile – $I_{f}(f,f_{0})$: Spectral profile as a function of frequency
(offset from central frequency), e.g. $\text{sinc}^{2}$ profile, Gaussian
profile

  4. 4.

bp_profile – $I_{bp}(f)$: Bandpass profile as a function of absolute frequency

These descriptors are parameters for add_signal and are Python functions by
type. A set of common functions are provided with setigen, and others can be
custom-written. The simplest and most ideal kind of narrow-band signal has a
constant intensity and drift rate; such signals can be created
straightforwardly through the wrapper function add_constant_signal.

For a pixel at $(t,f)$ in the time-frequency spectrogram, the intensity of a
synthetic signal is calculated as

$I(t,f)=I_{t}(t)I_{f}(f,I_{p}(t))I_{bp}(f).$ (21)

As such, Eq. 21 is computed for every pixel in the spectrogram, since there is
no robust way to constrain arbitrary intensity profiles. For example, even an
ideal Gaussian function is non-zero at all distances and defining a suitable
range depends on the experiment. For large spectrograms, it can be inefficient
to calculate intensities for pixels far from the main signal, so users can
provide a custom frequency range to limit the signal calculation.

The signal calculation is fully heuristic, in that the calculation is
completely user-specified and does not take other effects into account, such
as FFT leakage or spectral responses. Since intensity is treated as a function
of time and frequency, this process can overlook how intensities are
integrated in reality. As a partial solution, add_signal provides the option
to separately sub-integrate within each pixel in time and frequency
directions.

In a similar vein, a difficult effect to handle robustly is Doppler smearing,
in which a highly drifting signal will have its power spread into multiple
frequency channels within the same time channel (Sheikh et al., 2019). While
an analytical form exists for the spectral profile of a linearly drifting
cosine signal, the smearing effect will naturally apply to more complex
signals. Variable spectral profiles are not yet supported in setigen, but from
a user standpoint, it would be tedious to manually construct custom smearing
profiles that change at each timestep. Using a similar process to numerical
integration, add_signal has the option to approximate Doppler smearing by
computing and averaging a given number of copies of the signal, spaced evenly
between signal center frequencies in adjacent timesteps. For instance, for the
$i$th time channel at $t=t_{i}$, copies of the signal centered at even
spacings between $I_{p}(t_{i})$ and $I_{p}(t_{i+1})$ are averaged together to
get the $i$th spectral profile. This is done for all time channels, so that
channels with smaller signal drifts will be brighter than those with larger
signal drifts by the correct ratio, as long as the number of copies gives
enough coverage over the channel with the largest signal drift.

Sometimes it can be difficult or unwieldy to wrap up a desired signal property
into a separate function, or perhaps there is existing external code that
produces such properties. In these cases, we can instead use NumPy arrays to
describe these signals, rather than functions. As of now, the path, t_profile,
and bp_profile arguments can be arrays.

####  3.1.3 Common Frame Operations

Figure 1: Radio spectrogram plots created from setigen frames. A: Frame with
only synthetic chi-squared noise. B: Frame from panel A with an injected
synthetic signal at SNR=30. C: “Real” GBT observation of Voyager I carrier
signal at X-band. D: Frame from panel C with an injected synthetic signal at
SNR=1000, with the same drift rate as the injected signal in panel B.

Besides supporting noise and narrow-band signal injection, setigen comes with
a set of tools for radio spectrogram analysis. These range from convenience
functions for parameter calculations to frame-level data transformations.

For instance, estimating the SNR of a signal in an integrated spectrum is a
common step in radio analysis. This can be done through a frame’s integrate
function, which can also be used along the frequency axis to produce an
intensity time series array.

To inject a signal at a desired SNR, the get_intensity function calculates the
requisite signal level as

$I_{t}=\text{SNR}\cdot\frac{\sigma_{b}}{N_{t}^{1/2}},$ (22)

assuming that the frame has background noise with standard deviation
$\sigma_{b}$ and that the SNR is measured by dividing the integrated signal
maximum by the integrated noise deviation. As discussed in Section 3.1.1, each
frame tracks an estimate of $\sigma_{b}$ calculated using iterative sigma
clipping and updates it when synthetic noise is injected.

It can be convenient to define signals in terms of the pixels they traverse
rather than the frequencies. To convert between these for a given frame, one
can use the get_frequency and get_index functions. We define the unit drift
rate for a given spectrogram resolution to be the drift rate given by

$\dot{\nu}_{1}=\frac{\Delta f}{\Delta t},$ (23)

which can be accessed with the unit_drift_rate attribute. For a linearly-
drifting signal passing through the top and bottom of the frame, the
corresponding drift rate can be calculated using the get_drift_rate function.

Given a frame with a linearly-drifting signal, we can “de-drift” the frame
using setigen.dedrift. This shifts each spectrum an appropriate amount along
the frequency direction so that such a signal would, on average, appear to
have zero frequency drift, making it simpler to calculate the SNR. In
practice, empirical drift rates are not generally multiples of the unit drift
rate, so de-drifted signals will not be perfectly aligned.

We can create a “slice” of a frame by specifying left and right frequency
indices, analogous to NumPy array slicing, by using the frame’s get_slice
function. This results in a new frame with a truncated range, which can be
helpful for isolating signals in time-frequency space for further analysis.

If one is interfacing with other BL or astronomy codebases, outputting setigen
frames to filterbank or HDF5 format can be very useful. These are done via the
save_fil and save_hdf5 functions. Frame objects can also be written and loaded
with pickle, a convenient serialization method that can keep data and user-
provided metadata together.

####  3.1.4 Demonstration: Spectrogram Module

We present a minimal working example of creating a data frame with synthetic
noise and a drifting signal. First, we construct an empty frame with the
desired resolution; here, we use parameters that match those of BL’s high
frequency resolution data product:

⬇

from astropy import units as u

import setigen as stg

\parframe = stg.Frame(fchans=256,

tchans=16,

df=2.7939677238464355*u.Hz,

dt=18.253611008*u.s,

fch1=6095.214842353016*u.MHz)

Then, we add chi-squared noise with a desired mean, such as 10:

⬇

frame.add_noise(x_mean=10, noise_type=’chi2’)

Finally, we add a simple drifting signal through our frame at SNR=30 and plot
the result in decibels (dB). The inputs to add_signal shown below are pre-
written library functions that themselves return the functions described in
Section 3.1.2. Since they are indeed Python functions by type, the signal
parameters allow for much more flexibility beyond this basic example.

⬇

frame.add_signal(

stg.constant_path(

f_start=frame.get_frequency(index=100),

drift_rate=2*u.Hz/u.s

),

stg.constant_t_profile(

level=frame.get_intensity(snr=30)

),

stg.gaussian_f_profile(width=10*u.Hz),

stg.constant_bp_profile(level=1)

)

\parframe.bl_plot()

The frames after adding noise and after adding the signal are shown in Figures
1A and 1B.

We also show an example with a signal detected from Voyager I in an X-band
observation using the GBT, in Figure 1C. Injecting a signal into the Voyager
frame with the same drift rate as in the example (Figure 1B), now at SNR=1000,
we get Figure 1D.

###  3.2 Raw Voltage Module

Figure 2: Basic layout of a voltage pipeline written using setigen.voltage.

The raw voltage module is designed for synthesizing complex voltage data,
providing a set of classes that models the signal processing pipeline
described in Section 2. Instead of directly synthesizing spectrogram data, we
can produce real voltages, pass them through a virtual pipeline, and record to
file in GUPPI raw format. As this process models actual hardware used by BL
for recording raw voltages, this enables lower level testing and
experimentation.

The basic signal flow is shown in Figure 2. At the lowest level, a DataStream
can accept noise and signal sources (as Python functions) and generate real
voltages on demand. An Antenna models an antenna or dish used in radio
telescopes and has one or two DataStream objects, corresponding to linear
polarizations that are unique and not necessarily correlated. As described in
Section 2, the sampled real voltages are passed to a processing pipeline which
consists, at its core, of a digitizer, a polyphase filterbank (PFB), and a
requantizer. In hardware, processing is done in fixed point arithmetic on an
FPGA, but for simplicity, we use floating point. The digitizer quantizes input
voltages to a specified number of bits and a target full width at half maximum
(FWHM) in the quantized voltage space. The filterbank implements a software
PFB, coarsely channelizing input voltages. The requantizer takes the resulting
complex voltages and quantizes each component to either 8 or 4 bits, suitable
for saving into GUPPI raw format.

The RawVoltageBackend object wraps around these elements and connects the full
pipeline together. Given an observation length in seconds or a number of data
recording “blocks,” the main function record retrieves real voltage samples as
needed and passes them through each backend element, finally saving the
quantized complex voltages out to disk.

Since voltage data is taken with very high sample rates, e.g. Gigasamples/sec
(Gsps), the voltage module is much more computationally expensive than the
spectrogram module. To increase efficiency, most of the data manipulations are
done with matrix operations, allowing for GPU acceleration with CuPy (Okuta et
al., 2017).

####  3.2.1 Antennas and DataStreams

The DataStream class represents a stream of real voltage data for a single
polarization and antenna. A data stream has an associated sample rate $f_{s}$,
such as 3 GHz for the BL DR. As of now, the voltage module does not implement
heterodyne mixing or bandpass filtering. Instead, data streams use a reference
frequency fch1 and frequency sign (ascending or descending from fch1) for
voltage calculations.

The Antenna class is similarly defined by a sample rate, reference frequency,
and frequency sign. For two linear polarizations, an Antenna’s data streams
are available via the x and y attributes. For one polarization, only the
former is available. For convenience, the streams attribute gets the list of
available data streams for an antenna. One can add noise and signal sources to
these individual data streams.

Real voltage noise is modeled as ideal Gaussian noise and added through the
add_noise function. Note that this actually stores a Python function to the
data stream that is only evaluated when get_samples is called. It also updates
the data stream’s noise_std attribute, which keeps track of the standard
deviation of the voltages in that data stream. This is useful for injecting
signals at target spectrogram SNRs.

Drifting cosine signals can be added to a data stream using
add_constant_signal. For more complex signals, one can write custom voltage
functions to add using add_signal. Voltage signal sources are Python functions
that accept an array of timestamps and output a corresponding sequence of real
voltages. Here is a simple example that adds a non-drifting cosine signal with
frequency f_start:

⬇

def cosine_signal(ts):

delta_f = f_start - antenna.x.fch1

return np.cos(2 * np.pi * delta_f * ts)

\parantenna.x.add_signal(cosine_signal)

As custom signals are added, the noise_std parameter may no longer accurately
reflect the background noise. In these cases, one can run the data stream’s
update_noise function to estimate noise empirically. This is not done by
default to save computation, especially when there are multiple well-behaved
voltage sources (e.g. Gaussian noise, cosine signals).

####  3.2.2 Quantization

The quantization process takes a continuous input voltage distribution and
scales it to a target distribution that can be described by $N_{\text{bits}}$
bits. Since real voltage noise can be modeled by a Gaussian process, we can
define this scaling in terms of the standard deviation or FWHM.

For real voltages $\\{v\\}$, target bit size $N_{\text{bits}}$, target mean
$\mu_{q}$ (ideally 0), and target standard deviation $\sigma_{q}$, the
quantized voltages $v_{q}$ are given by:

$\displaystyle v_{s}$
$\displaystyle=\left\lfloor\frac{\sigma_{q}}{\sigma_{v}}(v-\langle
v\rangle)+\mu_{q}\right\rfloor$ (24) $\displaystyle v_{q}$
$\displaystyle=\min(\max(-2^{N_{\text{bits}}-1},v_{s}),2^{N_{\text{bits}}-1}-1)$
(25)

We can define quantizers in terms of a target FWHM $w_{q}$, in which case
$\sigma_{q}=\frac{w_{q}}{2\sqrt{2\ln 2}}$.

The digitizer quantizes real voltages, while the requantizer receives complex
voltages and quantizes per complex component. Quantization is run per
polarization and antenna, and background statistics can be cached to save
computation in subsequent calls. This is facilitated by the RealQuantizer and
ComplexQuantizer classes.

####  3.2.3 Polyphase Filterbank

The PolyphaseFilterbank class implements and applies a PFB to quantized input
voltages. Instead of directly applying a $P$-point FFT, a PFB first splits
incoming voltages between $P$ branches and lets $M$ samples accumulate in each
branch (Price, 2021). A windowing function is applied over the $M\times P$
samples, the samples are summed over the $M$ so-called polyphase taps, and
finally a $P$-point FFT is taken of the result to get complex raw voltages in
$N_{\text{coarse}}=P/2$ coarse channels. Further samples are read in groups of
$P$ and split between the PFB branches; accumulated samples step forward to
the next tap to make room. PFBs have a better channel response than standard
FFTs, especially as $M$ increases, and are common in high spectral resolution
radio backends (Price, 2021).

The two main parameters for a PolyphaseFilterbank are the number of taps $M$
and the number of branches $P$. Since the PFB works on $MP$ samples at once,
the object continuously caches samples for on-demand computation. The PFB also
accepts a symmetric windowing function as an argument (Hamming, by default)
and generates $MP$ coefficients up front (Blackman & Tukey, 1958).

####  3.2.4 Combining Components and Recording Data

The RawVoltageBackend class contains the full machinery to collect, process,
and write complex voltage data to GUPPI raw files, as in the standard pipeline
shown in Figure 2. Nevertheless, since the individual signal processing
components are all exposed as part of the voltage module, custom pipelines can
be written by chaining them in different ways.

A RawVoltageBackend takes in components external to the data recording process
as parameters, such as the antenna, digitizer, PFB, and requantizer. All other
parameters and functions are specific to data recording and actually obtaining
data from the external components.

As described by Lebofsky et al. (2019), the block size $N_{\text{blocksize}}$
refers to the number of bytes in a single block of data in GUPPI format. Each
data block has an associated header with observing metadata, such as target
and frequency information. The number of blocks per file also must be
specified to size individual raw files; multiple raw files may be associated
with a single pointing. For standard 5 minute GBT observations, BL DR uses
$N_{\text{blocksize}}=134217728$ with 128 blocks per file.

To specify the coarse channels that are actually recorded to disk, we can set
the starting index and the number of consecutive channels $N_{\text{chan}}$ to
ultimately save. Purely for computational efficiency, we always perform a full
FFT and truncate to obtain the desired coarse channels, instead of directly
doing the transform operation on the subset of coarse channels. Especially
when using a GPU to accelerate synthesis, this can fill up memory rather
quickly, potentially to the point of overflow. Therefore, the
RawVoltageBackend has an additional option to divide individual data blocks
into a given number of sub-blocks, such that each sub-block will fully fit in
memory.

For a single antenna, the number of bytes $N_{\text{blocksize}}$ in a block
can be related to the number of time channels $N_{t,\text{block}}$
corresponding to a single block in (non-integrated) spectrogram format as

$\displaystyle N_{\text{blocksize}}$
$\displaystyle=2N_{\text{pol}}\left(\frac{N_{\text{bits,r}}}{8}\right)N_{\text{chan}}N_{t,\text{block}}$
(26)
$\displaystyle=\frac{1}{4}N_{\text{pol}}N_{\text{bits,r}}N_{\text{chan}}N_{t,\text{block}},$
(27)

based on the structure of raw files as described by Lebofsky et al. (2019).

####  3.2.5 Multi-Antenna Support

To simulate voltage data for interferometric pipelines, it can be useful to
synthesize raw voltage data from multiple antennas. setigen supports
synthesizing multi-antenna output through the MultiAntennaArray class, which
creates a list of $N_{\texttt{ant}}$ antennas each with an associated integer
delay (in time samples). In addition to the individual data streams that allow
the user to add noise and signals to each antenna, there are “background” data
streams bg_x and bg_y in MultiAntennaArray, representing correlated noise or
RFI that is detected at each antenna, subject to the (relative) delays.
Signals and noise can therefore be added to the background across all array
elements as well as to individual antennas.

The only difference in the pipeline is instead of supplying a Antenna as input
to a RawVoltageBackend, one would supply a MultiAntennaArray. Then, the output
is saved as a multi-antenna extension of the GUPPI raw format.

####  3.2.6 Creating Signals at a Target Spectrogram SNR

During the course of the full signal processing pipeline, an injected cosine
signal passes through multiple quantization and FFT steps. In many SETI
experiments, a signal’s SNR in spectrogram data is used for thresholding and
analysis, so it is important to be able to estimate this SNR given pipeline
parameters.

Suppose that we have a cosine signal with amplitude $A$ at a frequency
corresponding to the center of a fine spectral channel, and that this signal
is injected onto a background of Gaussian noise
$\mathcal{N}(0,\sigma_{v}^{2})$. Since the voltage data is real-valued, the
signal magnitude becomes $A/2$ in frequency space. As the voltages pass
through the coarse and fine channelization steps, the signal magnitude picks
up factors of $P$ and $N_{\text{fine}}$, respectively, compared to the
background noise.

The background noise will follow a chi-squared distribution with
$\text{DOF}=2N_{\text{pol}}N_{\text{int}}$ (Section 3.1.1), scaled by
multiplicative factors arising from quantization and FFT calculations. Since
the input voltage noise has variance $\sigma_{v}^{2}$, the standard deviation
of the noise power $\sigma_{b}$ will be proportional to the standard deviation
$\sigma_{b,0}$ of a chi-squared distribution with mean $\sigma_{v}^{2}$. The
time integration step to get the SNR will reduce this noise by a factor of
$N_{t}^{1/2}$.

To get an expression for $N_{t}$ given observation parameters, suppose our
synthetic observation has $N_{\text{block}}$ total blocks and that the time
covered by a single block is $\tau_{\text{block}}$. Then, we have the
following equations:

$\displaystyle\Delta t$ $\displaystyle=\frac{N_{\text{int}}}{\Delta
f}=\frac{P}{f_{s}}N_{\text{fine}}N_{\text{int}}$ (28)
$\displaystyle\tau_{\text{block}}$ $\displaystyle=N_{t,\text{block}}\Delta t$
(29) $\displaystyle N_{t}$
$\displaystyle=\frac{N_{\text{block}}\tau_{\text{block}}}{N_{\text{int}}\Delta
t}=\frac{N_{\text{block}}N_{t,\text{block}}}{N_{\text{int}}}.$ (30)

Combining all of these factors, we can express the final SNR of the signal as
the ratio between the integrated (mean) signal power and the integrated
background noise standard deviation as

$\displaystyle\sigma_{b,0}$
$\displaystyle=\sigma_{v}^{2}\left(\frac{2}{\text{DOF}}\right)^{1/2}$ (31) SNR
$\displaystyle=\frac{I}{\sigma_{b}}=\frac{(A/2)^{2}PN_{\text{fine}}}{\sigma_{b,0}/N_{t}^{1/2}}.$
(32)

This yields the amplitude or signal level in terms of the target SNR:

$A=\left(\text{SNR}\cdot\frac{4\sigma_{b,0}}{PN_{\text{fine}}N_{t}^{1/2}}\right)^{1/2}$
(33)

Notice that $A$ has a linear dependence on the standard deviation $\sigma_{v}$
of the real voltage noise in a data stream, which can arise from multiple
sources, especially in a multi-antenna array. Given pipeline parameters, the
get_level function can be used to calculate $A/\sigma_{v}$.

For a non-drifting cosine signal, we can also approximate the effect of
spectral leakage between fine channels by comparing the signal frequency to
the nearest channel central frequency. A signal with amplitude $A$ centered at
a frequency $\delta f$ away from the center of the closest fine spectral
channel will have its power $I$ attenuated by999$\text{sinc }x=\sin{x}/x$

$\frac{I^{\prime}}{I}=\text{sinc}^{2}\left(\frac{|\delta f|}{\Delta
f}\right).$ (34)

Since intensity goes as voltage squared, we provide a function
get_leakage_factor to calculate an amplitude adjustment factor $f_{l}$ to
easily scale from $A$ to a new amplitude $A^{\prime}$ that corresponds to the
non-attenuated intensity:

$\displaystyle f_{l}$ $\displaystyle=\frac{1}{\text{sinc}\left(\frac{|\delta
f|}{\Delta f}\right)}$ (35) $\displaystyle A^{\prime}$ $\displaystyle=f_{l}A.$
(36)

Finally, for a linearly-drifting cosine signal, if the drift rate $\dot{\nu}$
exceeds the unit drift rate $\dot{\nu}_{1}$, signal power will be smeared
across multiple frequency bins in spectrogram data. This is a linear effect in
spectrogram data, so cosine amplitudes should be increased by a factor of
$(\dot{\nu}/\dot{\nu}_{1})^{1/2}$ to counter-act the apparent loss in power.

####  3.2.7 Injecting Synthetic Signals into Raw Voltage Data

In addition to creating fully synthetic complex voltage data from scratch, the
RawVoltageBackend supports injecting or adding synthetic data to existing
observational GUPPI raw data. The pipeline remains mostly the same, except for
a few important differences that we detail below.

In order to get meaningful results, we must know and match details about the
specific signal processing pipeline that produced the existing raw data.
setigen provides a helper function called get_raw_params to extract header
information from the raw data file, but other information must be provided
separately by the user, such as the sampling rate and PFB parameters.

Since recorded voltage data has already gone through multiple quantization
steps, we cannot directly add time series voltages together (i.e. at the
original ADC sampling rate). Instead, we choose to synthesize complex voltage
data separately, add it to the recorded voltage data, and apply a final
quantization step to match the initial distribution as best as possible.

However, this process requires that we create and process signals that are not
necessarily embedded in noise. In typical narrow-band signal injection
scenarios, we wish to synthesize and inject signals whose distributions are
non-Gaussian (e.g. a cosine signal). Since the quantization steps assume that
the input and output voltage distributions are both Gaussian, attempting to
quantize bare narrow-band signals will cause distortion and introduce clipping
artifacts. Furthermore, without a reference noise distribution, quantization
can scale the magnitude of processed signals in undesired ways, making SNR
estimation difficult.

To address these issues, we approach the quantization steps differently. If
there is already a synthetic noise source, we proceed normally through all
steps in the pipeline. Otherwise, we skip the initial digitization step before
the PFB, and instead treat the input voltages as if they followed a zero-mean
Gaussian distribution with variance 1. Using a reference distribution allows
us to set signal magnitudes with the get_level function to achieve target SNR
levels. We then estimate the post-PFB mean and standard deviation of the
reference Gaussian voltages and quantize the synthetic voltages based on these
values instead of those from the “real” synthetic distribution. This way, if
the synthesized voltages were actually embedded in $\mathcal{N}(0,1)$ noise,
the resulting signal quantization would be very similar.

For each data block in the recorded raw file, the RawVoltageBackend will set
requantizer statistics (target mean $\mu_{q}$ and target standard deviation
$\sigma_{q}$) calculated from the existing data for each combination of
antenna, polarization, and complex component. The synthetic voltages are
requantized to the corresponding standard deviations in each complex
component, but instead of centering to the target mean, they are centered to
zero mean. This is so that when we add the quantized synthetic data to the
existing data, we do not change the overage voltage mean. After these are
added together, we finally requantize once more to the target mean and target
standard deviation to match the existing data statistics and magnitudes as
best as possible.

####  3.2.8 Demonstration: Voltage Module

Figure 3: Spectrogram derived from synthetic raw voltages, showing the edge of
the coarse channel bandpass shape and a bright, slightly drifting cosine
signal. The top panel shows an integrated profile, showing PFB scalloping loss
towards the left and the synthetic signal towards the right.

Here, we present a simple pipeline created with the raw voltage module to
inject a drifting cosine signal in Gaussian noise. First, we create the signal
processing elements:

⬇

from astropy import units as u

from setigen.voltage import *

\pard = RealQuantizer(target_fwhm=32,

num_bits=8)

\parf = PolyphaseFilterbank(num_taps=8,

num_branches=1024)

\parr = ComplexQuantizer(target_fwhm=32,

num_bits=8)

Then, we create the antenna, setting the sampling rate and reference
frequency. With two polarizations, we can add Gaussian noise and a constant
amplitude, Doppler drifting cosine signal to both data streams:

⬇

a = Antenna(sample_rate=3*u.GHz,

fch1=6000*u.MHz,

ascending=True,

num_pols=2)

\parfor s in a.streams:

s.add_noise(v_mean=0,

v_std=1)

\pars.add_constant_signal(f_start=6002.1*u.MHz,

drift_rate=-2*u.Hz/u.s,

level=0.004)

We connect these components through the recording backend, defining the
dimensions and size of the final raw voltage data product, and record a block
of data to file.

⬇

rvb = RawVoltageBackend(a,

digitizer=d,

filterbank=f,

requantizer=r,

start_chan=0,

num_chans=64,

block_size=134217728,

blocks_per_file=128,

num_subblocks=32)

\parrvb.record(output_file_stem=’example_1block’,

num_blocks=1,

length_mode=’num_blocks’,

header_dict={’TELESCOP’: ’GBT’},

verbose=True)

After saving the raw voltages to disk, we reduce using rawspec with
$N_{\text{fine}}=1024$ and $N_{\text{int}}=4$. A snippet of the resulting
spectrogram output is shown in Figure 3, where intensities are plotted on a
decibel scale. The signal is readily apparent, as is the frequency bandpass
shape arising from the PFB.

##  4 Discussion

###  4.1 Limitations

While setigen is a flexible library that enables quick narrow-band dataset
generation, it is important to discuss the limitations when using it for
science.

First and foremost, setigen relies on heuristic, user-defined signals, rather
than simulations from first principles. The search for technosignatures is
necessarily informed by human bias, specifically applied via our assumptions
about a technosignature’s potential characteristics and morphology. It is
possible that radiation from an extraterrestrial intelligence will exist in a
form that we have not considered or designed searches for. Even when we
consider only the problem of excision of anthropogenic RFI, we have to be
careful when applying algorithms developed using the simplest of narrow-band
signals. Although there might never be a way of fully emulating the breadth
and variety of the RFI environment, setigen can still be used to generate
labeled, complex signals to test the efficacy of new and existing algorithms.

In a similar vein, the spectrogram module enables users to quickly generate
signals that “look” like the narrow-band signals we see in observations.
However, since spectrogram signal injection does not have access to phase
information, it is impossible to replicate the “correct” intensity statistics
when adding a signal to integrated Stokes I noise. For example, adding a
perfect cosine signal to zero-mean Gaussian noise in the voltage domain
results in a non-central chi-squared intensity distribution in Stokes I data,
but adding a signal with constant intensity directly to chi-squared noise in a
spectrogram does not result in the same distribution (over the pixels occupied
by that signal; McDonough & Whalen, 1995). While this effect is negligible for
high SNR signals, algorithms developed to target low SNR signals may suffer
from intrinsic inaccuracies in the intensity statistics.

Signal injection in the complex voltage domain also has limitations since we
are not able, in software, to directly add signals in the real (analog)
voltage stage. Raw data is quantized multiple times in hardware, so the
injection step has to take place using complex voltages that are quantized in
a similar way. While fundamental steps in the pipeline are linear, such as PFB
operations (Eq. 10), quantization inherently breaks this linearity. Because of
this, summing real and synthetic voltages that are independently processed can
lead to artifacts and intensity discrepancies that would not arise if we could
inject at the start of the signal processing pipeline.

###  4.2 Future Directions

setigen is written and developed with the needs of SETI researchers in mind,
so new functionality and improvements are constantly being added. Here, we
describe some potential enhancements that may be added in the near future.

As it stands, the spectrogram module is especially targeted at producing small
frames with synthetic signals rather than injecting into large, broadband
observations. While this suffices in many cases, it may be useful to inject
within large data files in which frequency bandpass shapes significantly
change the background intensities. For instance, for use in SNR estimation,
setigen calculates background noise statistics over an entire frame rather
than localized around the target signal injection frequency. For a large
enough frame, this is both an inefficient and inaccurate calculation due to
variable bandpass shapes. An improvement would be to localize the noise
calculation to a window around the target injection site, as well as to
similarly localize the signal injection calculation to prevent unnecessary
computation.

The spectrogram module is also currently designed expressly to synthesize
narrow-band signals. There are many similarities in both signal processing and
experimental design between technosignature searches and searches for time-
varying phenomena such as pulsars and fast radio bursts (FRBs); setigen could
thus be expanded to include broadband signal injection (Zhang et al., 2018;
Gajjar et al., 2021).

An exciting potential addition is to use parameterized ML methods to create
labeled, realistic signals. By taking ideas from style transfer, a synthetic
RFI signal could be created by specifying heuristic parameters and having an
ML model generate such a signal with RFI-like properties (Gatys et al., 2016;
Dai et al., 2017). While generative adversarial networks (GANs) have been used
before to create radio spectrograms (Zhang et al., 2018), conditional GANs
that accept input parameters might help produce more specific, labeled
signals, which can be better for certain SETI experiments. Furthermore, better
RFI modeling could help improve ML-based searches for astrophysical phenomena
like FRBs in the presence of different classes of RFI.

Some of these enhancements may use a lot more computational power than the
current synthesis process, so the option to GPU-accelerate the standard
spectrogram module would be critical. Some of these enhancements may require a
more careful look at file input/output methods when reading and writing large
observational data files to avoid unnecessary or slow operations.

The raw voltage module can also be expanded to support alternate radio
telescope configurations and backends, such as those behind interferometers
like MeerKAT (Jonas, 2009). While setigen already has basic multi-antenna
functionality, it could be helpful to build on this with general-use
utilities, such as routines that predict how a given injected signal would
appear across multiple antennas or beams. The voltage module could also
support additional requantization and recording modes, such as 2 and 16-bit.
As interferometer usage in modern radio SETI continues to grow, setigen
capabilities can be extended to help test signal detection in commensal and
beam-formed observations (Czech et al., 2021).

##  5 Summary

In this paper, we presented setigen, an open-source Python library for the
creation and injection of synthetic narrow-band radio signals. setigen can
produce both finely channelized spectrogram data and coarsely channelized
complex voltage data. The spectrogram module is designed to be intuitive and
quick to use to facilitate the construction of synthetic datasets for SETI
experiments and testing. While the voltage module is more complex and
computationally intensive, it enables analysis of signals that pass through a
software-defined pipeline, which can be helpful in understanding the
implications of the instrumentation pipeline itself in SETI searches.

setigen is constantly being improved with the needs of SETI research in mind.
As open-source software, the library is freely available, and we encourage the
SETI community to use and contribute to it.

##  6 Acknowledgements

Breakthrough Listen is managed by the Breakthrough Initiatives, sponsored by
the Breakthrough Prize Foundation. The Green Bank Observatory is a facility of
the National Science Foundation, operated under cooperative agreement by
Associated Universities, Inc. We thank the staff at the Green Bank Observatory
for their operational support.

## References

  * Bellanger et al. (1976) Bellanger, M., Bonnerot, G., & Coudreuse, M. 1976, IEEE Transactions on Acoustics, Speech, and Signal Processing, 24, 109
  * Blackman & Tukey (1958) Blackman, R. B., & Tukey, J. W. 1958, Bell System Technical Journal, 37, 485
  * Brzycki et al. (2020) Brzycki, B., Siemion, A. P., Croft, S., et al. 2020, Publications of the Astronomical Society of the Pacific, 132, 114501
  * Collette et al. (2017) Collette, A., Tocknell, J., Caswell, T. A., et al. 2017
  * Czech et al. (2021) Czech, D., Isaacson, H., Pearce, L., et al. 2021, PASP, 133, 064502
  * Dai et al. (2017) Dai, B., Fidler, S., Urtasun, R., & Lin, D. 2017, in Proceedings of the IEEE International Conference on Computer Vision, 2970–2979
  * Drake (1961) Drake, F. D. 1961, Physics Today, 14, 40.  http://adsabs.harvard.edu/abs/1961PhT....14...40D
  * DuPlain et al. (2008) DuPlain, R., Ransom, S., Demorest, P., et al. 2008, in Advanced Software and Control for Astronomy II, Vol. 7019, International Society for Optics and Photonics, 70191D
  * Enriquez & Price (2019) Enriquez, E., & Price, D. 2019, ascl, ascl
  * Enriquez et al. (2017) Enriquez, J. E., Siemion, A., Foster, G., et al. 2017, ApJ, 849, 104
  * Gajjar et al. (2021) Gajjar, V., Perez, K. I., Siemion, A. P. V., et al. 2021, AJ, 162, 33
  * Gatys et al. (2016) Gatys, L. A., Ecker, A. S., & Bethge, M. 2016, in Proceedings of the IEEE conference on computer vision and pattern recognition, 2414–2423
  * Harp et al. (2019) Harp, G., Richards, J., Shostak, S., et al. 2019, arXiv preprint arXiv:1902.02426
  * Harris & Haines (2011) Harris, C., & Haines, K. 2011, Publications of the Astronomical Society of Australia, 28, 317
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hickish et al. (2016) Hickish, J., Abdurashidova, Z., Ali, Z., et al. 2016, Journal of Astronomical Instrumentation, 5, 1641001
  * Hunter (2007) Hunter, J. D. 2007, Computing in science & engineering, 9, 90
  * Jonas (2009) Jonas, J. L. 2009, Proceedings of the IEEE, 97, 1522
  * Lebofsky et al. (2019) Lebofsky, M., Croft, S., Siemion, A. P. V., et al. 2019, PASP, 131, 124505
  * MacMahon et al. (2018) MacMahon, D. H. E., Price, D. C., Lebofsky, M., et al. 2018, Publications of the Astronomical Society of the Pacific, 130, 044502.  https://doi.org/10.1088%2F1538-3873%2Faa80d2
  * Margot et al. (2018) Margot, J.-L., Greenberg, A. H., Pinchuk, P., et al. 2018, AJ, 155, 209
  * Margot et al. (2021) Margot, J.-L., Pinchuk, P., Geil, R., et al. 2021, The Astronomical Journal, 161, 55
  * McDonough & Whalen (1995) McDonough, R. N., & Whalen, A. D. 1995, Detection of signals in noise (Academic Press)
  * Nita et al. (2007) Nita, G. M., Gary, D. E., Liu, Z., Hurford, G. J., & White, S. M. 2007, Publications of the Astronomical Society of the Pacific, 119, 805
  * NRAO (2019) NRAO. 2019, GBT Proposer’s Guide, ,
  * Okuta et al. (2017) Okuta, R., Unno, Y., Nishino, D., Hido, S., & Loomis, C. 2017, in Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS).  http://learningsys.org/nips17/assets/papers/paper_16.pdf
  * Oliphant (2006) Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (Trelgol Publishing USA)
  * Pence et al. (2010) Pence, W. D., Chiappetti, L., Page, C. G., Shaw, R. A., & Stobie, E. 2010, Astronomy & Astrophysics, 524, A42
  * Pinchuk & Margot (2021) Pinchuk, P., & Margot, J.-L. 2021, arXiv preprint arXiv:2108.00559
  * Pinchuk et al. (2019) Pinchuk, P., Margot, J.-L., Greenberg, A. H., et al. 2019, The Astronomical Journal, 157, 122
  * Pinchuk et al. (2019) Pinchuk, P., Margot, J.-L., Greenberg, A. H., et al. 2019, AJ, 157, 122
  * Price et al. (2019) Price, D., Enriquez, J., Chen, Y., & Siebert, M. 2019, The Journal of Open Source Software, 4, 1554
  * Price (2021) Price, D. C. 2021, in The WSPC Handbook of Astronomical Instrumentation: Volume 1: Radio Astronomical Instrumentation (World Scientific), 159–179
  * Price et al. (2018) Price, D. C., MacMahon, D. H. E., Lebofsky, M., et al. 2018, Publications of the Astronomical Society of Australia, 35, 41
  * Price et al. (2020) Price, D. C., Enriquez, J. E., Brzycki, B., et al. 2020, AJ, 159, 86
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, Astronomy & Astrophysics, 558, A33
  * Sheikh et al. (2020) Sheikh, S. Z., Siemion, A., Enriquez, J. E., et al. 2020, arXiv preprint arXiv:2002.06162
  * Sheikh et al. (2019) Sheikh, S. Z., Wright, J. T., Siemion, A., & Enriquez, J. E. 2019, The Astrophysical Journal, 884, 14
  * Siemion et al. (2013) Siemion, A. P. V., Demorest, P., Korpela, E., et al. 2013, The Astrophysical Journal, 767, 94
  * Tarter (2001) Tarter, J. 2001, Annual Review of Astronomy and Astrophysics, 39, 511
  * Taylor (1974) Taylor, J. H. 1974, A&AS, 15, 367
  * Thompson et al. (2017) Thompson, A. R., Moran, J. M., & Swenson, George W., J. 2017, Interferometry and Synthesis in Radio Astronomy, 3rd Edition, doi:10.1007/978-3-319-44431-4
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature methods, 17, 261
  * Welch et al. (2009) Welch, J., Backer, D., Blitz, L., et al. 2009, Proceedings of the IEEE, 97, 1438
  * Werthimer et al. (1985) Werthimer, D., Tarter, J., & Bowyer, S. 1985, in Symposium-International Astronomical Union, Vol. 112, Cambridge University Press, 421–424
  * Worden et al. (2017) Worden, S. P., Drew, J., Siemion, A., et al. 2017, Acta Astronautica, 139, 98 .  http://www.sciencedirect.com/science/article/pii/S0094576517303144
  * Wright et al. (2014) Wright, S. A., Werthimer, D., Treffers, R. R., et al. 2014, in Ground-based and Airborne Instrumentation for Astronomy V, Vol. 9147, International Society for Optics and Photonics, 91470J
  * Zhang et al. (2018) Zhang, Y. G., Gajjar, V., Foster, G., et al. 2018, ApJ, 866, 149
  * Zhang et al. (2018) Zhang, Y. G., Won, K. H., Son, S. W., Siemion, A., & Croft, S. 2018, in 2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP), IEEE, 1114–1118
  * Zhang et al. (2019) Zhang, Y. G., Won, K. H., Son, S. W., Siemion, A., & Croft, S. 2019, arXiv e-prints, arXiv:1901.04636


# PKS 1424+240: yet another masquerading BL Lac object as a possible IceCube
neutrino source

P. Padovani1,2, B. Boccardi3, R. Falomo4, P. Giommi5,6,7
1European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching
bei München, Germany
2Associated to INAF - Osservatorio di Astrofisica e Scienza dello Spazio, Via
Piero Gobetti 93/3, I-40129 Bologna, Italy
3Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn,
Germany
4INAF - Osservatorio Astronomico di Padova, vicolo dell’Osservatorio 5,
I-35122, Padova, Italy
5Associated to Agenzia Spaziale Italiana, ASI, via del Politecnico s.n.c.,
I-00133 Roma, Italy
6Institute for Advanced Study, Technische Universität München,
Lichtenbergstrasse 2a, D-85748 Garching bei München, Germany
7Center for Astro, Particle and Planetary Physics, New York University Abu
Dhabi, United Arab Emirates
E-mail<EMAIL_ADDRESS>
(Accepted XX. Received YY; in original form ZZ)

###### Abstract

We show that the blazar PKS 1424+240, which has been recently associated by
IceCube with a neutrino excess at the $3.3\,\sigma$ level together with three
other sources, is similar to the first plausible non-stellar neutrino source,
TXS 0506+056, in being also a masquerading BL Lac object, i.e., intrinsically
a flat-spectrum radio quasar with hidden broad lines and a standard accretion
disk. We point out that these two sources share other properties, including
spectral energy distribution, high powers, parsec scale properties, and
possibly radio morphology. We speculate that the relatively rare combination
of proton-loaded jets, possibly typical of high-excitation sources, and
efficient particle acceleration processes, related to their relatively high
synchrotron peak frequencies, might favour neutrino production in these two
sources. GB6 J1542+6129, which has also recently appeared twice in a list of
IceCube associations, seems also to belong to this rare blazar sub-class,
which includes at most $\approx 20$ Fermi-4LAC blazars.

###### keywords:

neutrinos — radiation mechanisms: non-thermal — galaxies: active — BL Lacertae
objects: general — radio continuum: galaxies — gamma-rays: galaxies

††pubyear: 2022††pagerange: PKS 1424+240: yet another masquerading BL Lac
object as a possible IceCube neutrino source–PKS 1424+240: yet another
masquerading BL Lac object as a possible IceCube neutrino source

##  1 Introduction

Nine years ago the IceCube Neutrino Observatory111http://icecube.wisc.edu
detected the first high-energy astrophysical neutrinos of likely extragalactic
origin with energies up to $>1$ PeV ($10^{15}$ eV) and since then has produced
a steady list of events (e.g. Aartsen et al., 2020, and references therein).
At present, however, only two astronomical objects have been associated with a
significance larger than $\sim 3\,\sigma$ with these astrophysical neutrinos.
Namely, the blazar TXS 0506+056 (IceCube Collaboration et al., 2018a; IceCube
Collaboration, 2018b) at $z=0.3365$ and the local ($z=0.004$) Seyfert 2 galaxy
NGC 1068. Aartsen et al. (2020), in fact, have reported on an excess of
neutrinos at the $2.9\,\sigma$ level from the direction of NGC 1068 and a
$3.3\,\sigma$ excess in the northern sky due to significant p-values in the
directions of NGC 1068 and three blazars: TXS 0506+056, PKS 1424+240, and GB6
J1542+6129.

Padovani et al. (2019) showed that TXS 0506+056, despite appearances, is not a
blazar of the BL Lac type but instead a masquerading BL Lac object, namely a
flat-spectrum radio quasar (FSRQ222Based on optical spectroscopy blazars are
classified into FSRQs and BL Lac objects, with the former displaying strong,
broad, quasar-like emission lines and the spectra of the latter being often
completely featureless and sometimes exhibiting weak absorption and emission
lines (e.g. Urry & Padovani, 1995).) whose emission lines are swamped by a
very bright, Doppler-boosted jet, unlike “real” BL Lacs, which are instead
intrinsically weak-lined. This is extremely relevant for two reasons: (1)
“real” BL Lacs and FSRQs belong to two very different physical classes, i.e.,
objects without and with high-excitation emission lines in their optical
spectra, referred to as low-excitation (LEGs) and high-excitation galaxies
(HEGs), respectively (e.g. Padovani et al., 2017, and references therein); (2)
masquerading BL Lacs, being HEGs, benefit from several radiation fields
external to the jet (i.e., the accretion disc, photons reprocessed in the
broad-line region (BLR) or from the dusty torus), which, by providing more
targets for the protons might enhance neutrino production as compared to LEGs.

The aims of this Letter are to: (1) investigate if the second blazar in the
Aartsen et al. (2020)’s list, PKS 1424+240, also qualifies as a masquerading
BL Lac object; (2) check if TXS 0506+056 and PKS 1424+240 (and GB6 J1542+6129)
share other relevant properties. We use a $\Lambda$CDM cosmology with Hubble
constant $H_{0}=70$ km s-1 Mpc-1, matter density $\Omega_{\rm m,0}=0.3$, and
dark energy density $\Omega_{\Lambda,0}=0.7$. Spectral indices are defined by
$S_{\nu}\propto\nu^{-\alpha}$ where $S_{\nu}$ is the flux at frequency $\nu$.

##  2 The source

###  2.1 Main astronomical data

PKS 1424+240 is one of the most distant TeV-detected blazars, ranked number
seven in TeVCat333http://tevcat.uchicago.edu/ at the time of writing, and “a
rare example of a luminous HBL444Blazars are sub-divided on the basis of the
rest-frame frequency of their low-energy (synchrotron) hump ($\nu_{\rm
peak}^{S}$) into low- (LBL/LSP: $\nu_{\rm peak}^{S}$ $<10^{14}$ Hz [$<$ 0.41
eV]), intermediate- (IBL/ISP: $10^{14}$ Hz$~{}<$ $\nu_{\rm peak}^{S}$
$<10^{15}$ Hz [0.41 – 4.1 eV)], and high-energy (HBL/HSP: $\nu_{\rm peak}^{S}$
$>10^{15}$ Hz [$>$ 4.1 eV]) peaked sources respectively (Padovani & Giommi,
1995; Abdo et al., 2010).” (Cerruti et al., 2017). Its $\nu_{\rm peak}^{S}$
has been estimated to be around $10^{15}-10^{16}$ Hz, depending on its state
(Abdo et al., 2010; Archambault et al., 2014).

Figure 1: The optical spectrum of PKS 1424+240 obtained at the Gran Telescopio
Canarias with the OSIRIS spectrograph (Paiano et al., 2017). The continuum is
well described by a power law with spectral index $\alpha=0.9$ (equivalent to
F${}_{\lambda}\propto\lambda^{-1.1}$). Two weak emission lines of [O ii] 3727
Å and [O iii] 5007 Å (see insets) are detected at z = 0.6047. Main telluric
lines are indicated by $\oplus$.

The redshift of PKS 1424+240 has a confusing history. The issue has been
settled by Paiano et al. (2017), who detected two faint emission lines at 5981
and 8034 Å identified with [O ii] 3727 Å and [O iii] 5007 Å, with equivalent
widths 0.05 and 0.10 Å respectively, at $z=0.6047$ (Fig. 1). This redshift is
also consistent with the presence of a group of galaxies at $z\sim 0.6$ very
likely associated with the source (Rovero et al., 2016). These detections
imply line powers $L_{\rm[\text{O\,{ii}}]}\sim 4\times 10^{41}$ erg s-1 and
$L_{\rm[\text{O\,{iii}}]}\sim 10^{42}$ erg s-1.

As no sign of the host galaxy is apparent in the optical spectrum we cannot
decompose it into a non-thermal power law and an elliptical galaxy template to
estimate the black hole mass, $M_{\rm BH}$ (e.g. Padovani et al., 2022,
hereafter P22; see their Appendix for details). However, we can place an upper
limit on $L_{\rm\text{Mg\,{ii}}}<5\times 10^{42}$ erg s-1 assuming this line
has a full width at half maximum of 10,000 km s-1. By using eq. (5) and Table
2 of Shaw et al. (2012) we then derive $M_{\rm BH}<8\times 10^{8}~{}\rm
M_{\odot}$ (where $\rm M_{\odot}$ is one solar mass). Note that this is a very
robust upper limit, only slightly above the value one would get assuming the
host galaxy to be a typical giant elliptical ($M_{\rm BH}\sim 6.3\times
10^{8}M_{\odot}$: P22), which translates into an Eddington power $L_{\rm
Edd}<10^{47}$ erg s-1, with $L_{\rm Edd}=1.26\times 10^{46}~{}(M/10^{8}\rm
M_{\odot})$ erg s-1.

###  2.2 Source characterization

Figure 2: $P_{\rm 1.4GHz}$ vs. $L_{\rm[\text{O\,{ii}}]}$ for the P22’s sample
(black filled circles), with masquerading sources highlighted (larger empty
circles). Sources for which $L_{\rm[\text{O\,{ii}}]}$ has been estimated from
$L_{\rm[\text{O\,{iii}}]}$ are denoted by black empty squares. The magenta
triangle is PKS 1424+240, a brown circle indicates TXS 0506+056, while the
green filled circle denotes the single FSRQ in the sample. The solid blue line
is the locus of jetted (radio-loud) quasars, with the two dotted lines
indicating a spread of 0.5 dex, which includes most of the points in Fig. 4 of
Kalfountzou et al. (2012). Arrows denote upper limits on
$L_{\rm[\text{O\,{ii}}]}$. Adapted from Fig. 1 of P22.

We claim that PKS 1424+240 is yet another example of a masquerading BL Lac
object linked to IceCube neutrinos. Following Padovani et al. (2019) and P22
this classification is based on four criteria:

  1. 1.

its location on the radio power – emission line power, $P_{\rm 1.4GHz}$ –
$L_{\rm[\text{O\,{ii}}]}$, diagram (Fig. 4 of Kalfountzou et al., 2012), which
defines the locus of jetted (radio-loud) quasars. With a catalogued 1.4 GHz
flux density of 0.43 Jy and a radio spectral index $\alpha_{\rm r}\sim 0.4$
(NASA/IPAC Extragalactic Database) $P_{\rm 1.4GHz}\sim 5\times 10^{26}$ W
Hz-1, which, together with $L_{\rm[\text{O\,{ii}}]}\sim 4\times 10^{41}$ erg
s-1, puts this source exactly on the locus, as shown by Fig. 2;

  2. 2.

a radio power $P_{\rm 1.4GHz}>10^{26}$ W Hz-1, typical of HEGs;

  3. 3.

an Eddington ratio $L/L_{\rm Edd}>0.03$, hence $>0.01$, which is also one of
the defining characteristics of HEGs. This is the ratio between the accretion-
related observed luminosity and the Eddington luminosity. The estimation of
the former is derived from $L_{\rm[\text{O\,{ii}}]}$ and
$L_{\rm[\text{O\,{iii}}]}$, as detailed in Padovani et al. (2019). We obtain
$L\sim 3\times 10^{45}$ erg s-1 and a BLR luminosity $L_{\rm BLR}\sim 10^{44}$
erg s-1;

  4. 4.

a $\gamma$-ray Eddington ratio $L_{\gamma}/L_{\rm Edd}>2$ and therefore above
the proposed dividing line between “real” BL Lacs and FSRQs (0.1: Sbarrato et
al. 2012). In fact, $L_{\gamma}=1.7\times 10^{47}$ erg s-1 (0.1 – 100 GeV), as
derived from the Fermi-LAT Fourth Source Catalog (4FGL-DR2) energy flux and
photon index (Abdollahi et al., 2020; Ballet et al., 2020). We stress that,
due to its very large $L_{\gamma}$, PKS 1424+240 falls outside the boundaries
of Fig. 2 of P22 (the $L_{\gamma}-P_{\rm 1.4GHz}$ plot) and is quite extreme
in terms of its location in Fig. 3 (the $\nu_{\rm peak}^{S}$ – $L_{\gamma}$
plot) in the same paper.

In short, PKS 1424+240 meets all the above criteria to be classified as a
masquerading BL Lac object.

Figure 3: The archival SED of PKS 1424+240 (red points); for comparison the
SED of TXS 0506+056 is overlaid (light grey points). The overall SEDs of the
two sources are quite similar.

###  2.3 Similarities between PKS 1424+240 and TXS 0506+056

Apart from being both masquerading BL Lacs, there are quite a few other
commonalities between PKS 1424+240 and TXS 0506+056, which make them somewhat
special and unusual blazars. Namely:

  1. 1.

they are both very powerful IBLs/HBLs, as TXS 0506+056 has $P_{\rm 1.4GHz}\sim
1.5\times 10^{26}$ W Hz-1 and $L_{\gamma}\sim 2.9\times 10^{46}$ erg s-1,
$\sim 3-6$ times smaller than the PKS 1424+240 values;

  2. 2.

the overall SEDs of the two sources are quite similar, both in normalisation
and shape. This is clear from Fig. 3, which shows the SED of PKS 1424+240
(red) together with that of TXS 0506+056 (grey). The low-energy parts of the
SEDs are almost identical, while there are some difference in the X-ray and
$\gamma$-ray bands, mostly because of the higher average $\nu_{\rm peak}^{S}$
of PKS 1424+240;

  3. 3.

their parsec scale properties, inferred from Very Long Baseline Interferometry
(VLBI) studies, are much more typical of HBLs than of FSRQs. The observed
brightness temperature of the centimetre-VLBI core, $T_{\rm B}$, is in both
cases relatively low, with median values $\sim 6\times 10^{10}$ $\rm K$ in PKS
1424+240 and $\sim 7\times 10^{10}$ $\rm K$ in TXS 0506+056 (Homan et al.,
2021), i.e. close to the equipartition value (${\sim}5\times 10^{10}$ $\rm K$:
Readhead 1994) and therefore indicative of a modest Doppler boosting. This is
also supported by proper motion measurements, revealing low apparent speeds of
$v_{\rm app}/c\equiv\beta_{\rm app}=2.83\pm 0.89$ (Lister et al., 2019) and
$1.07\pm 0.14$ (Lister et al., 2021) for PKS 1424+240 and TXS 0506+056
respectively. FSRQs, especially the $\gamma$-ray powerful ones, are instead
characterized by median core $T_{\rm B}$ approaching and often exceeding the
inverse Compton catastrophe limit555The inverse Compton catastrophe is the
rapid cooling through inverse Compton scattering of a synchrotron-emitting
region, which implies a threshold temperature $\sim 10^{12}$ $\rm K$
(Kellermann, Pauliny-Toth, & Williams, 1969). due to strong Doppler boosting
(Kovalev et al., 2009), as well as by much higher maximum apparent speeds
($\beta_{\rm app}\sim 10-30$, e.g. Jorstad et al. 2017). Homan et al. (2021)
estimate relatively low Doppler and Lorentz factors, $\delta\sim 1.8$ and
$\Gamma\sim 1.5$ for TXS 0506+056, and $\delta\sim 1.4$ and $\Gamma\sim 4$ and
for PKS 1424+240, as typically found for HBLs (Piner & Edwards, 2018, and
references therein).666Homan et al. (2021) derive $\delta$ by assuming that
all sources have the same intrinsic median brightness temperature equal to the
sample median, $T_{\rm B,int}=10^{10.609\pm 0.067}$ K; hence $\delta=T_{\rm
B,obs}/T_{\rm B,int}$. This population-wide approach, while good for average
values, might not be the most appropriate one for individual sources. A
relatively low Lorentz factor is also consistent with the large apparent jet
opening angles in both sources (28∘ and 65∘ for TXS 0506+056 and PKS 1424+240,
respectively, above the median value of 22.6∘ found for the $\gamma$-ray
sources in the MOJAVE sample, see Pushkarev et al., 2017). Indeed, the
intrinsic aperture of the boosting cone is found to be inversely proportional
to $\Gamma$, meaning that slower jets will appear less collimated. In
conclusion, although these two sources have a standard (Shakura - Sunyaev)
disk, as FSRQs, the VLBI properties of their jets are not FSRQ-like;

  4. 4.

their extended radio power, $P_{\rm ext}$, is relatively high, placing the
sources in an intermediate regime between Fanaroff-Riley (FR) I and II objects
(Fanaroff & Riley, 1974), or rather closer to the FRII range ($\log P_{\rm
ext}\gtrsim 25.5$ W Hz-1 at 1.4 GHz). Very Large Array (VLA) imaging of PKS
1424+240 at 1.4 GHz reveals a symmetric two-sided structure, with an extended
flux density ${\sim}130\,\rm mJy$ (Fig. 12, Rector, Gabuzda, & Stocke, 2003),
which implies $\log P_{\rm ext}=26.3$ W Hz-1 ($\alpha=0.8$), at the lower end
of FRII sources. No extended radio emission has been resolved in TXS 0506+056
on kiloparsec (kpc) scales, but an estimate of the extended radio flux density
can be obtained by comparing VLBI data with almost contemporaneous single-dish
measurements obtained with the Owens Valley Radio Observatory telescope. Li et
al. (2020) have shown that the extended emission contributes between 1 and 18
per cent of the total, from which, since the smallest percentage is tied to
the largest total flux density (and viceversa), we estimate an extended flux
density in the range 15 – 54 mJy i.e. $\approx 35$ $\rm mJy$ at 15 GHz. This
implies $\log P_{\rm ext}=25.9^{+0.2}_{-0.4}$ W Hz-1 at 1.4 GHz for
$\alpha=0.8$, which is close to the FRI/FRII break. Sources classified as BL
Lacs but showing intermediate or FRII-like extended radio powers and
morphologies are common among bright blazar samples (Rector & Stocke, 2001;
Kharb, Lister, & Cooper, 2010). However, this is mainly true for LBLs, while
the properties of HBLs match those of their expected parent population, the
FRI radio galaxies (Rector & Stocke, 2001; Giroletti et al., 2004). For
example, the HBLs in the Einstein Extended Medium Sensitivity Survey,
characterised by $0.045\leq z\leq 0.638$ ($\langle z\rangle\sim 0.3$), have
$\log P_{\rm ext}\sim 23.7$ (and $\leq 25$) W Hz-1 at 1.4 GHz (Rector &
Stocke, 2001). Therefore, PKS 1424+240 and TXS 0506+056 appear to be rare
examples of IBLs/HBLs with FRII-like radio powers. If the extended radio power
is a good proxy of the total jet power (e.g. Willott et al., 1999), the latter
should also be in the FRII regime. High jet powers $\sim 10^{45}-10^{46}$ erg
s-1 are also derived for both sources based on SED modeling (e.g. Aleksić et
al. 2014; Ansoldi et al. 2018), much larger than the values, for example, of
the low-redshift HBLs MKN 421 and MKN 501 ($\sim 6\times 10^{43}$ erg s-1:
Potter & Cotter 2015);

  5. 5.

the issue of the radio morphology is, instead, difficult to address. TXS
0506+056 is unresolved on kpc scales but the properties of the VLBI jet may be
consistent with the development of an FRI morphology. In fact, the rapid
increase of the jet opening angle with distance observed in VLBI images at 43
GHz has been interpreted by Ros et al. (2020) as possible evidence for jet
deceleration taking place at $\sim 70-140$ pc from the BH, a process known to
lead to the formation of an FRI jet (Laing & Bridle, 2014). The VLA image of
PKS 1424+240 shown by Rector, Gabuzda, & Stocke (2003) presents mixed
features, which would need to be better examined at higher resolution: the
structure is reminiscent of an edge-brightened FRII, but the rough symmetry of
the two-sided jet may be indicative of deceleration as observed in FRIs. In
conclusion, one interesting possibility that needs to be further explored is
the fact that PKS 1424+240 and TXS 0506+056 might belong to the rare class of
efficiently accreting sources, which develop an FRI jet morphology, i.e., FRI-
HEGs (e.g. Heckman & Best, 2014). A sample of 1.4 GHz selected extended radio
sources from the Sloan Digital Sky Survey with $0.03<z<0.1$ contains, for
example, only 5/245, i.e. 2 per cent, of FRI-HEGs (Miraghaei & Best, 2017).
Even if that were the case, however, these would have to be sources with FRI
morphology and FRII-like jet powers, as discussed above.

###  2.4 GB6 J1542+6129

IceCube Collaboration et al. (2021) have performed a search for flare neutrino
emission in the 10 years of IceCube data and reported a cumulative time-
dependent neutrino excess in the northern hemisphere at the 3 $\sigma$ level
associated with only four sources: a radio galaxy, M87, two blazars, TXS
0506+056 and GB6 J1542+6129, and the Seyfert 2 galaxy NGC 1068, these last
three sources being also associated with neutrinos by Aartsen et al. (2020).
GB6 J1542+6129, an IBL (Giommi & Padovani, 2021), has then appeared twice in
two years in a list of $3-3.3\,\sigma$ IceCube associations and is therefore
also worth investigating. The redshift of this object is not available: all we
know is that $0.34\leq z\leq 1.76$ (Shaw et al., 2013), which translates into
$P_{\rm 1.4GHz}=10^{26.4\pm 0.7}$ W Hz-1 and $L_{\gamma}/L_{\rm
Edd}=10^{0.3\pm 0.9}$ (for $M_{\rm BH}=6.3\times 10^{8}M_{\odot}$, i.e.
assuming the host galaxy to be a typical giant elliptical: Section 2.1). Both
values are perfectly consistent with a masquerading BL Lac classification
(Section 2.2). From the data given by Homan et al. (2021) we derive $T_{\rm
B}=1.7-3.6\times 10^{11}$ K, $\beta_{\rm app}=0.8-3.1$, $\delta=4.2-8.7$, and
$\Gamma=2.3-5.0$ for the adopted redshift range. $\delta$ and $\Gamma$ are
therefore also relatively small, which suggests that GB6 J1542+6129 shares
similar radio properties with the other two sources.

##  3 Discussion and Summary

The properties just described highlight the peculiarities of PKS 1424+240 and
TXS 0506+056, which may be relevant in the context of neutrino emission
models. In particular, a fundamental question concerns the jet particle
composition, since the emission of neutrinos requires the presence of
relativistic protons. The loading of heavy particles in the jet may occur
either in the jet formation process or during jet propagation, as a
consequence of entrainment of gas and stars from the environment. Jets formed
through the Blandford-Znajeck mechanism (Blandford & Znajek, 1977) are thought
to be light, being mainly composed by electron-positron pairs, while disk-
launched jets (e.g. Blandford & Payne, 1982) are directly loaded with disk
material, and are likely to be heavier and slower (see also Hawley & Krolik,
2006; McKinney, 2006; Broderick & Tchekhovskoy, 2015). If both mechanisms take
place at the same time in the AGN, mass exchange, possibly favoured by the
development of magneto-hydrodynamic instabilities, may occur between the
central relativistic jet and the surrounding heavy disk wind (see e.g. the
discussion by Sikora, Stawarz, & Lasota, 2007). Mixing and entrainment of
heavy particles from the ambient medium are also possible, particularly on
larger scales when a jet decelerates and disrupts. Indeed, Croston, Ineson, &
Hardcastle (2018) were able to constrain the composition of the radio-lobe
plasma in samples of FRI and FRII sources, showing that the decelerating jets
in FRIs require a stronger proton content with respect to FRIIs. These authors
also showed that this phenomenon is independent of the AGN accretion mode,
being strictly related to the strong interplay between the jet and the
environment, which characterizes FRI sources.

However, the most energetic phenomena in the jet, leading to the production of
$\gamma$-rays up to the TeV regime, are known to occur at smaller distances
from the BH, on sub-parsec and parsec scales (e.g. Madejski & Sikora, 2016,
and references therein). Therefore, it is the jet composition closer to the
central engine that may be most relevant in relation to neutrino production.
Contrary to the scenario described above for the larger scales, modeling of
the jet broadband emission (Celotti & Fabian, 1993; Sikora et al., 2005) and
observational studies of the jet circular polarization (Homan et al., 2009)
suggest that the dynamics of powerful jets is dominated by protons on VLBI
scales. The parsec scale properties of the FRI jet in the LEG M 87, on the
other hand, are consistent with the dominance of a pair plasma (Reynolds et
al., 1996; Broderick & Tchekhovskoy, 2015) and, in general, the assumption of
light jets appears adequate for the modeling of jet deceleration in FRIs
(Laing & Bridle, 2002, 2014). This possible difference in the jet composition
may reflect a difference in the jet formation mechanism in LEGs (FRIs and
FRIIs) and HEGs (mostly FRIIs and a few FRIs). Based on the analysis of the
jet collimation profiles, Boccardi et al. (2021) have suggested that a more
extended and prominent disk wind surrounds the relativistic spine in HEGs,
while the jet expansion profiles in most LEGs are consistent with a jet origin
in the vicinity of the ergosphere, as directly observed in M 87. These results
support the idea that jets in HEGs are more heavily loaded with protons on
sub-parsec and parsec scales.

It then follows that the production of neutrinos might be favoured in sources
like PKS 1424+240 and TXS 0506+056. If jets produced by HEGs are indeed more
heavily loaded with protons, and if protons and electrons are accelerated in
the same regions and through the same mechanisms, then jets from IBL/HBL-HEGs
may be those where the two conditions necessary for neutrino production are
best met. Namely: 1) a high proton loading; 2) the occurrence of extremely
efficient particle acceleration processes. In this scenario FSRQs would not
fulfil the second requirement. In addition, as we mentioned above, the
contribution from several external radiation fields typical of HEGs might
further enhance neutrino production in these masquerading BL Lac objects.

Any physical mechanism, which may be invoked to explain the emission up to TeV
energies in PKS 1424+240 and TXS 0506+056 will have to be reconciled, however,
with the low $\Gamma$ and $\delta$ factors observed on parsec scales. Most
solutions to the so-called Doppler factor crisis have been proposed so far for
the case of “classic” low-power HBLs. These include the near-core decelerating
jet models (Piner & Edwards, 2004, and references therein) and the shock-in-
jet models, in which the low apparent speeds are proposed to reflect the
formation of quasi-stationary shocks in the jet (Hervet, Boisson, & Sol, 2016,
and references therein). It is a subject for future investigation whether such
solutions are also applicable to the case of the high-power jets in PKS
1424+240 and TXS 0506+056.

In summary, PKS 1424+240, as TXS 0506+056, is another example of a
masquerading BL Lac object associated with an IceCube neutrino excess. The two
sources share also other properties, including spectral energy distribution,
high powers, parsec scale properties, and possibly radio morphology. We
suggest that the combination of proton-loaded jets, which might be typical of
high-excitation sources, and efficient particle acceleration, implied by their
relatively high synchrotron peak frequencies, might favour neutrino
production. Finally, we note that GB6 J1542+6129, which has been also recently
associated to neutrinos, seems also to belong to this very rare blazar sub-
class, although a redshift determination would help to confirm this.

Based on their very high radio and $\gamma$-ray powers we estimate that
similar sources make up only 2.1 per cent of the total IBL and HBL
population777This is the fraction of sources with radio and $\gamma$-ray
powers larger than those of TXS 0506+056 in the IBL/HBL sample put together by
Giommi & Padovani (2021). GB6 J1542+6129 has $P_{\rm 1.4GHz}=10^{26.4\pm 0.7}$
W Hz-1 and $L_{\gamma}=10^{47.2\pm 0.9}$ erg s-1.. Since IBLs plus HBLs
comprise $\sim 42$ per cent of the Fermi-4LAC (clean sample) blazars with SED
classification (Ajello et al., 2020; Lott, Gasparrini, & Ciprini, 2020),
blazars of the type IceCube has associated with neutrinos constitute at most 1
per cent of the $\gamma$-ray selected population ($\lesssim 20$ sources),
since we have not included in our calculation their other peculiar properties.

## Acknowledgments

We thank Maria Petropoulou for her comments on the paper. This work is
supported by the Deutsche Forschungsgemeinschaft through grant SFB 1258
“Neutrinos and Dark Matter in Astro- and Particle Physics”.

## Data Availability

The flux-calibrated and de-reddened spectrum is available in the online
database ZBLLAC (http://web.oapd.inaf.it/zbllac/).

## References

  * Aartsen et al. (2020) Aartsen M. G., Ackermann M., Adams J., Aguilar J. A., Ahlers M., Ahrens M., Alispach C., et al., 2020, PhRvL, 124, 051103
  * Abdo et al. (2010) Abdo A. A., Ackermann M., Ajello M., et al., 2010, ApJ, 716, 30
  * Abdollahi et al. (2020) Abdollahi S., Acero F., Ackermann M., Ajello M., Atwood W. B., Axelsson M., Baldini L., et al., 2020, ApJS, 247, 33
  * Ajello et al. (2020) Ajello M., Angioni R., Axelsson M., Ballet J., Barbiellini G., Bastieri D., Becerra Gonzalez J., et al., 2020, ApJ, 892, 105
  * Aleksić et al. (2014) Aleksić J., Ansoldi S., Antonelli L. A., Antoranz P., Babic A., Bangale P., Barres de Almeida U., et al., 2014, A&A, 567, A135
  * Ansoldi et al. (2018) Ansoldi S., Antonelli L. A., Arcaro C., Baack D., Babić A., Banerjee B., Bangale P., et al., 2018, ApJL, 863, L10
  * Archambault et al. (2014) Archambault S., Aune T., Behera B., Beilicke M., Benbow W., Berger K., Bird R., et al., 2014, ApJL, 785, L16
  * Ballet et al. (2020) Ballet J., Burnett T. H., Digel S. W., Lott B., The Fermi-LAT Collaboration, 2020, arXiv:2005.11208
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  * Boccardi et al. (2021) Boccardi B., Perucho M., Casadio C., Grandi P., Macconi D., Torresi E., Pellegrini S., et al., 2021, A&A, 647, A67
  * Broderick & Tchekhovskoy (2015) Broderick A. E., Tchekhovskoy A., 2015, ApJ, 809, 97
  * Celotti & Fabian (1993) Celotti A., Fabian A. C., 1993, MNRAS, 264, 228
  * Cerruti et al. (2017) Cerruti M., Benbow W., Chen X., Dumm J. P., Fortson L. F., Shahinyan K., 2017, A&A, 606, A68
  * Croston, Ineson, & Hardcastle (2018) Croston J. H., Ineson J., Hardcastle M. J., 2018, MNRAS, 476, 1614
  * Fanaroff & Riley (1974) Fanaroff B. L., Riley J. M., 1974, MNRAS, 167, 31P
  * Ghisellini, Tavecchio, & Chiaberge (2005) Ghisellini G., Tavecchio F., Chiaberge M., 2005, A&A, 432, 401
  * Giommi & Padovani (2021) Giommi P., Padovani P., 2021, Universe, 7(12), 492
  * Giroletti et al. (2004) Giroletti M., Giovannini G., Taylor G. B., Falomo R., 2004, ApJ, 613, 752
  * Hardcastle & Croston (2020) Hardcastle M. J., Croston J. H., 2020, NewAR, 88, 101539
  * Hawley & Krolik (2006) Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103
  * Heckman & Best (2014) Heckman T. M., Best P. N., 2014, ARA&A, 52, 589
  * Hervet, Boisson, & Sol (2016) Hervet O., Boisson C., Sol H., 2016, A&A, 592, A22
  * Homan et al. (2009) Homan D. C., Lister M. L., Aller H. D., Aller M. F., Wardle J. F. C., 2009, ApJ, 696, 328
  * Homan et al. (2021) Homan D. C., Cohen M. H., Hovatta T., Kellermann K. I., Kovalev Y. Y., Lister M. L., Popkov A. V., et al., 2021, ApJ, in press (arXiv:2109.04977)
  * IceCube Collaboration et al. (2018a) IceCube Collaboration, 2018, Science, 361, 147
  * IceCube Collaboration (2018b) IceCube Collaboration et al., 2018, Science, 361, eaat1378
  * IceCube Collaboration et al. (2021) IceCube Collaboration, Abbasi R., Ackermann M., Adams J., Aguilar J. A., Ahlers M., Ahrens M., Alispach C., et al., 2021, ApJL, 920, L45
  * Jorstad et al. (2017) Jorstad S. G., Marscher A. P., Morozova D. A., Troitsky I. S., Agudo I., Casadio C., Foord A., et al., 2017, ApJ, 846, 98
  * Kalfountzou et al. (2012) Kalfountzou E., Jarvis M. J., Bonfield D. G., Hardcastle M. J., 2012, MNRAS, 427, 2401
  * Kellermann, Pauliny-Toth, & Williams (1969) Kellermann K. I., Pauliny-Toth I. I. K., Williams P. J. S., 1969, ApJ, 157, 1
  * Kharb, Lister, & Cooper (2010) Kharb P., Lister M. L., Cooper N. J., 2010, ApJ, 710, 764
  * Kovalev et al. (2009) Kovalev Y. Y., Aller H. D., Aller M. F., Homan D. C., Kadler M., Kellermann K. I., Kovalev Y. A., et al., 2009, ApJL, 696, L17
  * Laing & Bridle (2002) Laing R. A., Bridle A. H., 2002, MNRAS, 336, 1161
  * Laing & Bridle (2014) Laing R. A., Bridle A. H., 2014, MNRAS, 437, 3405
  * Li et al. (2020) Li X., An T., Mohan P., Giroletti M., 2020, ApJ, 896, 63
  * Lister et al. (2019) Lister M. L., Homan D. C., Hovatta T., Kellermann K. I., Kiehlmann S., Kovalev Y. Y., Max-Moerbeck W., et al., 2019, ApJ, 874, 43
  * Lister et al. (2021) Lister M. L., Homan D. C., Kellermann K. I., Kovalev Y. Y., Pushkarev, A. B., Ros, E., Savolainen, T., 2021, ApJ, in press (arXiv:2108.13358)
  * Lott, Gasparrini, & Ciprini (2020) Lott B., Gasparrini D., Ciprini S., 2020, arXiv:2010.08406
  * Madejski & Sikora (2016) Madejski G., Sikora M., 2016, ARA&A, 54, 725
  * McKinney (2006) McKinney J. C., 2006, MNRAS, 368, 1561
  * Miraghaei & Best (2017) Miraghaei H., Best P. N., 2017, MNRAS, 466, 4346
  * Padovani & Giommi (1995) Padovani P., Giommi P., 1995, ApJ, 444, 567
  * Padovani et al. (2017) Padovani P., et al., 2017, A&AR, 25, 2
  * Padovani et al. (2019) Padovani P., Oikonomou F., Petropoulou M., Giommi P., Resconi E., 2019, MNRAS, 484, L104
  * Padovani et al. (2022) Padovani P., Giommi P., Falomo R., Oikonomou F., Petropoulou M., Glauch T., Resconi E., Treves A., Paiano S., 2022, MNRAS, 510, 2671 (P22)
  * Paiano et al. (2017) Paiano S., Landoni M., Falomo R., Treves A., Scarpa R., Righi C., 2017, ApJ, 837, 144
  * Piner & Edwards (2004) Piner B. G., Edwards P. G., 2004, ApJ, 600, 115
  * Piner & Edwards (2018) Piner B. G., Edwards P. G., 2018, ApJ, 853, 68
  * Potter & Cotter (2015) Potter W. J., Cotter G., 2015, MNRAS, 453, 4070
  * Pushkarev et al. (2017) Pushkarev A. B., Kovalev Y. Y., Lister M. L., Savolainen T., 2017, MNRAS, 468, 4992
  * Readhead (1994) Readhead A. C. S., 1994, ApJ, 426, 51
  * Rector & Stocke (2001) Rector T. A., Stocke J. T., 2001, AJ, 122, 565
  * Rector, Gabuzda, & Stocke (2003) Rector T. A., Gabuzda D. C., Stocke J. T., 2003, AJ, 125, 1060
  * Reynolds et al. (1996) Reynolds C. S., Fabian A. C., Celotti A., Rees M. J., 1996, MNRAS, 283, 873
  * Ros et al. (2020) Ros E., Kadler M., Perucho M., Boccardi B., Cao H.-M., Giroletti M., Krauß F., et al., 2020, A&A, 633, L1
  * Rovero et al. (2016) Rovero A. C., Muriel H., Donzelli C., Pichel A., 2016, A&A, 589, A92
  * Sbarrato et al. (2012) Sbarrato T., Ghisellini G., Maraschi L., Colpi M., 2012, MNRAS, 421, 1764
  * Shaw et al. (2012) Shaw M. S., Romani R. W., Cotter G., Healey S. E., Michelson P. F., Readhead A. C. S., Richards J. L., et al., 2012, ApJ, 748, 49
  * Shaw et al. (2013) Shaw M. S., Romani R. W., Cotter G., Healey S. E., Michelson P. F., Readhead A. C. S., Richards J. L., et al., 2013, ApJ, 764, 135
  * Sikora et al. (2005) Sikora M., Begelman M. C., Madejski G. M., Lasota J.-P., 2005, ApJ, 625, 72
  * Sikora, Stawarz, & Lasota (2007) Sikora M., Stawarz Ł., Lasota J.-P., 2007, ApJ, 658, 815
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, PASP, 107, 803
  * Willott et al. (1999) Willott C. J., Rawlings S., Blundell K. M., Lacy M., 1999, MNRAS, 309, 1017