COSMIC: An Ethernet-based Commensal, Multimode Digital Backend on the Karl G. Jansky Very Large Array for the Search for Extraterrestrial Intelligence

Each antenna in the VLA outputs data via the VLA DTS, as described in depth in Durand (2001) and Freund (2002). For each antenna, digitized voltage data for the two polarizations and multiple frequency tunings (i.e., intermediate frequencies, IFs, or LO tunings) are framed, augmented with timing and error-check metadata, and transmitted over a single fiber. Each DTS fiber carries 120 Gb s⁻¹ of data—consisting of 96 Gb s⁻¹ of digitized sample data plus protocol overload—over 12 wavelength-multiplexed 10 Gb s⁻¹ lanes, which utilize ∼1550 nm lasers spaced at 200 GHz. The DTS protocol is nonstandard but has the following features.

• 1.  
  Synchronization patterns are built into the protocol to allow a downstream receiver to perform clock recovery and read the data stream without additional reference clocks.
• 2.  
  The underlying data transport bit clock is synchronous with the VLA samplers, enabling a receiver to recover the original sampling clock.
• 3.  
  A timing pulse every 50 ms is embedded in the DTS streams, allowing streams from multiple antennas to be deterministically aligned in time.
• 4.  
  Checksums (the number of bits in a transmitted message) are transmitted within the data streams to enable a receiver to monitor transmission errors.

Since the DTS protocol has a custom design, it is readable by a low-level programmable logic chip, such as an FPGA. The role of the COSMIC optical interfacing system is to split the DTS fibers and manipulate the COSMIC copy of the DTS stream to interface with an off-the-shelf FPGA processing platform. At each antenna in the VLA, two IF bands exist, with a total bandwidth of up to 2 GHz in the 8 bit mode and 8 GHz in the 3 bit mode. Currently, COSMIC operates only when the telescope is operating in the 8 bit mode, utilizing up to 1024 MHz of bandwidth for each IF.

2.1.1. DTS Splitting

Each of the VLA antenna's DTS fibers is split with inexpensive optical components. However, it is of paramount importance that the fiber splitting (and resulting ∼3 dB attenuation of the DTS signal at each of the splitter outputs) does not adversely affect the ability of WIDAR or COSMIC to receive the data streams. To maintain the DTS power levels after splitting and, further, to allow flexible tuning of these levels to support the VLA configuration-dependent input power levels, the COSMIC deployment includes a tunable variable attenuator and fixed amplification of the DTS signals upstream of the fiber split.

After a period of testing and qualification, the Fiberstore M6200-25PA erbium-doped fiber amplifier (EDFA) was chosen for COSMIC and installed by the NRAO (Figure 4). This amplifier is a cost-effective, off-the-shelf product providing 25 dB of amplification and designed to support long-range wavelength-multiplexed data links. The M6200-25PA ¹⁷ EDFA also supports the addition of a remote-controllable variable attenuator. One amplifier is required for each VLA antenna, and these are housed in rack-mounted enclosures, with each enclosure supporting up to seven individual amplifiers (Figure 4).

Zoom In Zoom Out Reset image size

2.1.2. Electrical Conversion

The role of the COSMIC station processors is to receive the DTS data streams, provide station-level data processing, and output data as a stream of User Datagram Protocol/Internet Protocol (UDP/IP) packets to the downstream array processing system. One of the driving goals of the COSMIC implementation is to use as much off-the-shelf equipment as possible to minimize nonrecurring engineering costs and risks. For the COSMIC station processors, this meant finding a cost-effective commercial platform offering the following features.

• 1.  
  A powerful FPGA capable of interfacing with the custom DTS data format and providing sufficient computing resources for digital signal processing.
• 2.  
  Many QSFP+ connections to facilitate interfacing with the 336 input data fibers (12 per antenna), as well as further connections for outputting processed data.
• 3.  
  A form factor based on rack-mountable enclosures to minimize the need for thermal and mechanical engineering.

The AlphaData ADM-PCIe-9H7 ¹⁸ meets all of these requirements (Figures 5 and 7). It is based on an AMD Virtex Ultrascale+ xcvu37p FPGA (AMD 2023) and provides substantial signal processing resources (including 9024 multiplier cores, 340 Mb of on-die memory, and 8 GB of on-chip memory) and QSFP+ connectors (four onboard, plus eight available via AD-PCIE-FQSFP expansion cards/daughter boards). The ADM-PCIe-9H7 comes in an industry-standard PCIe form factor and is installable in most standard rack-mountable computer servers that support GPUs.

Zoom In Zoom Out Reset image size

The COSMIC system uses one ADM-PCIe-9H7 card, paired with two ADM-PCIe-9H7 quad-QSFP+ daughter cards, to process data from a pair of VLA antennas. In this configuration, a set of cards receives 24 lanes of DTS data via six QSFP+ interfaces and transmits up to 400 Gb s⁻¹ of data to the downstream processing system via a further four QSFP+ outputs, as show in Figure 6.

Zoom In Zoom Out Reset image size

The COSMIC system comprises of 15 cards to facilitate the processing of 28 antenna inputs, including the provision of two "hot spares." The FPGA cards are hosted on Tyan B7119F77V10E4HR-2T55-N servers, ¹⁹ which are designed to support up to 21 single-width PCIe cards in a 4U form factor. A single server is able to support five FPGA processor cards and their QSFP+ expansions (see Figure 7).

Zoom In Zoom Out Reset image size

Control of the FPGA cards is exposed to the rest of the COSMIC system via a REST (a set of architectural constraints used by developers) interface running on the FPGA host servers, with communications to the FPGAs running over the PCIe bus by means of Linux drivers provided by the FPGA vendor, Xilinx/AMD.

2.2.1. Station Processing

The job of the COSMIC data interconnect system is to facilitate the transfer of antenna data from the station boards to the compute cluster. This interconnect system enables the "corner-turn" operation frequently employed in radio astronomy, in which data streams that are arranged in parallel by antenna are converted to instead be parallel by frequency. After this conversion, the downstream processing nodes can receive data from all antennas in the array. For COSMIC, however, only a subset of the total observation band is directed to each processing node. These data may be correlated, beamformed, or otherwise combined to leverage the distributed nature of the VLA's observing aperture.

As is typical of many modern radio telescopes, COSMIC uses off-the-shelf Ethernet switches as data interconnects (Hickish et al. 2016). Such hardware allows the corner-turn operation to be accomplished by simply addressing data frames entering the network to an appropriate destination processor such that data frames containing common frequency channels end up at a common processing node. With each COSMIC station processor outputting data over a pair of 100 GbE links per antenna, the full COSMIC data interconnect system comprises 58,100 GbE input streams and a potentially similar number of output streams.

For the greatest cost-effectiveness and availability of the switches, it is preferable to limit the COSMIC data network so that no more than 64 nodes in total need to connect to any one switch. COSMIC achieves this by dividing the data over two separate 100 GbE networks by ensuring that the station processors direct half of their total output bandwidth down each of their two available 100 GbE outputs (one switch for each LO tuning).

With this division, the COSMIC data interconnect system can be implemented as two completely independent networks, each with 28 inputs and a similar number of outputs. Each of these networks is built around a single 64-port 100 GbE switch. Such switches are widely available; specifically, COSMIC utilizes the N9K-C9364C switch, manufactured by Cisco.

The specifications of the compute cluster that processes the signals after COSMIC digital signal processing are shown in Table 1 and the bottom of Figure 7. The high data rate of the VLA and the real-time calibration and technosignature search goals set the main requirements for the compute cluster design. The data are transmitted through a pair of 100 GbE optical transceivers to each compute node in the cluster, where there are two nodes per GPU server. To facilitate real-time processing, each GPU node is allocated 32 × 1 MHz channels out of the total 1.024 GHz bandwidth for further channelizing, beamforming, and search processes.

In the current incarnation of COSMIC, the GPU compute cluster consists of 22 GPU servers, where each node is fitted with two network interface cards (NICs), 8 TB of nonvolatile memory (NVMe) storage, CPUs, and GPUs. Including a small amount of overhead, each NIC on each node ingests data at a rate of ∼3.4 Gb s⁻¹ during on-the-fly mapping scans, with 32 MHz of bandwidth distributed to each of the 22 servers. In the technosignature beamforming and search mode (Section 3.7), data are transferred from the FPGAs to the NVMe storage access buffers via a 100 GbE switch, as the NVMe drives provide fast read and write speeds. This design is based on benchmarks from commissioning, which indicate that a network bandwidth of at least 75% of line-rate 100 GbE can be ingested, making this a compute-limited SETI experiment.

Sets of 32 × 1 MHz coarse channels are distributed across multiple downstream compute nodes in full polarization from each of the VLA antennas to perform further data processing. The 100 GbE switch is used to transfer the data through the FPGA optical fiber interface to the GPU compute nodes. The data flow within the pipeline is managed using the hashpipe software package (an application that helps create threads and shared memory buffers between them; MacMahon et al. 2018).

When the VLA is either pointing toward a target, performing on-the-fly mapping, or calibrating, hashpipe will capture UDP packets from the FPGA, arrange them into the Green Bank Ultimate Pulsar Processing Instrument (GUPPI) raw format (Ford & Ray 2010), and write the data to the NVMe buffer on the GPU node. The GUPPI raw data format consists of a plain text header, loosely based on the FITS format (Pence et al. 2010), followed by a block of binary data. Each data block contains 2¹⁷ time samples, 32 coarse frequency channels, two polarizations, and up to 27 antennas' worth of data, amounting to an effective block size of 108 MB. The number of time samples is chosen to ensure that each batch processing action can be performed in its most efficient configuration. However, this is flexible and could be adjusted in the future.

The metadata received from the VLA contain information on the "intent" of the observation, specifying whether it is a target field or a calibration source. Moreover, the VLA offers a collection of well-defined observation intents, ²³ which are used as criteria within the YAML observation configuration file. This observation configuration file specifies a pair of destination hashpipe instances that perform the calibration process when the observation intent is indicated as "calibration." The target observation "intent" launches a different pair of hashpipe instances that record the data stream into GUPPI raw files. The observational YAML file also specifies different postprocessing procedures; data collected for the calibration process are collated, and gain solutions are updated, whereas data from target observations are beamformed and searched for technosignatures. The overall data flow is shown in Figure 8. Depending on the intent, the hashpipe automation pipeline directs the data flow to follow either the correlation and calibration pathway or the beamforming and search process.

The software across the multiple GPU servers is maintained by means of a centralized read-only operating system from which all GPU nodes boot using NetBoot, ²⁴ as described in MacMahon et al. (2018). This ensures that all of the GPU servers are running identical versions of the software and that all user access privileges are consistent across the processing systems.

COSMIC as a system is intended to run autonomously. YAML files specifying observation criteria are submitted by users to trigger the COSMIC system into action when metadata from Redis match those criteria. These YAML files contain information regarding observation frequencies and bandwidths (as these can be adjusted and do not need to incorporate the full output bandwidth by the VLA, depending on the science goals), observation intent, and scan duration. The YAML files also detail what sort of observation pipelines need to be engaged postrecording.

Accordingly, the only human intervention required on the recording end is the submission of observational YAML files to the system for monitoring.

The delay tracking and calibration processes (described in Section 3.3.1) operate autonomously when the delay corrections are <500 ns as expected, as in normal day-to-day operation. Delay phase tracking is chosen as the VLA updates the sky position via the multicast system. This will engage the source finder to automatically search the nearby sky, generate R.A. and decl. values for a point of interest, and feed this into the delay engine, which computes the required delays and phases for the F-engines ²⁵ to apply to the relevant data with a time to load. In this way, we phase up the data streams in the F-engines to that point of interest. This calibration process engages when a calibration target appears, triggering a correlation observation from which gains are derived, and the calibration delays and phases are produced and loaded.

If the delay corrections are greater than 500 ns, such as when the antennas move to a new configuration, human intervention is needed. In this scenario, an observation for calibration and cross-correlation will not yield the frequency resolution required to produce new delay values. Therefore, it is necessary to submit a special YAML file that instructs the hashpipe instance to record an observation of a bright calibrator source and save the GUPPI raw files without searching and then to use a bespoke Python correlator to upchannelize and produce calibration delay values. We then manually inspect the generated delays and submit them as the new default for the calibration process, if they are deemed appropriate. If the values are not appropriate, the data will be further investigated, and the system will not be used for science until new appropriate values are obtained.

When the VLA is pointed at and the observer has marked the observation intent as a flux density, bandpass, phase, or gain calibrator source, the raw voltages from each antenna are cross-correlated using the xGPU ²⁶ (Clark et al. 2013) software correlator to produce visibility data products in four polarizations (RR, LL, RL, and LR, where R and L refer to right and left circular polarization, respectively). The software correlator writes these data products into a data file in the UVH5 format, ²⁷ which contains all of the information required for real-time calibration but can be used in imaging applications as well.

For calibration, we follow standard calibration procedures utilized for interferometric telescopes, but these procedures are executed through a real-time and autonomous process. The NRAO provides detailed information about VLA calibration on their website for users of the telescope, ²⁸ which we follow for the general purpose of calibration. For technosignature detection, the primary motivation for COSMIC, two types of real-time calibration are implemented: delay calibration and gain calibration. Currently, we are not conducting any bandpass or amplitude calibrations in real time, although the correlated data of amplitude and bandpass calibrators are saved in case we need to calibrate the bandpass responses and amplitudes following the detection of extraterrestrial intelligence (ETI). However, amplitude equalization is conducted in the FPGAs at the time of configuration to normalize the signal levels from all antennas. This is done to ensure that the beamformed output is not dominated by signals from certain antennas. We have determined that through this process, we obtain an equivalent system flux density within 10% of the value reported by the NRAO for the VLA.

3.3.1. Delay Calibration

3.3.2. Fixed Delays

3.3.3. Delay Tracking and the Geometric Delay Model

As the Earth rotates, the projection of each baseline changes with respect to the phase tracking center. The phase tracking centers (the placement of the coherent beams) are chosen by the target selector (Section 3.5) to be within the full width at half-maximum (FWHM) of the primary beam.

The phase-centered controller can operate in two modes. In mode 1, the R.A. and decl. values are sent out at specific intervals, usually corresponding to updates from the multicast system through Redis, and the delay engine is left to derive phase solutions for the current pointing location. In mode 2, at time intervals that align with the GUPPI raw file boundaries (8 s for the B, C, and D configurations of the VLA and 2 s for the A configuration) during data recording, the phase-center controller updates the phase pointing such that each file is phased to only a single pointing location. This precision is achieved by associating load times (t_load) with the pointing coordinates and sending the pointing coordinates and load times out several seconds before t_load.

The delay model receives the pointing positions and associated t_load values via a Redis channel and, using the positions of the antennas along with the time of reception t_r, derives quadratic coefficients C₀, C₁, and C₂ such that the delay at t_r is calculated as follows:

$$\begin{eqnarray}&&\mathrm{Delay}({t}_{{\rm{r}}})={C}_{0}+{C}_{1}\times {t}_{{\rm{r}}}+\displaystyle \frac{{C}_{2}}{2}\times {t}_{{\rm{r}}}^{2}.\end{eqnarray} \tag{ 1 }$$

The delay model then sends the coefficients C₀, C₁, and C₂ out to each of the F-engines via separate Redis channels, along with the current scans of sideband and center frequency, t_r, and the received t_load.

Finally, the F-engines receive this information, with which they perform a quadratic interpolation to the provided t_load for both the delay,

$$\begin{eqnarray}&&\mathrm{Delay}({t}_{\mathrm{load}})={\tau }_{s}+{C}_{0}+{C}_{1}\times {t}_{\mathrm{diff}}+\displaystyle \frac{{C}_{2}}{2}\times {t}_{\mathrm{diff}}^{2},\end{eqnarray} \tag{ 2 }$$

and the delay rate,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}(\mathrm{Delay}({t}_{\mathrm{load}}))={C}_{1}+{C}_{2}\times {t}_{\mathrm{diff}}^{2},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{t}_{\mathrm{diff}}={t}_{\mathrm{load}}-{t}_{{\rm{r}}},\end{eqnarray} \tag{ 4 }$$

and τ_s is the fixed calibration delay derived from a GUPPI raw file for the FPGA streams ³⁰ in seconds.

Delay(t_load) has units of nanoseconds and is a decimal float value. The integer part is separated and applied as a coarse delay (prechannelization), while the fractional part is applied as a fine delay (postchannelization). $\tfrac{d}{{dt}}(\mathrm{Delay}({t}_{\mathrm{load}}))$ has units of ns s⁻¹ and is also loaded into the F-engines for the FPGAs to perform linear interpolation from the loaded fractional delay.

It is imperative to continuously update the delays on the F-engines, as linear interpolation does not provide the accuracy required to track sources accurately, and the fractional delays in the F-engines are susceptible to information overflow. For this reason, for each antenna on each F-engine node, a Python thread is spawned at 0.5 s intervals to quadratically reinterpolate the delay and delay rates to t_load (if no specific t_load is provided, t_load is set to 0.5 s into the future) and load the results into the F-engines at t_load.

Simultaneously, a phase and phase rate are derived from the computed delay and delay rate, which are also applied postchannelization, as follows:

$$\begin{eqnarray}&&{\rm{\Phi }}({t}_{\mathrm{load}})=2\pi \times {{lo}}_{\mathrm{cf}}\times {{sb}}_{s}\times {{lo}}_{\mathrm{cf}}\times {{sslo}}_{s}\times \mathrm{delay}({t}_{\mathrm{load}})\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}({\rm{\Phi }}({t}_{\mathrm{load}}))=2\pi \times {{sb}}_{s}\times {{sslo}}_{s}\times \displaystyle \frac{d}{{dt}}(\mathrm{delay}({t}_{\mathrm{load}})),\end{eqnarray} \tag{ 6 }$$

where sb_s is the per-stream sideband value, sslo_s is the per-stream effective LO, and lo_cf = 2π × sb_s ∗ f_shift is a phase-correction factor that needs to be applied due to the f_shift mixing performed for each FPGA stream before fine channelization.

3.3.4. Gain Calibration

The delay errors are treated as fixed values when generating the output of the calibration pipeline during delay tracking, but due to uncertainties (of less than 10 ns) in the fixed delay values, there are still some residual uncorrected delays. The first priority of COSMIC is the detection of technosignature signals, and the only calibration that we consider after delay calibration is the calibration of the complex gain; no absolute flux density calibration is performed.

Gain calibration is an antenna-based correction that accounts for time-varying factors associated with the instrument and the atmosphere. If these factors remain uncorrected, the differences in the gains for each antenna can impact the phasing of the incoming voltage streams, leading to decorrelation of the sum of the antenna signals.

To conduct the gain calibration, we use a modified version of the sdmpy ³¹ package, which was initially written to calibrate the VLA Science Data Model data sets for the realfast commensal system (Law et al. 2018). In this version of the package, the calibration utility in the sdmpy software has been modified to work with the COSMIC UVH5 data format.

The gain calibration solutions for each stream are written out as a Python dictionary. The derive_gains method in the calibrate_uvh5 class of calibrate_uvh5.py ³² is used to derive the antenna gains. Gain calibration is carried out independently on each 32 MHz bandwidth (32 × 1 MHz coarse channels) of the correlated UVH5-formatted data (for each LO tuning), which are distributed across multiple GPU compute nodes. This is done such that each of the 44 compute nodes (two compute nodes per GPU server) has 32 MHz of data to process. The resultant gain dictionaries are sent to the head node, where the gains across the multiple 32 MHz subbands are combined along the frequency axis for further processing. After all gain dictionaries have been collected for the 1 MHz coarse channels for both LO tunings, all polarizations, and each antenna, the head node sorts and orders the complex gains into a Python numpy matrix with dimensions of n_pol × n_tunings and n_frequencies.

This resulting gain matrix is fed into a calibration kernel that can perform either a linear or Fourier interpolation to calculate the residual delays and phases for each antenna and polarization. The calibration phases previously loaded into the FPGAs are also subtracted from the received gain matrix to calculate the new delay residuals and phases that are provided to the F-engines.

A plot of the residual delays and the per-antenna phase solutions for both IFs is also uploaded to a Slack ³³ interface channel in real time to allow scrutiny by all collaborating scientists. The amplitude plots are also uploaded for review even though no amplitude calibration is performed.

Figure 9 shows an example plot of the phases of the gain solutions plotted as a function of frequency for each antenna that is sent to the designated Slack channel for a single calibration observation. The current COSMIC system collects data only in the second half of each IF, spanning 512 MHz, ³⁴ due to computational resource limitations, which is why some of the band is grayed out and zero (see Section 5.1 for a discussion of upgrades to the recording bandwidth). The plot shows consistent flat phases close to zero for each operational antenna, as expected for well-behaved phase calibration. We note that radio frequency interference (RFI) from sources external or internal to the antenna data stream can introduce bad phase solutions, especially near the region of 1–2 GHz, which is heavily affected by satellite RFI.

Zoom In Zoom Out Reset image size

Along with the real-time calibration, we also calculate certain statistical parameters to determine the quality of the calibration.

• 1.  
  The inverse FFT of the correlated spectra for each antenna with respect to a reference antenna is used to estimate the residual delay. The signal-to-noise ratio (S/N) of the delay peak serves as a measure of the coherence between voltages for the corresponding baseline. A low S/N value for a baseline indicates decoherence between the voltages. If the S/N value is 4 or less, then the calibration is not applied to the F-engines.
• 2.  
  The standard deviation of the phases is calculated to understand their spread.
• 3.  
  An antenna grade,

  $$\begin{eqnarray}&&{G}_{\mathrm{ant}}=\displaystyle \frac{\left|{\sum }_{\nu =0}^{\nu =n\_\mathrm{chan}}\mathrm{Gains}(\nu ,\mathrm{ant})\right|}{{\sum }_{\nu =0}^{\nu =n\_\mathrm{chan}}| \mathrm{Gains}(\nu ,\mathrm{ant})| },\end{eqnarray} \tag{ 7 }$$

  is calculated for each stream. A value of G_ant = 1 indicates a flat zero phase across frequencies and thus a good calibration. In contrast, G_ant = 0 results from a large variation in phase across frequencies, which indicates a poor calibration.
• 4.  
  A frequency channel grade,

  $$\begin{eqnarray}&&{G}_{\nu }=\displaystyle \frac{\left|{\sum }_{\mathrm{ant}=0}^{\mathrm{ant}=n\_\mathrm{ant}}\mathrm{Gains}(\nu ,\mathrm{ant})\right|}{{\sum }_{\mathrm{ant}=0}^{\mathrm{ant}=n\_\mathrm{ant}}| \mathrm{Gains}(\nu ,\mathrm{ant})| },\end{eqnarray} \tag{ 8 }$$

  is calculated for each stream. A value of G_ν = 1 indicates that the antennas are properly phased to the intended position, implying a good calibration. G_ν = 0 indicates phase differences between the antennas resulting from a poor calibration.
• 5.  
  An overall grade,

  $$\begin{eqnarray}&&G=\displaystyle \frac{\left|{\sum }_{\mathrm{ant},\nu =0}^{\mathrm{ant},\nu =n\_\mathrm{ant},n\_s,n\_\mathrm{chan}}\mathrm{Gains}(\nu ,\mathrm{ant},s)\right|}{{\sum }_{\mathrm{ant},\nu =0}^{\mathrm{ant},\nu =n\_\mathrm{ant},n\_s,n\_\mathrm{chan}}| \mathrm{Gains}(\nu ,\mathrm{ant},s)| },\end{eqnarray} \tag{ 9 }$$

  is calculated across the full gain matrix. This overall grade tends toward unity when the calibration recording is phased to the intended position across all antennas, streams, and frequencies. The overall grade tends toward zero when the calibration recording is incorrectly phased across one or more of these variables.

The start of VLA observations is triggered when the antennas start to slew toward the designated source. However, the antenna information received through Redis will show that the antennas are not "on source"; therefore, they will be flagged accordingly until they arrive at the target location. Additionally, any antennas that are offline (stowed) and therefore never "on source" are flagged prior to correlation or raw data recording. The COSMIC system also flags any antenna that has unstable DTS information and thus does not record data from antennas with variable power.

RFI consists of human-made signals not related to the signal being sought, often generated by terrestrial technology or satellites. In our search for technosignatures, we aim to eliminate RFI—both external and intrinsic to the processing pipeline—to seek out signals associated with extraterrestrial technology. Although RFI detection and excision is a critical step in determining the difference between signals of terrestrial and extraterrestrial origin, the current COSMIC pipeline does not implement real-time RFI excision; therefore, early data can be used to evaluate the RFI environment and experiment with potential excision algorithms on representative data. Over the many years of operation of the VLA, the observatory has gathered a vast amount of data regarding the RFI problem. Therefore, when we initially evaluate the data for scientific merit, we rely on these historical data ³⁵ from the facility to avoid the regions of the observation band that are significantly affected by RFI.

We note that the RFI environment for extremely high-resolution spectroscopy, such as the technosignature data of interest to COSMIC, can be much more complex than the environment for the signals typically acquired for general observatory purposes. It is envisioned that in the future, a database of RFI will be created within the COSMIC data processing pipeline to enable us to mask these interfering signals before beginning a drift search. However, we can presently identify regions of strong RFI in which signals should be disregarded on a first pass through the data. An example of power spectra from one night of observation is shown in Figure 10.

Zoom In Zoom Out Reset image size

The NRAO's multicast system continuously provides the current pointing for each observation, from which the target selector software ³⁶ calculates the primary FoV ³⁷ (see Figure 16). Using this FoV, the target selector determines which stars are available for observation and calculates an optimized list of the coordinates of the highest interest at which to form coherent beams. The target selector draws these target stars and other objects (falling within the primary FoV) from several databases, including a 32 million star sample derived from Gaia Data Release II (Czech et al. 2021) and the Breakthrough Listen Exotica Catalog (Lacki et al. 2021). The targets are ranked in priority based on their distance. The goal is to search as many unique stars as possible, prioritizing our closest neighbors. The target coordinates are fed into a special HDF5-based "beamformer recipe file (BFR5)" used by the beamformer to form coherent beams.

COSMIC forms at least five coherent tracking beams, depending on the available computing resources and the amount of time spent recording. The idea is to form all of the beams at once over the full recorded FoV. The raw voltage data are stored in GUPPI raw binary files, which contain blocks of data and header information.

All of the time samples in a GUPPI raw file are acquired to perform an FFT to finely channelize the recorded data stream from the 1 MHz coarse channels to ∼8 Hz fine frequency channels. The beamforming phase coefficients are calculated as follows:

$$\begin{eqnarray}&&w(\phi )={e}^{j\phi },\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&\phi =2\pi f\tau +{\psi }_{f},\end{eqnarray} \tag{ 11 }$$

f is the center frequency of the coarse (1 MHz) channel, τ is the delay relative to the reference antenna, and ψ_f is the phase calibration solution for a particular coarse channel.

The BFR5 "beamformer recipe file" provides the list of target source coordinates and other relevant information to inform the calculation of the phase coefficients for the beam formation process, which is performed using the Breakthrough Listen Accelerated DSP Engine (BLADE; L. Cruz et al. 2023, in preparation), ³⁸ a software suite designed for the Allen Telescope Array in California, USA, and modified for use in the COSMIC pipeline. Currently, COSMIC performs a 131,027-point FFT and forms five coherent beams plus an incoherent beam. The resulting total of 64 time samples per spectrum is channelized to 7.6294 Hz given a 1 MHz coarse channel. Only one set of phase coefficients is used for beam formation for the recorded time span of 8.388 s (when the VLA is in the B, C, or D configuration) or ∼2 s (when the VLA is in the A configuration). However, these values can be flexibly tuned depending on the situation. This approach results in an average beamforming computation time on the GPU of approximately 50 s from file ingestion to the generation of output from BLADE (including memory transfers). Improvements to the BLADE software package to decrease this processing time are ongoing.

At the output of the beamformer, the raw voltages of the coherent beams are converted into power, and the polarizations are summed to a pseudo-Stokes I total intensity. These processes are all implemented on the GPU with CUDA. The output is then passed to the technosignature search algorithm discussed in Section 3.6.

To assess the performance of the beamformer and ensure that the phases and delays are appropriately accounted for, as well as to verify the general operation of the beamformer, a bright class II methanol maser was observed at 6.7 GHz. We obtained test observations of W51M, which, according to Etoka et al. (2012), has a brightness of ∼250 Jy and is close to a point source at the VLA B-configuration ³⁹ resolution of ∼2''. We imaged the VLA WIDAR-collected data using CASA (McMullin et al. 2007) to find the exact location of the source and independently verify the topocentric frequency of emission in comparison to COSMIC, as both data sets were recorded simultaneously. We determined the source position to lie at R.A.(J2000) = 19^h23^m4395 and decl.(J2000) = 14°30^m3434, with the topocentric frequency of the peak emission being at 6667.986 MHz. By looking at the time-averaged power spectrum of the coherent beamformed data from COSMIC, we found that the frequencies matched and that the signal was detected at a velocity that was well matched with the published results. Thus, by forming coherent beams at a 025 separation around the WIDAR position for the source, we achieved a pointing accuracy within half of the VLA point-spread function.

To test the sensitivity of the expected flux density recovery, with the VLA in the A configuration, we obtained two observations of the class II methanol maser associated with W3OH at R.A.(J2000) = 02^h27^m038192 and decl.(J2000) = 61°52^m25230. In one observation, the source was near the center of the primary beam (boresight), and it was at the half-power point of the primary beam in the other, approximately 2' from the phase center for the A configuration, each with ∼8 s of recorded data from the source. These two observations were recorded with 22 operational VLA antennas (the data were taken during maintenance, so the full array was not available), while COSMIC was configured with the raw data output and save mode turned on but the technosignature search turned off. As observed in Figure 11, the maser at the FWHM of the primary beam exhibits half the power level of the maser at the boresight (a ∼2 dB difference), as expected.

Zoom In Zoom Out Reset image size

3.7.1. Beamformer Efficiency

To measure the beamformer efficiency, observations of how the S/N changes with the addition of antennas were computed as an overall measure of the correlated sky noise (Kudale & Chengalur 2017). For this analysis, we formed coherent beams toward the W51M methanol maser, a strong emitter for which perturbations in the S/N can be reduced by using the frequency channel with the most intense signal and the time-averaged power over a 5 minute observation. The observations were calibrated as described in Section 3), with an evaluation of the phase stability before and after the maser observation. Using the BLADE beamforming code, we incrementally added 3–21 antennas in sets of three antennas in each beamforming round. The noise was calculated by taking the mean of 4000 RFI-free channels (each 8 Hz wide, for a total of 32 kHz) from the same 1 MHz coarse channel as the maser emission line but not including the maser emission. The S/N was then calculated by subtracting the peak power of the source and the noise, both in units of dB.

To determine the S/N improvement, which is shown on the y-axis in Figure 12, all values were then subtracted from the value obtained from only three antennas (the first data point). The ideal S/N improvement is calculated as a linear increase from 0 to 7.8 dB (10∗log10(18) − 10∗log10(3) = 7.8 dB), representing the theoretical expectation for a 100% operational system. The plot shows that as more antennas are sequentially added, the S/N increases linearly, as expected. However, the total power is approximately 20% lower than expected if all signal paths for all antennas are functioning at the theoretical total efficiency, with some contribution from phase calibration errors. This is within the expected tolerances for the system.

Zoom In Zoom Out Reset image size

3.7.2. Coherent versus Incoherent Beams

A coherent beam, in the receiving paradigm, is produced by altering the phase on the elements of an array such that parallel-plane electromagnetic radiation emanating from a particular source combines constructively (Section 3.7). This method trades off beam size for sensitivity when compared to an incoherent beam. An incoherent beam, produced by summing the total intensity of all receiving elements (disregarding phase information by definition), exposes sources across the FoV of each telescope, which is still important when performing a blind SETI search.

We thus compared the coherent and incoherent beam formation processes to ensure that the expected decrease in sensitivity for the W51M class II methanol maser signal was observed. The beamforming process was conducted on the same observation with the same integration time and number of antennas, with a phase calibration completed as described in Section 3.3. Because the data were not flux density calibrated, we evaluated the ratio based on the normalized power of the two signal strengths for the peak emission from the source.

For the incoherent summation, the noise (σ), source flux density (S_src), and equivalent system flux density (S_sys) are related as follows:

• 1.  
  S_src/S_sys is constant, and
• 2.  
  S_src/σ scales as $\sqrt{N}$, where N is the number of antennas.

For the coherent summation, the following relations hold:

• 1.  
  S_src/S_sys scales as N, and
• 2.  
  S_src/σ also scales as N.

Therefore, if the S/N is defined as S_src/σ, the coherent and incoherent summations on the same data set can be related by

$$\begin{eqnarray}&&\mathrm{Norm}\,{\mathrm{Power}}_{\mathrm{coherent}}/\mathrm{Norm}\,{\mathrm{Power}}_{\mathrm{incoherent}}=N.\end{eqnarray} \tag{ 12 }$$

Here, the Norm Power value is calculated by subtracting the power for a frequency range that is far from the source of the signal from the power for the spectral frequency range of the signal and dividing by the off-source power ((on–off)/off). For an observation of the W51M methanol maser, we performed coherent and incoherent beamforming operations on the same data set with 23 operational antennas. As shown in Figure 13, we found a maximum ratio of 20.5, which suggests an 89% efficiency, similar to the value computed in Figure 12. As seen from the phase calibration plots for this observation, four antennas had problems with their phase delays, contributing to this reduced efficiency. This is considered an acceptable level for an autonomous real-time system, but further improvements will be sought.

Zoom In Zoom Out Reset image size

Hypothesized radio technosignatures take many different forms in terms of both the frequency width of the signal and the cadence of the signal. Each signal form requires a different search approach, and multiple methodologies have been explored over the years. Many of these approaches and signal forms were discussed in Ekers et al. (2002), but other techniques have been discussed since that time (i.e., Houston et al. 2021; Luan et al. 2023; Suresh et al. 2023). However, publications on SETI have thus far placed a strong focus on drifting narrow frequency band (Hz-wide) signals (e.g., Enriquez et al. 2017; Price et al. 2020; Sheikh et al. 2021; Ma et al. 2023).

COSMIC, at least initially, is designed to passively look for electromagnetic radiation at radio frequencies of 0.75–50 GHz, which are observable by the standard VLA signal chain. ⁴⁰ To differentiate from the natural background of astrophysical emission, we search for continuous-wave (monochromatic) frequency-drifting signals in 1–8 Hz wide channels, where the drift comes from the difference in acceleration between the emitter and the receiver (the VLA telescope for COSMIC; Li et al. 2022).

As an example of a narrowband technological signal and COSMIC's ability to detect such signals, we observed the Voyager 1 downlink signal sent to Earth at 8.4 GHz and received by the VLA (Figure 14). Voyager 1, launched on 1977 September 5, is currently 159 au from Earth within the constellation of Ophiuchus (within the Oort Cloud). The coordinates were obtained from the Horizons astroquery database, and the expected frequency and Doppler shift were calculated using the Python astropy package. COSMIC detected the signal in both the incoherent and coherent beams near the expected frequency. The coordinates to form the coherent beam were obtained based on images of Voyager 1 by correlating the raw voltages. Figure 14 shows the signal that was detected when a coherent beam was formed at the corresponding location. The autonomous real-time search pipeline also detected the signal and recorded the information in the database (a process described in more detail in Section 3.9).

Zoom In Zoom Out Reset image size

The pipeline running in COSMIC uses the Doppler acceleration search algorithm seticore, ⁴¹ which is run on the data created via coherent and incoherent beamforming by means of the BLADE software package. seticore is a GPU implementation of TurboSETI (Enriquez et al. 2017), which is a Taylor tree search algorithm (Sheikh et al. 2019) used in many narrowband signal searches over the past 5 yr. The user-provided input parameters for seticore include an S/N threshold, which was set to 10 for the first 6 months of COSMIC's operation, and a Doppler drift rate, which is currently set to ±50 Hz s⁻¹. However, these parameters can be changed within the observational YAML file.

At the end of a search, a series of information regarding the signals found (hits) is logged in an SQL database, and raw voltage "stamp" files for each antenna are saved in the storage nodes under a directory name that is a combination of the NRAO scheduling block ID, the project ID, and the scan ID. These stamp files are segments of the raw antenna voltages around the brightest signals, containing at least 200 Hz of frequency data and all time samples available in the recording. These voltages can later be plotted using seticore and/or correlated and imaged near signals of interest to confirm the signal characteristics.

In the future, machine-learning approaches similar to the search described in Ma et al. (2023) could be applied to COSMIC-generated data. However, currently, the data are manually reviewed through statistical analyses, especially as we learn more about these data products. These analyses will be discussed in a future paper.

The COSMIC pipeline is designed to identify signals of interest by detecting narrowband emission ("hits"), and then various postprocessing pipelines may be created to filter those hits. In particular, two consecutive filters can be applied, a spatial filter (to determine whether a signal appears in one coherent beam but not in the others) and a drift rate filter (to determine whether the drift rate is between ±50 Hz s⁻¹ but not identically zero, indicating a source that is accelerating relative to the receiver). In fact, the maximum drift rate is flexible and can be changed at any time depending on scientific interest. However, the current value of ±50 Hz s⁻¹ is set larger than that for a typical de-Doppler search for technosignatures (Sheikh et al. 2019) due to our increased computing capabilities, allowing us to widen the restrictions on where we assume a transmitter to be located and how fast it is assumed to be rotating while also covering the drift rates determined by Li et al. (2022) for the planetary systems nearest to our solar system.

Any real astronomical signal of interest would be a point source centered within a single coherent beam or with minimal leakage into closely packed nearby coherent beams. Terrestrial signals would, in contrast, "bleed" into multiple beams within a single FoV. The Doppler search kernel produces "postage stamps" of raw voltage files for at least 200 Hz of bandwidth around the signal of interest and stores them in a series of folders for further evaluation, while the metadata regarding the hits contained within these stamp files are held in an SQL database. An example of the information that can potentially be obtained from such postage stamps is shown in Figure 15, where a drifting carrier signal is detected in the waterfall plots for each antenna. We would expect an astronomical signal to be detected by all operational antennas.

Zoom In Zoom Out Reset image size

Given access to both coherent and incoherent beam data, we can assume that any signal present in a coherent beam will also appear in the incoherent beam if the S/N is sufficiently high. Thus, when we compare the signal from the incoherent beam with the signal in the coherent beam, the incoherent beam should contain a signal, as follows:

$$\begin{eqnarray}&&{\rm{coherent}}{\rm{\_}}{\rm{S}}/{\rm{N}}\gt {\rm{sqrt}}({N}_{{\rm{antennas}}})\ast {\rm{incoherent}}{\rm{\_}}{\rm{S}}/{\rm{N}}.\end{eqnarray} \tag{ 13 }$$

If this is not true, or if a signal is detected in all coherent beams, then the signal is likely to be RFI and should be ignored as an ETI signal of interest. However, as we investigate all information being recorded by COSMIC and processed through BLADE and seticore, all hits are currently recorded into a database without filtering. Instead, signal rejection based on this type of comparison between the coherent and incoherent beams is performed is a postprocessing step, and the details will be provided in future work.

For any signal that passes these criteria, a check for known RFI or astronomical sources is completed, and a visual inspection of the characteristics of the drift profile is performed by means of waterfall plots. A total of 10 steps are adopted by Breakthrough Listen to verify a signal, as explained in Sheikh et al. (2021), including checking the telescope and digital signal processing system for potential problems, examining the drift rate evolution and evaluation of the signal of interest, and searching for other instances of the signal of interest in archival data. If a signal passes all of these criteria, reobservation of the region from which the signal originated, with either the VLA or other telescopes, is critical as a check for redetection. If such a signal were to be confirmed, this would be revolutionary, allowing the field to blossom into a new era of understanding the possible opportunities for life outside our own planet.

Accordingly, we have successful target-of-opportunity proposals in place with the VLA, the NRAO Very Long Baseline Array, the Allen Telescope Array, and the Parkes 64 m telescope in Australia. Additionally, through collaboration with Breakthrough Listen, we can utilize available time on the GBT. These collaborative efforts allow us to cover most, if not all, available frequencies that can be recorded with COSMIC and to cover the sky while utilizing telescopes with independent RFI environments and data collection systems. However, we make no assumption that a signal detected from an astronomical origin will be repeated. A repeated signal detected by another facility and/or hardware setup would be an ideal result.

COSMIC, as a system, spans many compute nodes, devices, and services and is reliant on all of them working in sync with the pipeline. Throughout the development of the software system, it was imperative for each process to log its state, errors, and interactivity, and this will continue to be important as additional functionalities are incorporated into COSMIC. These process logs allow faster debugging and shorten the time during which the system is not processing data. Where software services fail due to firmware or hardware issues, it becomes apparent that monitoring and logging the hardware and firmware status is also important.

The results of COSMIC delay (Section 3.3.1), F-engine (Section 2.2), and calibration (Section 3.3) state monitoring (Section 3.4) are stored in an Influx database, ⁴² which supports the nanosecond-scale logging necessary when debugging delay tracking and F-engine problems. Data are retained in this error database for only 30 days, since the high read/write polling rates require it to be situated on the head node main drive, which has limited memory.

The COSMIC GPU compute node information is polled via Prometheus ⁴³ exporters but is not stored for longer than a few hours due to storage resource limitations. systemd ⁴⁴ services are responsible for antenna control monitoring, while the delay engine, calibration, and hashpipe are also polled by different Prometheus exporters, as failure of these services will cause the system to stop commensal observations. It is not necessary to retain the compute node history, as this polling is primarily intended for alerting of compute problems requiring immediate attention.

Both the Influx content and Prometheus content are displayed via a Grafana ⁴⁵ dashboard. From this dashboard, an overview of the antenna, F-engine, delay, calibration, and compute node processing states may be accessed. In addition, Grafana allows custom alerts to be configured to trigger on specific thresholds for any state change, and messages can be sent to a designated Slack channel.

The COSMIC system is flexibly designed to allow for future upgrades to hardware and software. By adding additional compute clusters, we can coherently beamform toward more targets, thus covering a larger number of stars or other astronomical sources of interest or increasing the total bandwidth searched. With our current computing capabilities, we are searching five coherent beams when observing along with VLASS and between 4 and 64 beams for other observations, plus one incoherent beam. We can also record and search up to 1 GHz of total bandwidth during VLASS and up to 1.4 GHz of bandwidth along with other observations, although at the cost of reducing the number of coherent beams. In the future, we could increase the number of beams and the total bandwidth recorded and searched by expanding the computing infrastructure.

Currently, COSMIC is designed to work with the VLA's 8 bit system with a maximum total simultaneous bandwidth of 2 GHz (Perley et al. 2011). Future software upgrades will allow data to be recorded from the VLA when the 3 bit system is turned on, which offers 8 GHz of simultaneous bandwidth. As mentioned in Section 2, the channelization setup through the VLA control system does not impact COSMIC, so we can maintain flexibility regarding the resolutions in frequency and time with which we generate and record data. This flexibility will allow us to broaden the scientific output of COSMIC beyond SETI.

The current system is also designed to search for signals using the Stokes I total intensity. However, COSMIC ingests all four polarization data products and thus, with a software upgrade, could search for linear and circular polarization signals.

As with other programs where SETI is the main motivation for the design concept, other scientific endeavors can also be pursued with either a copy of the data or a change to the processing setup. In particular, COSMIC can enable searches for transients with submillisecond temporal resolution, such as fast radio bursts (e.g., Diermyer 2021; Faber et al. 2021), or the data can enable spectral line science and axionic dark matter searches (e.g., Foster et al. 2022). Moreover, with COSMIC's flexible design, many other scientific capabilities may be explored.

Another benefit of designing COSMIC as an Ethernet-based system is the ability to set up multicasting Ethernet technology for multimode commensal recording and processing systems. This means that other commensal systems could tap into the COSMIC digital processing rack to create other real-time scientific data outputs. Such initiatives could include offering additional functionalities for the VLA Low Band Ionospheric and Transient Experiment (Clarke et al. 2016), realfast (Law et al. 2015, 2018), or other commensal systems currently operating on the VLA.