# Prospects for Wideband VLBI Correlation in the Cloud

Ajay Gill  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada  Lindy Blackburn  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Arash Roshanineshat  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Department of Electrical and Computer Engineering, University of Arizona, 1230 E. Speedway Blvd., P.O. Box 210104, Tucson, AZ.  Department of Astronomy and Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065  Chi-Kwan Chan  Department of Astronomy and Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065  Sheperd S. Doeleman  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Michael D. Johnson  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Alexander W. Raymond  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Jonathan Weintroub  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

This paper proposes a cloud architecture for the correlation of wide bandwidth
VLBI data. Cloud correlation facilitates processing of entire experiments in
parallel using flexibly allocated and practically unlimited compute resources.
This approach offers a potential improvement over dedicated correlation
clusters, which are constrained by a fixed number of installed processor nodes
and playback units. Additionally, cloud storage offers an alternative to
maintaining a fleet of hard-disk drives that might be utilized intermittently.
Here, we describe benchmarks of VLBI correlation using the DiFX-2.5.2 software
on the Google Cloud Platform to assess cloud-based correlation performance. In
our analysis, the number of virtual CPUs (vCPUs) per Virtual Machine was
varied to determine the optimum configuration of cloud resources. The number
of stations was varied to determine the scaling of correlation time with VLBI
arrays of different sizes. Data transfer rates from Google Cloud Storage to
the Virtual Machines performing the correlation were also measured. Based on
the results, we present an example cloud correlation configuration. Current
cloud service and equipment pricing data is used to compile cost estimates
allowing an approximate economic comparison between cloud and cluster
processing. We note that the economic comparisons are based on cost figures
which are a moving target, and are highly dependent on factors such as the
utilization of cluster and media, which are a challenge to estimate. Our model
suggests that shifting to the cloud is an alternative path for high data rate,
low duty cycle wideband VLBI correlation that should continue to be explored.
In the production phase of VLBI correlation, the cloud has the potential to
significantly reduce data processing times and allow the processing of more
science experiments in a given year for the petabyte-scale data sets
increasingly common in both astronomy and geodesy VLBI applications.

##  1 Introduction

The technique of Very Long Baseline Interferometry (VLBI) links together many
radio telescopes to make high angular resolution observations of distant
astronomical sources, or in geodesy to precisely measure the shape and
orientation of the Earth (Clark et al.,, 1967; Moran et al.,, 1967; Thompson
et al.,, 2017). The sensitivity for VLBI instruments is determined by the
collecting area (Matthews et al.,, 2017), by the integration time, and by the
recorded bandwidth. For a given antenna size and atmospheric coherence time,
the sensitivity improves as $B^{1/2}$, where $B$ is the total recorded
bandwidth (Whitney et al.,, 2013). The millimeter-wave VLBI array known as the
Event Horizon Telescope (EHT), for example, has steadily increased its
bandwidth from less than 0.5 to 16 GHz (Event Horizon Telescope Collaboration
et al., 2019b, ) or a data recording rate for 2-bit samples of 64 Gbps
(gigabits-per-second). Such a large bandwidth yields petabyte-scale
recordings, which is an enormous data volume and therefore a challenge to
correlate.

In typical VLBI arrays, each station simultaneously observes the astronomical
target and records the data onto hard-disk drives that are then sent to a
correlation facility, where hardware such as Field-Programmable Gate Arrays
(FPGAs), Application-Specific Integrated Circuits (ASICs), Central Processing
Units (CPUs), and/or Graphical Processing Units (GPUs) perform the
correlation. Instead of transporting disk drives, some arrays transfer the
data directly over network with a parallel file system from a station or a
nearby staging area to a storage pool at the correlation facility. In this
study, we primarily consider disconnected-element interferometers with
stations separated by long baselines ($\gtrsim$ 100 km), where the data are
recorded on disk drives and transported to the correlation facility. Doubling
the received bandwidth to improve sensitivity means that the recording rate
and the data volume also doubles for a given integration time. Advances in
high-speed recording made possible by industry-supported enhancements in
analog-to-digital conversion (Patel et al.,, 2014) and improvements in network
data transfer mean that data volumes for modern VLBI arrays operating at wide
bandwidths can be many petabytes per month. The trend to use more bits per
sample to improve the signal-to-noise ratio and RFI robustness also translates
to increasing data rates and data volumes.

To address the difficulty of processing these large data sets, we consider the
possibility of using decentralized computational resources available over the
Internet. These resources, referred to as the cloud, enable the realization of
Virtual Machines (VMs) consisting of multiple virtual Central Processing Units
(vCPUs). Cloud-based VMs could process VLBI data in parallel at speeds
commensurate with the increasing recording rates. A benefit of the cloud
approach to VLBI correlation is the ability to process more science
experiments in a given year. A large dedicated supercluster, where the data is
transferred in parallel over network from the stations or a nearby staging
area provides similar improvements in data processing times, but it is not
always possible for every research group to access such a supercluster.

There have been previous studies of distributed computing for VLBI. In early
work, (van Langevelde,, 2006) employed a precursor to the cloud called the
grid in a project called FABRIC, where distributed computing was applied to
the Merlin array observing in an e-VLBI mode with correlation in real time via
the Internet. (Weston et al.,, 2017) performed benchmarks for VLBI correlation
using the Catalyst IT Ltd commercial cloud in New Zealand. They found that
performing correlation on the cloud provides identical fringe outputs compared
to cluster correlation, and they benchmarked a speed-up of a factor of $\sim$
2.5, correlating data rather narrower in bandwidth compared to the present
study. On the software side, CorrelX is a correlation software package
currently under development to leverage the cloud for VLBI
correlation111https://github.com/MITHaystack/CorrelX. After experimenting with
CorrelX and consulting with its authors, we elected instead to use the widely
used DiFX-2.5.2 VLBI software correlator (Deller et al.,, 2011) exclusively
for this study.

In this paper, we propose an architecture for shifting wideband VLBI
correlation to the cloud. To assess the efficiency of cloud correlation,
benchmarks across a parameter space of possible configurations were carried
out on the Google Cloud Platform (GCP)222https://cloud.google.com. Performance
and scaling trends were characterized by varying the number of vCPUs per VM,
the number of stations in the array, and testing different data transfer rates
from Google Cloud Storage (GCS) to the VMs. Although the GCP was used for this
study, the proposed cloud correlation architecture is generalizable to other
major cloud platforms such as Amazon AWS333https://aws.amazon.com and
Jetstream444https://jetstream-cloud.org (Stewart et al.,, 2015; Towns et al.,,
2014).

##  2 VLBI applications

The two major applications of VLBI are astronomical research and geodesy.

###  2.1 Astronomy

The angular resolution of a VLBI radio interferometer is
$\theta\simeq\lambda\,/{B_{\rm max}}$ in radians or $\theta\simeq 2\times
10^{5}\lambda\,/{B_{\rm max}}$ in arcseconds, where $\lambda$ is the
observation wavelength, and ${B_{\rm max}}$ is the length of the maximum
projected baseline. At a wavelength of 2 cm with intercontinental baselines of
10,000 km, an angular resolution of 0.4 mas can be achieved, whereas at the
shorter wavelength of 1.3 mm with baselines of 10,000 km, an angular
resolution of 20 $\mu$as is possible (Event Horizon Telescope Collaboration et
al., 2019c, ), a resolution that is unparalleled in Earth-based diffraction
limited observations.

The massive compute and storage resources provided by the cloud could handle
the increased data rates and data volumes resulting from wider bandwidth VLBI
observations. The scientific opportunities enabled by expanding the bandwidth
of VLBI have been considered in various studies (Doeleman et al.,, 2004;
Walker et al.,, 2007; Fish et al.,, 2013). This includes microarcsecond
astrometry of high mass star forming regions within the Milky Way (Reid et
al.,, 2009), stellar mass black holes (Miller-Jones et al.,, 2009; Reid et
al.,, 2011), pulsars and tests of general relativity (Deller et al.,, 2007),
and satellites within the Local Group. Increasing VLBI bandwidths can also
enhance the number of detectable targets, transients, and astrometric
reference sources, since the number of detectable sources scales with
bandwidth as N $\propto$ B3/4 (Reid and Honma,, 2014). Improvements in imaging
gravitational lens systems to study the properties of the relativistic jets in
distant Active Galactic Nuclei (Koopmans and Treu,, 2002) and small scale
structure near the cores of lensing galaxies can be expected (Mao et al.,,
2001; Winn et al.,, 2004; Zackrisson and Riehm,, 2010). More sensitive
observations could also allow constraining the properties of supernovae
(Bartel et al.,, 2000) and gamma ray bursts (Taylor et al.,, 1998; Bietenholz
et al.,, 2014; Abbott et al.,, 2017) by resolving their afterglows as well as
constraining the energetics of the explosions and the nature of the
circumburst environment.

At short wavelengths ($\lambda$ $\leq$ 1.3 mm), VLBI can resolve and image the
emission near the event horizon of nearby supermassive black holes, presenting
an opportunity to study general relativity in a strong gravity regime
(Doeleman et al.,, 2009). Recent VLBI observations at 1.3 mm by the EHT have
captured the first image of the 6.5 $\times$ 109 $M_{\odot}$ supermassive
black hole at the center of the giant elliptical galaxy M87 (Event Horizon
Telescope Collaboration et al., 2019a, ). Previous observations have detected
$\sim$ 5 Schwarzschild radii structure at the base of the relativistic jet
produced at the core of M87 (Doeleman et al.,, 2012; Akiyama et al.,, 2015).
Observations at 1.3 mm have also resolved structures on the scale of a few
Schwarzschild radii in the 4 $\times$ 106 $M_{\odot}$ black hole at the center
of the Milky Way, Sagittarius A∗ (Doeleman et al.,, 2008; Fish et al.,, 2011;
Johnson et al.,, 2015; Lu et al.,, 2018).

###  2.2 Geodesy

Geodesy is the study of the geometric shape, rotation and orientation in space
of Earth and its gravitational field. Since VLBI observations can be used to
determine antenna positions with millimeter accuracy, the technique is also
well suited for geodetic studies. VLBI observations for geodesy are performed
by the International VLBI Service for Geodesy and Astrometry (Schlüter and
Behrend,, 2007), together with the satellite range network and the Global
Navigation Satellite System antenna network (Whitney et al.,, 2013).

VLBI for geodesy is carried out by observing many extragalactic broadband
radio sources that are uniformly distributed across the sky (Schuh and
Behrend,, 2012). Since fast slewing antennas are required to quickly observe
different sources around the sky, geodetic VLBI antennas typically have
smaller diameters. To compensate for the decrease in sensitivity due to the
smaller antenna diameters, higher bandwidths are needed to reach an equivalent
sensitivity. To observe a uniform source distribution, it is often necessary
in geodesy to observe weaker sources, which in turn also requires higher
bandwidths.

Under the current VLBI Global Observing System (VGOS) operation, VLBI
correlation is distributed to multiple computer clusters, where each cluster
is constrained by the number of available playback units and the percentage of
time allocated for
geodesy555https://ivscc.gsfc.nasa.gov/about/org/components/co-list.html. In
contrast, shifting the correlation of geodetic VLBI to the cloud could provide
greater user flexibility by leveraging the massive I/O and compute resources
for a short period of time.

##  3 Cloud correlation architecture

We propose a cloud correlation architecture that leverages the parallel nature
of the correlation workload and the massive parallel I/O provided by
distributed cloud storage. Figure 1 highlights the architectural differences
between a typical VLBI dedicated cluster correlator and the proposed cloud
correlation architecture.

Figure 1: Comparison of standard cluster correlation architecture (left) with
the proposed cloud correlation architecture (right). The proposed data flow is
highlighted with blue arrows, where new data products are shaded green, and
major new components in the development of the cloud pathway are shaded pink.
A typical cluster correlates the data streamed from recorders, where each
cluster is limited by the number of available processor nodes and playback
units as well as the percentage of cluster time allocated to different
experiments. On the other hand, the cloud correlation architecture begins with
the data separated by scans distributed across cloud storage staged using the
Google Transfer Appliances and correlates the entire experiment in parallel
using independent jobs. The ability to process data in parallel is equally
true for a dedicated correlation facility consisting of a large supercluster
that transfers data over network using a parallel file system, resulting in
similar improvements in data processing times. Since not all research groups
have access to a supercluster, however, the cloud correlation framework
provides greater accessibility and flexibility that enables researchers to
leverage massive compute and storage resources for a short period of time. The
cloud also allows international collaborators to monitor and process the data
without the need to physically mount hard drives onto a local computing
cluster.

VLBI observations are normally separated into short integration scans, during
which all antennas in the array simultaneously observe the same target source
and typically record the data onto hard-disk recorders. An example of such
hard-disk recorders are the Mark 6 recorders, which are disk-based data
capture and record systems based on commercial off-the-shell (COTS) hardware
and open-source software with a target sustained data capture rate of greater
than 16 Gbps (Cappallo et al.,, 2013; Whitney et al.,, 2013). Google also
provides similar bulk data recorders based on COTS hardware called Google
Transfer Appliances (GTAs)666https://cloud.google.com/transfer-appliance/. We
consider two possibilities for recording the data at the telescopes. The data
could be directly recorded onto the GTAs. Alternatively, if the current GTA
technology is not robust and reliable enough at present for recording at high
altitudes and extreme temperatures, the data could initially be recorded onto
the existing recorders and then transferred onto the GTAs. Recording onto the
GTAs was not tested in this study, but it is not unreasonable to expect that
the next generation GTAs could function with a similar capacity, speed, and
reliability as the existing VLBI recorders.

Once on the GTAs, the data are shipped to a data center, where they are staged
onto the Google Cloud Storage (GCS). Multiple VMs can now access the raw data
independently, enabling the correlation of many scans to proceed in parallel.
The parallel processing of data is equally true for a dedicated correlation
facility that uses a large parallel file system as a staging area.

The Workflow Manager in Figure 1 represents a suite of potential management
software that can orchestrate and monitor the correlation and post-processing
of data on the cloud. This includes pipeline generation scripts that
initialize a specific correlation or post-processing task in terms of a
sequence of processes and dependencies, and a job scheduler that executes the
pipeline and dynamically allocates VM instances as needed. Although a Workflow
Manager was not used in this study, Kubernetes777https://kubernetes.io is an
example of an open-source, industry-wide system that could be used as a
Workflow Manager platform for cloud correlation. For the benchmarks in this
study, Docker888https://docs.docker.com containers were used to create an ad-
hoc scheme that packaged correlation and post-processing software along with
their dependencies for deployment to the GCP.

##  4 Benchmark correlation setup

We performed benchmarks based on two test data sets derived from a single
recording made on an EHT backend system (Vertatschitsch et al.,, 2015). This
recording was captured at a 2 GHz bandwidth and 2-bits per sample, which when
critically sampled corresponds to an 8 Gbps data rate per polarization. Two
orthogonal polarizations were recorded, and the data processing produces full
polarization correlation products. The two benchmark data sets consist of
copies made from the single recording to create synthetic arrays, one of 10
virtual stations, and the other of 20 virtual stations. Thus, all the virtual
stations in the test data sets have an effective 2 GHz bandwidth, dual-
polarization, 2-bit sampling, and 20-second duration subscans, and all
stations are copies of the same data. To ensure realistic processing and
geometric calculations, each station was given an Earth-centered coordinate
corresponding to a current EHT VLBI site location. For each 20-second duration
subscan, the 10-station data set has a total data volume of 400 GB; the
20-station data set is 800 GB. The 20-second duration subscans were used to
easily accommodate the data size into the memory of the VMs. Since VLBI
observation integrations are typically of longer duration, the 20-second
benchmark results can be scaled linearly to longer scans.

The benchmarks carried out in this work could have been performed on either
synthetically generated pure noise or recordings of noise in a laboratory
setting, but it was convenient to use EHT recordings to verify the cluster
versus cloud agreement for data taken under nominal astronomical observing
conditions. For verification, benchmarks using DiFX-2.5.2 were also performed
on the cluster at the Smithsonian Astrophysical Observatory consisting of
eight 10-core Intel Xeon E5-2640 v4 processors and one 16-core Intel Xeon
E5-2620 v4 processor. The results of running the correlation on the GCP and on
the cluster on identical datasets provided excellent agreement.

The GCP contains various types of VM instances based on the number of vCPUs
and memory. For the n1 series of machine types, a vCPU is implemented as a
single hardware hyper-thread on one of the available CPU platforms. For the
benchmarks, the n1-highmem-16, n1-highmem-32, n1-highmem-64, n1-highmem-96,
and the n1-megamem-96 VMs were used. The specifications and price per hour for
these VMs999https://cloud.google.com/compute/pricing#machinetype are given in
Table 1. VMs with large memory were used to easily accommodate the data volume
of the test data sets in memory. The CPU platforms used were the Intel Xeon E5
(Sandy Bridge), Intel Xeon E5 v2 (Ivy Bridge), Intel Xeon E5 v3 (Haswell), and
Intel Xeon E5 v4 (Broadwell E5) for the highmem VMs up to n1-highmem-64, and
either of the four (Platinum, Gold, Silver, and Bronze) types of the Intel
Xeon Scalable Processor (Skylake) for the n1-highmem-96 and n1-megamem-96
VMs101010https://cloud.google.com/compute/docs/cpu-platforms.

The correlation was performed on the DiFX-2.5.2 software correlator. All
available vCPUs on a VM were used as the DiFX core data processing nodes, even
the ones on which the Datastream or the FXManager processes are run. One
thread per core data processing node was used. The correlation parameters were
set up with a Fast Fourier Transform (FFT) spectral resolution of 0.015625
MHz, an FFT window size of 64 $\mu$s, an FFT size of 262144, an integration
time of 0.512 seconds, and an output spectral resolution of 0.5 MHz.

##  5 Results

The benchmarks were run varying the number of vCPUs and number of stations,
and data transfer rates were also investigated.

Machine type | Number of vCPUs | Memory | Price/hr (USD) | Preemptible price/hr (USD)
---|---|---|---|---
n1-highmem-16 | 16 | 104 GB | $0.94 | $0.20
n1-highmem-32 | 32 | 208 GB | $1.89 | $0.40
n1-highmem-64 | 64 | 416 GB | $3.79 | $0.80
n1-highmem-96 | 96 | 624 GB | $5.68 | $1.20
n1-megamem-96 | 96 | 1433.6 GB | $10.67 | $2.26
Table 1: Specifications and price per hour of the Virtual Machines on the
Google Cloud Platform used in the benchmarks. The quantities are accurate as
of the time of writing.

###  5.1 Number of vCPUs

To estimate the optimal number of vCPUs per VM to use for computation, we used
the n1-highmem-16, n1-highmem-32, n1-highmem-64, and n1-highmem-96 VMs, with
the number of vCPUs varied from 16, 32, 64, and 96, respectively, as well as
the number of stations varied from 2 to 10 stations. The 96-vCPU maximum was
set by the ready availability of VMs up to that size on the GCP.

The results of computational time and cost with number of stations for
different VM vCPUs are shown in Figure 2. The 16 and 32-vCPU benchmarks were
only done up to 5 and 6 stations, respectively, since the computation for even
more stations in those cases would have been too time-consuming. Figure 3
shows the computational time with the number of vCPUs for a 5-station dataset,
suggesting an approximately inverse relationship between computational time
and the number of vCPUs.

To estimate the uncertainty in the measurements, the 96-vCPU, 10-station
measurement was performed twice. The difference in correlation time was less
than 1%. For the 10-station test data set, these results suggest that using 96
vCPUs per VM is the most efficient for data processing.

Figure 2: Computational time in seconds (left) and computational cost in vCPU-
hour (right) with number of stations on the n1-highmem-16, n1-highmem-32,
n1-highmem-64, and n1-highmem-96 VMs. The results suggest that using 96 vCPUs
per VM is the most efficient for data processing. Figure 3: Computational time
in seconds with number of vCPUs for a 5-station dataset on n1-highmem-16,
n1-highmem-32, n1-highmem-64, and n1-highmem-96 vCPU VMs. The result suggests
an approximately inverse relationship between computational time and number of
vCPUs.

###  5.2 Number of stations

To study the scaling of correlation time and thus computational cost with the
number of stations, N, we used the 96-vCPU, n1-megamem-96, 1433.6 GB VM at the
us-west1 data center in Oregon. We varied the number of stations from 2 to 20.
The results are shown in Figure 4. To estimate the uncertainty in the
measurements, the 20-station measurement was performed twice, and the
difference in correlation time between the two measurements was less than
$0.3$%. The measured computational time in seconds as a function of the number
of stations was fit with a quadratic with a R2 of 0.998.

$t\simeq 1003\bigg{(}\frac{\rm N}{10}\bigg{)}^{2}+1060\bigg{(}\frac{\rm
N}{10}\bigg{)}+285.3$ (1)

Figure 4: Correlation time with number of stations of a 20-second, 20-station,
dual-polarization data set with a recorded bandwidth of 2 GHz per polarization
and a data volume of 800 GB. The correlation was performed on the
n1-megamem-96 1433.6 GB VM launched at the us-west1 data center in Oregon. The
linear term is expected due to the FFT per station, and the quadratic term is
expected due to the channel-by-channel cross multiplications per baseline
following the FFT performed by the DiFX-2.5.2 software. The crossover point
between the linear and non-linear terms is N $\sim$ 11.

The processing time is expected to exhibit a linear dependence on station
number because many processing steps in DiFX-2.5.2 are performed once per
station. A quadratic term in station number is also expected because the
channel-by-channel cross-multiplications that follow the FFT are done for each
baseline. The non-linear term begins to dominate at large N with a crossover
point at N $\sim$ 11\. The number of full-polarization correlation products is
4N2, including auto-correlations.

The deviation from purely quadratic dependence at large N could potentially be
due to each Datastream and the FXManager processes taking a significant
fraction of one vCPU’s capacity, effectively limiting the number of total
vCPUs performing the core processing. The set up time, which is the time from
when DiFX-2.5.2 is instantiated until the first bits begin to get correlated,
and the tear down time, which is the time from when the last bits are
correlated to when all the processes are terminated, could be substantial for
the relatively short 20-second subscans used for the benchmarks. This effect
is almost independent of the number of stations. The slowdown at large N could
also be due to the extra memory required because of CPU cache misses described
in (Deller et al.,, 2011). The performance at large N could potentially be
improved by further optimizing the DiFX-2.5.2 correlation parameters.

###  5.3 Data transfer rate

We transferred 10-station VLBI data consisting of a single 20-second subscan
as 20 VLBI Data Interchange Format (VDIF) files each with a size of 20 GB
making a total size of 400 GB, from the GCS to n1-highmem-96, 624 GB VM to
measure the data transfer rate. We then transferred the same data to two
different VMs in parallel to study its effect on the data transfer rate. The
results are shown in Figure 5.

Figure 5: Data transfer rate and transfer time for transferring of VLBI Data
Interchange Format (VDIF) files with a total size of 400 GB from the Google
Cloud Storage to a single n1-highmem-96, 624 GB VM (left) and two of the same
VMs in parallel. The aggregate data transfer rate increases as the number of
files increases until about 12 files (240 GB), after which it saturates at
$\sim$ 6 Gbps. The data transfer rate is not affected when transferring the
same data to two VMs in parallel.

First, the aggregate data transfer rate to a single VM increases as the number
of files increases until about 12 files (240 GB), after which it saturates at
$\sim$ 6 Gbps. Second, the data transfer rate is not affected when
transferring the 10-station data to two VMs in parallel, as shown in the right
panel of Figure 5. These results suggest that it is not unreasonable to expect
that VLBI data separated by scans can be independently transferred from the
GCS to several VMs without a decrease in the data transfer rate, since the
cloud fabric is designed to support independent resource scaling while
simultaneously catering a large collection of users at any given time.

###  5.4 Example correlation configuration

Based on the benchmark results, we describe an example cloud correlation
configuration in Figure 6. In all sections that follow, we consider the
processing of a 10-station dual-polarization VLBI observation consisting of a
total of 200, 5-minute individual scans, constituting a total observation time
of $\sim$ 16.7 hours and a total data size of 1.2 PB recorded at a bandwidth
of 2 GHz per polarization. Each 5-minute scan or 20-second subscan referred to
in Figure 6 already contains data from all 10 stations.

Figure 6: An example cloud correlation configuration of a 10-station, dual-
polarization VLBI data set with 2 GHz bandwidth per polarization, a total size
of 1.2 PB, and a total observation time of $\sim$ 16.7 hours. A total of 200,
5-minute scans are subdivided into 3000, 20-second subscans each with a size
of 400 GB. Each 5-minute scan and 20-second subscan already consists of data
from all 10 stations. Each 20-second subscan is transferred from the Google
Cloud Storage to 3000 individual n1-highmem-96 VMs. Each individual VM
correlates a 20-second subscan in $\sim$ 2400 seconds, which corresponds to a
total correlation time of $\sim$ 4 hours ($\sim$ 25% of total exposure time)
if processed with 6 consecutive instances of using 500 VMs in parallel at a
time. For this example correlation configuration, the total data transfer cost
is $540, and the total computational cost is $2400.

The first step consists of staging the data from the GTAs onto the GCS. Each
5-minute scan from all 10 stations consists of 6 TB of data volume, which is
too large to transfer onto any individual VM.

The second step involves pre-processing. The correlation parameters are
initially set up, after which the data are split into smaller time segments or
multiple frequency channels or transformed into the frequency domain for
faster correlation. For instance, each 5-minute scan could be divided into
smaller time segments that are small enough to load onto individual VMs. The
5-minute scans can be split into 15, 20-second subscans, each with a size of
400 GB, making a total of 3000 subscans for the entire dataset. These subscans
can then be transferred from the GCS to 3000 independent n1-highmem-96, 624 GB
VMs in parallel. The data transfer benchmark results in Figure 5 suggest that
it is not unreasonable to expect that the parallel data transfer rate will not
be affected when transferring to multiple VMs with the assumption that
bottlenecks are not encountered. The transfer of 3000 scans should take $\sim$
530 seconds per VM. Future optimization could allow the transfer of data as it
being processed, as further discussed in Section 6.3.

The third step consists of the data processing. To make the computation more
manageable, the correlation could be performed using six consecutive instances
of simultaneously using 500 subscans or 500 VMs in parallel. From Figure 2,
the expected DiFX-2.5.2 correlation time of a 10-station, 20-second subscan
with a size of 400 GB subscan is $\sim$ 2400 seconds per VM. The total
expected correlation time for all 3000 subscans is $\sim$ 4 hours ($\sim$ 25%
of total exposure time).

The total time for which the data is stored in the cloud will be dominated by
the time required to set up the optimal correlation parameters in the pre-
processing step and the scrutinization of the correlated products afterwards.
It is sometimes necessary to perform multiple trial correlations to determine
the optimal correlation parameters, troubleshoot bugs in the software
correlator, and fix errors in the recorded data. An advantage of cloud
correlation is the faster turn-around time on running trial correlations,
which could make the process of determining the optimal correlation parameters
much more efficient for the correlator operator. The additional time for
setting up the correlation parameters and scrutinizing the correlated data
products will generally be greater for part-time arrays compared to dedicated
full-time arrays like the Very Long Baseline Array and geodetic arrays such as
VGOS.

##  6 Cloud cost estimate

Based on the benchmark results, the costs for processing the example
10-station observation on the cloud can be estimated. The various components
that contribute to the overall cost including the rental of recording units,
shipping, data storage, data transfer, and computation are discussed in this
section.

###  6.1 Recording and shipping

Here, we consider recording directly onto to the GTAs at the telescopes with
the assumption that either the current or the next generation GTAs are
reliable enough to record data at high altitudes. The alternative is to first
record the data onto existing VLBI recorders and then transfer the data to the
GTAs. Recording the 1.2 PB data would require twelve 100 TB GTAs, each with a
usage fee of $300 with 10 free days of on-site usage. The total cost would be
12 $\times$ $300 $=$ $3,600. The two-way shipping cost of the GTAs for a 16
Gbps observation is estimated at $100 per site, or $\sim$ $1,000.

###  6.2 Storage

Data storage dominates the cost of the cloud correlation architecture. On the
GCS, there are four types of storage
classes111111https://cloud.google.com/storage/docs/storage-classes#regional:
regional, multi-regional, nearline, and coldline.

Regional storage consists of data storage in only one specific region with no
redundancy. If storage of data only in a single region with no redundancy is
acceptable for the experiment, which can be expected to be the case for most
VLBI experiments, then regional storage could be used. Regional storage costs
$0.02 per GB per month. The storage cost model on the GCS is not monthly
based, and per-day pricing is also available.

Multi-regional storage is designed for frequently accessed data, is geo-
redundant (stored in at least two different geographic locations), and ensures
maximum data availability, but it is more expensive than regional storage. It
costs $0.026 per GB per month.

Nearline and coldline storage are low-cost storage options for infrequently
accessed data, once a month for nearline and once a year for coldline storage.
The cost is $0.01 per GB per month for nearline storage and $0.007 per GB per
month for coldline storage. However, both nearline and coldline storage have a
retrieval cost of $0.01 per GB and $0.05 per GB, respectively.

Since many radio telescopes are located at remote places such as the South
Pole and on top of mountains, it is not unreasonable to expect delays in
acquiring the GTAs from some stations in the VLBI array. Therefore, the
cheaper nearline or coldline storage could be used for longer term storage
until the GTAs from all stations have been acquired, with the caveat of
retrieval costs. We consider the storage of 1.2 PB in regional storage for one
month for a cost estimate of $25,166. This represents a best-case lower bound,
and assumes no major delays from telescopes at remote locations. Where such
delays exist, costs for GTA rental and cloud storage may be substantially
higher, and nearline or coldline storage needs to be specified at a minimum.

###  6.3 Data transfer to Virtual Machines

In our benchmark operational model, correlation of data does not begin until
raw data from all antennas has been transferred in full to the VMs. The data
transfer of 3000 scans to the VMs in parallel should take $\sim$ 530 seconds
per VM, which would cost $0.18 per VM with a total cost of $0.18 $\times$ 3000
$=$ $540 for all 3000 subscans, considering a preemptible price of $1.20 per
hour for the n1-highmem-96, 624 GB VM. With relatively minor changes to the
underlying software and data model, however, it would be possible instead to
stream data as it is being processed and remove any separate cost for the data
transfer time. The DiFX-2.5.2 software correlator supports a streaming mode in
which the data can be accepted via a network socket. Since data transfer
bandwidth exceeds the rate of correlation, this would not impact correlation
run-time. Because we did not use a streaming data model for the model,
however, we include the costs of transfer time in Table 2.

Virtual Machine type | Record & ship | Storage | Transfer | Computation | Total | Total
---|---|---|---|---|---|---
(16 Gbps) | (16 Gbps) | (16 Gbps) | (16 Gbps) | (16 Gbps) | (64 Gbps)
Preemptible VMs | $4,600 | $25,166 | $540 | $2,400 | $32,706 | $130,824
Non-preemptible VMs | $4,600 | $25,166 | $2,510 | $11,370 | $43,646 | $174,584
Table 2: Cloud cost estimate for the DiFX-2.5.2 software correlation of a
10-station, dual-polarization VLBI dataset recorded at 16 Gbps with 2 GHz
bandwidth per polarization and an observation time of $\sim$ 16.7 hours on the
Google Cloud Platform. The data storage cost is for the storage of 1.2 PB of
data for one month in regional storage in South Carolina (us-east1). A typical
EHT VLBI experiment is recorded at 64 Gbps, which quadruples the 16 Gbps
costs, as shown in the rightmost column and is used as the comparable to
cluster cost in Section 7.

###  6.4 Computation

The DiFX-2.5.2 correlation time of a 400 GB subscan is expected to take $\sim$
2400 seconds, which would cost $0.80 per VM with a total cost of $0.80
$\times$ 3000 $=$ $2,400 for all 3000 subscans, considering a preemptible
price of $1.20 per hour for the n1-highmem-96, 624 GB VM. Here, we do not
consider some up front computational cost of running trial correlations in the
estimate.

Preemptible VMs are affordable, short-lived VMs, which last up to 24 hours.
The disadvantage of using a preemptible VM is that it could be terminated at
any time by the GCP if other higher priority compute instances require that
particular VM. However, since the results suggest that the correlation time is
$\sim$ 4 hours, and the benchmarks were performed using preemptible VMs, we
expect that preemptible VMs would be appropriate for cloud correlation. Non-
preemptible VMs could also be used if necessary. The correlation cost for a
400 GB subscan would be $3.79 per VM with a total cost of $3.79 $\times$ 3000
$=$ $11,370 for all 3000 subscans, considering a non-preemptible price of
$5.68 per hour for the n1-highmem-96, 624 GB VM.

We estimate the total cost for the cloud correlation of the full example
observation in Table 2. For the 16 Gbps recorded benchmark data set, using
preemptible VMs, the total cost is $32,706, whereas using non-preemptible VMs,
the total cost is $43,646. A typical EHT experiment, however, is recorded at
64 Gbps, which quadruples the 16 Gbps costs to $130,824 and $174,584 for
preemptible and non-preemptible VMs, respectively.

##  7 Cluster cost estimate

To provide an economic comparison between the cloud and a physical cluster, we
consider the cost of a 1000-core, CPU-based cluster. Such a cluster is capable
of correlating a 16 Gbps, 10-station VLBI experiment in a time period four
times longer than the corresponding observation cumulative scan time.

###  7.1 Cluster hardware

Equipment for this hypothetical cluster includes: 25 high-performance rack
mount Linux machines, each with a dual-node with each node based on an Intel
Xeon Processor with 20 cores; fast Ethernet switches and cables to tie the
nodes together; 10 Mark6 playback machines each capable of delivering data at
a rate of 16 Gbps; and physical infrastructure such as racking,
uninterruptable power supplies (UPS) and network and power distribution
cables. We have summed quoted costs for this equipment and arrive at an
estimated capital investment of about $500,000. Assuming a five-year lifetime
for the cluster, the cost is amortized linearly at $100,000 per year.

###  7.2 Media

In the cluster case, GTAs are not rented, so 1.2 PB of hard disk drive space
needs to be purchased. We round up slightly given that the Mark6 disks have to
be installed in 8-pack disk modules. At the time of writing, 10 terabyte
enterprise quality Helium disk drives are selling for $330, so this scales to
$42,240. Allowing for 8-pack disk modules for the Mark6, 16 modules at $500
each, the media cost is estimated at $50,000. This is likewise linearly
amortized over five years or $10,000 per year.

Cost type | Cluster | Media | Cluster | Media | IT + | Total | Total
---|---|---|---|---|---|---|---
(16 Gbps) | (16 Gbps) | (64 Gbps) | (64 Gbps) | Elect. | (16 Gbps) | (64 Gbps)
Capital | $500,000 | $50,000 | $500,000 | $200,000 | $700,000 | $1,250,000 | $1,400,000
(no prorate)
Annual | $100,000 | $10,000 | $100,000 | $40,000 | $140,000 | $250,000 | $280,000
(no prorate)
Annual | N/A | N/A | $8,000 | $40,000 | $11,200 | N/A | $59,200
(8% cluster)
Annual (8% cluster | N/A | N/A | $8,000 | $10,000 | $11,200 | N/A | $29,200
25% media)
Table 3: Summary of cluster cost estimates of the correlation of a 10-station,
dual-polarization, 2 GHz bandwidth per polarization, 1.2 PB data set recorded
at 16 Gbps as well as 4.8 GB data set recorded at 64 Gbps with an observation
time of $\sim$ 16.7 hours. Both the capital and annual costs are given
assuming a five-year replacement for the cluster and disks, linear
amortization, disk cost scaling from 16 to 64 Gbps, and prorating of costs
assuming 8% utilization of cluster and 25% utilization of media. For
simplicity, the five-year cumulative costs of electricity and IT services (0.2
FTE per year or one FTE at an assumed loaded $200,000) are considered
“capital” costs.

###  7.3 Power, cooling, and IT support

The power consumption of such a cluster is estimated as 30 kW, and assuming an
energy cost of $0.20 per kWh, the power to run the cluster 24/7/365 amounts to
$50,000 per year. A widely quoted rule-of-thumb for data centers suggests
roughly equal power is needed to cool the cluster. With cooling, the total
annual electric bill is $100,000 per year. We assume 0.2 FTE per year of a
computer systems technician (“IT” support) to assist with the assembly and
maintenance of the cluster, and estimate $40,000 for this service.

###  7.4 Amortized total cost

Under these assumptions, the annualized cost of running the cluster sums to
$250,000 assuming cost of media sufficient for a 16 Gbps experiment. Again, a
typical EHT VLBI experiment is run at a bandwidth of 64 Gbps. Although scaling
from 16 to 64 Gbps in the cloud quadruples the cost as per Table 2, the
cluster would be more heavily utilized for 64 Gbps, but at no greater marginal
cost. This is with the exception of media, for which an annualized cost of
$40,000 to record four times the amount of data must be assumed. This brings
the annual cluster cost to $280,000 for a 64 Gbps experiment. Thus, the
annualized cost of the cluster is more than double the $\sim$ $131,000 cost
for a 64 Gbps cloud correlation.

###  7.5 Prorating the cluster for utilization

The cluster cost estimates assume that it is paid for and run for the full
year. However, it is available to be utilized on other experiments and
computation in general when not running VLBI correlations. Making the
assumption of perfect utilization at other times — which assumes excellent
operations management — it is appropriate to prorate the annual cost of the
cluster for the time it is actually used during the year.

We estimate that a 64 Gbps, 10-station, 4.8 PB experiment would require about
11 days of continuous run time for correlation on the hypothetical 1000-core
cluster. Allowing for setup and disk changes, as the cluster only reads and
correlates 16 Gbps per run, this might reasonably translate into about a month
of wall clock time or about 8% of the year. If cluster costs are prorated for
8% of the year — though not prorating the costs of media — the cluster is then
the more economical choice at $59,200.

###  7.6 Prorating the media

We also consider the possibility of prorating media, which can be recycled
through multiple observations, as is common in VLBI. However, the 8% factor
used to prorate the cluster for a month of use is not appropriate in the case
of media. The Mark6 disk packs, for instance, are shipped to the sites well
ahead of the observation and are conditioned with read-write cycles there,
which is a time-consuming process. After the observation, they are packed and
shipped back to the correlator for processing. The data continues to reside on
the Mark6 media modules until the correlation job is completed, for which we
have assumed one month of clock time.

We estimate two additional months on each side of the observation, to pre-
check, pack and ship the disks to the site, including clearing customs; to
condition them and then run the observation; and then pack, ship, unpack and
stage back at the correlator. In principle, the cluster can be used for other
processing jobs during the observation and shipping time, but the disks are
not available to recycle into other observations. With the one month of
processing in addition to the two months for shipping and staging, we prorate
disks to three months, or 25% of annual. Of course, other observation
opportunities are available at exactly the right times for the media to be
reused four times a year. Under this assumption, the cluster cost is lower
still than the cloud for a 64 Gbps observation.

###  7.7 Cost comparison summary

Table 3 provides a summary of the three types of scaling relations used for
the cluster cost estimates in the prior discussion:

  * •

amortization: the capital investment considered over an assumed five-year
lifetime of the equipment.

  * •

bandwidth scaling: our benchmarks are for a 16 Gbps recording while a typical
EHT wideband experiment is recorded at 64 Gbps.

  * •

prorating: the sharing of capital investment to the benefit of multiple
projects.

| Cloud | Cluster | Cluster | Cluster
---|---|---|---|---
| (no prorate) | (8% cluster prorate) | (8% cluster, 25% media prorate)
Total (64 Gbps) | $130,824 | $280,000 | $59,200 | $29,200
Table 4: Cloud and cluster cost comparison of the correlation of a 10-station,
dual-polarization, 4.8 PB dataset recorded at 64 Gbps with 2 GHz bandwidth per
polarization with an observation time of $\sim$ 16.7 hours. Three cases are
considered for the cluster cost: no prorating, 8% prorating of the cluster,
and 8% prorating of the cluster as well as 25% prorating of the media. The
prorated cluster costs are lower than the cloud cost. In practice, however,
prorating the cluster is almost certainly overly optimistic, especially in the
case of media.

Table 4 provides a cloud and cluster cost comparison for the correlation of
4.8 GB data set recorded at 64 Gbps, where the cluster cost is considered
under three cases: no prorating, 8% cluster prorating, and 8% cluster and 25%
media prorating. We note that the primary benefit of the cloud is not the
cost, but rather the ability to deploy much larger computational resources
than practically possible on a cluster to run the computation much faster and
thereby shorten the time-to-science. This in turn implicitly improves the
effective utilization of the entire astronomical instrument, not merely the
correlator. Deploying more cores for a shorter time in the cloud does not
significantly affect the cost of the computation. To the contrary, running the
computation faster with a greater number of cores reduces the costs of
regional storage, which are the dominant cloud costs. However, it is not
unreasonable to expect the storage costs on the cloud to decrease with time.

Prorating the cluster and the disk packs is only valid if the rest of the year
the equipment is effectively utilized on other compute jobs (cluster) and
observations (media). In practice, the prorating is almost certainly overly
optimistic, especially in the case of media, where observations have to become
available on precisely a three-month cadence for the 25% prorating to be
achieved. The prorating calculation is useful as a guide to what is possible
in principle. Further, the cluster estimates presented in this subsection are
intended only as a rough basis for comparison, with a proper cost comparison
being highly assumption dependent, rather complex to put together, and
dependent on cost figures which vary greatly over time.

Viewing utilization actually achieved as an economic parameter, there exists a
crossover point between the cloud and cluster cost. For low duty cycle
wideband observations where substantial cloud resources are strategically
deployed, the dominant storage costs for cloud are reduced. We note that when
not utilized, cloud resources are simply released back to the provider,
whereas for the cluster, efficient utilization is a significant logistical
problem for the operating institution.

In practice, corporate cloud service providers may provide services at a
significant discount, or even gratis, to scientific research groups and non-
profits. This practical factor could dramatically affect the comparative
retail costs assumed here, substantially affecting the cost crossover point
between the cloud and the cluster.

##  8 Summary

The trend for VLBI observations to record at greater bandwidths has led to
increased data volumes. This paper demonstrates the viability of porting wide
bandwidth VLBI correlation processing to cloud platforms as a promising
alternative to cluster correlation.

Benchmarks run on 10-station and 20-station synthetic data sets using the
DiFX-2.5.2 software correlator indicate that full wideband VLBI observations
with data volumes exceeding a petabyte can be efficiently processed in elapsed
times that are much shorter than the time spent during on-sky observing.
Similar to some VLBI arrays that transfer data over network and process the
correlation in a large computer cluster, the cloud allows the processing of
large data sets in parallel. Clusters are often limited by the number of
installed compute nodes and playback units as well as the percentage of time
allocated for processing different science experiments. Since not every
research group has access to a large computer cluster, the advantage of the
cloud is its greater accessibility and flexibility that allows researchers to
leverage massive compute and storage resources for a short period of time.

Cloud-based computation and data storage enables distributed operation of
correlation, allowing international collaborators to efficiently monitor and
process VLBI data without the need for physically mounting hard disk media
onto localized computing clusters.

The cloud architecture also ensures built-in upgrades of the rented recording
appliances, storage, and compute hardware that can be expected to keep up with
industry trends and takes advantage of new commodity equipment as vendors
naturally adopt next-generation hardware and software in their cloud
platforms. An additional benefit is that all maintenance and operational costs
of computing resources are outsourced to the cloud vendor. This study found
that the storage of data dominates the overall cost of cloud correlation.

An approximate comparison of cloud to cluster costs is provided. Current cloud
service and equipment pricing data is used to compile cost estimates which
allow an approximate economic comparison between cloud and cluster processing.
Cluster cost is highly dependent on efficient utilization of the physical
cluster and associated media. Practically speaking, utilization of cluster and
dedicated media inventory is a challenge to optimize, whereas the cloud
resources and transfer appliances have costs strictly proportional to time in
use.

We intend to continue this research, including looking into increased data
rates on the GTAs, and upgrading them to perform reliably at high elevation
sites. A test observation is envisaged as an acid test of the operational
framework. This further research aims to establish whether correlation on the
cloud can provide economical extensible path for wideband VLBI arrays to
achieve breakthrough science possible through the efficient analysis of big-
data.

##  9 Acknowledgments

This work was supported by grants from the National Science Foundation
(AST-1440254, AST-1743747 & AST-1828513) and the Gordon and Betty Moore
Foundation (GBMF-5278). This work was supported in part by the Black Hole
Initiative at Harvard University, which is funded by a grant from the John
Templeton Foundation. We thank Geoffrey B. Crew, Michael Titus, and Victor
Pankratius for assistance with the CorrelX and DiFX-2.5.2 software packages,
and Robert W. Wilson for helpful discussions. We gratefully acknowledge
support from the Google Cloud Platform to carry out benchmarks in this work.
Early developments and benchmarks were also done on the CyVerse Atmosphere
Cloud (supported by the National Science Foundation under Award Numbers
DBI-0735191 and DBI-1265383) and the Jetstream Cloud (Stewart et al.,, 2015;
Towns et al.,, 2014). The VLBI software correlation was performed using the
DiFX-2.5.2 software correlator (Deller et al.,, 2011).

## References

  * Abbott et al., (2017) Abbott, B. P. et al. (2017).  A gravitational-wave standard siren measurement of the Hubble constant.  Nature, 551(7678):85–88.
  * Akiyama et al., (2015) Akiyama, K. et al. (2015).  230 GHz VLBI observations of M87: event-horizon-scale structure during an enhanced very-high-energy $\gamma$-ray state in 2012.  ApJ, 807(2):150.
  * Bartel et al., (2000) Bartel, N. et al. (2000).  The Changing Morphology and Increasing Deceleration of Supernova 1993J in M81.  Science, 287:112–116.
  * Bietenholz et al., (2014) Bietenholz, M. F., De Colle, F., Granot, J., Bartel, N., and Soderberg, A. M. (2014).  Radio limits on off-axis GRB afterglows and VLBI observations of SN 2003gk.  MNRAS, 440(1):821–832.
  * Cappallo et al., (2013) Cappallo, R., Ruszczyk, C., and Whitney, A. (2013).  Mark6: Design and Status.  Haystack Mark 6 Memo Series, Memo 005.
  * Clark et al., (1967) Clark, B. G., Cohen, M. H., and Jauncey, D. L. (1967).  Angular Size of 3c 273B.  ApJL, 149:L151.
  * Deller et al., (2011) Deller, A. T., Brisken, W. F., Phillips, C. J., Morgan, J., Alef, W., Cappallo, R., Middelberg, E., Romney, J., Rottmann, H., Tingay, S. J., and Wayth, R. (2011).  DiFX-2: A more flexible, efficient, robust, and powerful software correlator.  PASP, 123(901):275–287.
  * Deller et al., (2007) Deller, A. T., Tingay, S. J., Bailes, M., and West, C. (2007).  DiFX: A software correlator for very long baseline interferometry using multiprocessor computing environments.  PASP, 119(853):318–336.
  * Doeleman et al., (2008) Doeleman, S. et al. (2008).  Event-horizon-scale structure in the supermassive black hole candidate at the Galactic Centre.  Nature, 455:78.
  * Doeleman et al., (2009) Doeleman, S. et al. (2009).  Imaging an Event Horizon: submm-VLBI of a Super Massive Black Hole.  In Astro2010: The Astronomy and Astrophysics Decadal Survey, volume 2010, page 68.
  * Doeleman et al., (2012) Doeleman, S. S. et al. (2012).  Jet-Launching Structure Resolved Near the Supermassive Black Hole in M87.  Science, 338(6105):355.
  * Doeleman et al., (2004) Doeleman, S. S., Hough, D. H., Kulkarni, S. R., Lonsdale, C. J., Marscher, A. P., O’Dea, C. P., Taylor, G. B., Wilner, D. J., and Wrobel, J. M. (2004).  Mapping the Future of VLBI Science in the US.
  * (13) Event Horizon Telescope Collaboration et al. (2019a).  First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole.  ApJ, 875(1):L1.
  * (14) Event Horizon Telescope Collaboration et al. (2019b).  First M87 Event Horizon Telescope Results. II. Array and Instrumentation.  ApJ, 875(1):L2.
  * (15) Event Horizon Telescope Collaboration et al. (2019c).  First M87 Event Horizon Telescope Results. IV. Imaging the Central Supermassive Black Hole.  ApJ, 875(1):L4.
  * Fish et al., (2013) Fish, V. et al. (2013).  High-Angular-Resolution and High-Sensitivity Science Enabled by Beamformed ALMA.  arXiv e-prints, page arXiv:1309.3519.
  * Fish et al., (2011) Fish, V. L. et al. (2011).  1.3 mm Wavelength VLBI of Sagittarius A*: Detection of Time-variable Emission on Event Horizon Scales.  ApJ, 727(2):L36.
  * Johnson et al., (2015) Johnson, M. D. et al. (2015).  Resolved Magnetic-Field Structure and Variability Near the Event Horizon of Sagittarius A*.  Science, 350(6265):1242–1245.
  * Koopmans and Treu, (2002) Koopmans, L. V. E. and Treu, T. (2002).  The Stellar Velocity Dispersion of the Lens Galaxy in MG 2016+112 at z=1.004.  ApJL, 568:L5–L8.
  * Lu et al., (2018) Lu, R. S. et al. (2018).  Detection of Intrinsic Source Structure at $\sim$ 3 Schwarzschild Radii with Millimeter-VLBI Observations of SAGITTARIUS A∗.  ApJL, 859:60.
  * Mao et al., (2001) Mao, S., Witt, H. J., and Koopmans, L. V. E. (2001).  The influence of central black holes on gravitational lenses.  MNRAS, 323:301–307.
  * Matthews et al., (2017) Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. (2017).  The ALMA phasing system: A beamforming capability for ultra-high-resolution science at (sub)millimeter wavelengths.  PASP, 130(983):015002.
  * Miller-Jones et al., (2009) Miller-Jones, J. C. A., Jonker, P. G., Dhawan, V., Brisken, W., Rupen, M. P., Nelemans, G., and Gallo, E. (2009).  The First Accurate Parallax Distance to a Black Hole.  ApJ, 706(2):L230–L234.
  * Moran et al., (1967) Moran, J. M., Crowther, P. P., Burke, B. F., Barrett, A. H., Rogers, A. E. E., Ball, J. A., Carter, J. C., and Bare, C. C. (1967).  Spectral line interferometry with independent time standards at stations separated by 845 kilometers.  Science, 157(3789):676–677.
  * Patel et al., (2014) Patel, N. A., Wilson, R. W., Primiani, R., Weintroub, J., Test, J., and Young, K. (2014).  Characterizing the Performance of a High-Speed ADC for the Sma Digital Backend.  Journal of Astronomical Instrumentation, 3:1.
  * Reid et al., (2009) Reid, M. J. et al. (2009).  Trigonometric Parallaxes of Massive Star Forming Regions: VI. Galactic Structure, Fundamental Parameters and Non-Circular Motions.  ApJ, 700(1):137–148.
  * Reid and Honma, (2014) Reid, M. J. and Honma, M. (2014).  Microarcsecond Radio Astrometry.  Annu. Rev. Astronomy and Astrophysics, 52:339–372.
  * Reid et al., (2011) Reid, M. J., McClintock, J. E., Narayan, R., Gou, L., Remillard, R. A., and Orosz, J. A. (2011).  The Trigonometric Parallax of Cygnus X-1.  ApJ, 742(2):83.
  * Schlüter and Behrend, (2007) Schlüter, W. and Behrend, D. (2007).  The International VLBI Service for Geodesy and Astrometry (IVS): current capabilities and future prospects.  J. Geodesy, 81(6):379–387.
  * Schuh and Behrend, (2012) Schuh, H. and Behrend, D. (2012).  VLBI: A fascinating technique for geodesy and astrometry.  Journal of Geodynamics, 61:68 – 80.
  * Stewart et al., (2015) Stewart, C. A., Cockerill, T. M., Foster, I., Hancock, D., Merchant, N., Skidmore, E., Stanzione, D., Taylor, J., Tuecke, S., Turner, G., Vaughn, M., and Gaffney, N. I. (2015).  Jetstream: A self-provisioned, scalable science and engineering cloud environment.  In Proceedings of the 2015 XSEDE Conference: Scientific Advancements Enabled by Enhanced Cyberinfrastructure, XSEDE ’15, pages 29:1–29:8, New York, NY, USA. ACM.
  * Taylor et al., (1998) Taylor, G. B., Frail, D. A., Kulkarni, S. R., Shepherd, D. S., Feroci, M., and Frontera, F. (1998).  The Discovery of the Radio Afterglow from the Optically Dim Gamma-Ray Burst of 1998 March 29.  ApJ, 502(2):L115–L118.
  * Thompson et al., (2017) Thompson, A. R., Moran, J. M., and Swenson, Jr., G. W. (2017).  Interferometry and Synthesis in Radio Astronomy, 3rd Edition.
  * Towns et al., (2014) Towns, J. et al. (2014).  XSEDE: Accelerating Scientific Discovery.  Computing in Science and Engineering, 16(5):62–74.
  * van Langevelde, (2006) van Langevelde, H. J. (2006).  Data processing software for radio astronomy.  In Proceedings of the 8th European VLBI Network Symposium, page 55.
  * Vertatschitsch et al., (2015) Vertatschitsch, L., Primiani, R., Young, A., Weintroub, J., Crew, G. B., McWhirter, S. R., Beaudoin, C., Doeleman, S., and Blackburn, L. (2015).  R2DBE: A Wideband Digital Backend for the Event Horizon Telescope.  PASP, 127(958):1226–1239.
  * Walker et al., (2007) Walker, C., Romney, J., Brisken, W., and Durand, S. (2007).  A proposal to increase the sensitivity of the VLBA.  NRAO VLBA Memo 15.
  * Weston et al., (2017) Weston, S., Gulyaev, S., Kim, B., Litchfield, A., Hall, D., Curtis, J., Yun, T. C., Ruthven, A., Davies, G., Lagos, B., and Christie, D. (2017).  Towards Cloud Correlation of VLBI Data.  In Haas, R. and Elgered, G., editors, 23rd European VLBI Group for Geodesy and Astrometry Working Meeting, pages 51–55.
  * Whitney et al., (2013) Whitney, A. R. et al. (2013).  Demonstration of a 16 Gbps Station-1 Broadband-RF VLBI System.  PASP, 125:196.
  * Winn et al., (2004) Winn, J. N., Rusin, D., and Kochanek, C. S. (2004).  The Central image of a gravitationally lensed quasar.  Nature, 427:613.
  * Zackrisson and Riehm, (2010) Zackrisson, E. and Riehm, T. (2010).  Gravitational Lensing as a Probe of Cold Dark Matter Subhalos.  Advances in Astronomy, 2010:478910.


# LHAASO discovery of highest-energy photons towards new physics

Chengyi Li  Bo-Qiang Ma  School of Physics and State Key Laboratory of Nuclear
Physics and Technology, Peking University, Beijing 100871, China  Center for
High Energy Physics, Peking University, Beijing 100871, China  Collaborative
Innovation Center of Quantum Matter, Beijing, China

###### Abstract

Ultrahigh-energy photons up to 1.4 peta-electronvolts have been observed by
new cosmic-ray telescope in China [1] – a hint that Lorentz invariance might
break down at the Planck-scale level [2].

###### keywords:

ultrahigh-energy gamma-ray, Lorentz invariance violation, quantum gravity,
string theory model, space-time foam

††journal: Science Bulletin, published as Sci. Bull. 66 (2021) 2254-2256,
https://doi.org/10.1016/j.scib.2021.07.030

As one of the most sensitive gamma-ray detector arrays currently operating in
the ultrahigh-energy (UHE) band, the Large High Altitude Air Shower
Observatory (LHAASO) recently reported large numbers of gamma-ray photons with
energies larger than $100$ TeV from twelve cosmic accelerators within the
Milky Way [1], including the most high-energetic $\gamma$-ray detected at
about $1.4$ PeV. The observation of such high-energy lights from space can be
used not only to investigate the long-lasting puzzle of the origin of UHE
cosmic rays and their acceleration mechanism in the extreme regions of the
Universe, but can provide also the unique and crucial opportunity in testing
fundamental concepts of physics, such as one basic symmetry underlying
Einstein’s relativity – Lorentz invariance [2].

Set more than 4.4 km above sea level on the Haizi Mountain in Sichuan, China,
the LHAASO is a new-generation mountain observatory [3] with unprecedented
capability to spot UHE particles that rain down on Earth from cosmos. Instead
of tracking cosmic rays that carry electric charges and hence are always bent
by cosmological magnetic fields, LHAASO aims to detect high-energy
$\gamma$-radiation photons, as they travel through the Universe in straighter
lines and are relatively easier to be traced. The observatory consists of
three interconnected arrays of detectors, the Kilometer Square Array (KM2A),
the Water Cherenkov Detector Array (WCDA), and the Wide Field-of-view (air)
Cherenkov Telescope Array (WFCTA), sprawling over a total area of about 1.36
km2. Amongst these three components of the LHAASO instrument, KM2A was
designed to measure ultrahigh-energy cosmic gamma-ray photons of energies
above $10$ TeV. The construction of the observatory started in November 2017.
It then began operations two years later, when LHAASO was not even half-
complete. Specifically, the partly completed KM2A detection array has been
operated stably to continuously collect the data since the end of 2019.
However the LHAASO-KM2A instrument is still under construction and can be
fully operational by July 2021.

Recently, the LHAASO collaboration [1] announced the detection of more than
530 UHE photons coming from twelve astrophysical gamma-ray sources gathering
around the Galactic disc (Figure 1), using the KM2A half-array, with
statistical significance of the photon signals from each source greater than
$7\sigma$ over the surrounding background. These gamma-ray photons possess
energies above 100 TeV up to the highest exceeding $1$ PeV, about a few dozen
to a hundred times the particle energies that the Large Hadron Collider (LHC)
can produce, by noting that LHC reaches energies of $6.5$ TeV per particle in
its second operational run. Amongst these LHAASO events, the most energetic
photon with energy at $1.4$ PeV was found from the source LHAASO J2032+4102,
detected with a $10.5\sigma$ significance, and located in a very active star-
forming region in the constellation of Cygnus. This record-breaking event
turns out to be the highest energy photon ever observed by human being. These
new observational results are of broad attention to the scientific society and
are extremely encouraging for those studies aiming to search for the origins
of UHE cosmic rays. Despite having several potential particle factories (i.e.,
PeVatrons) within our Galaxy, including pulsar wind nebulae, remnants of
supernova explosions, and star-forming regions [1], the cosmological objects
responsible for these UHE gamma-rays have not yet been firmly localized and
identified, leaving open the origin of these extreme cosmic accelerators.

Figure 1: (Color online) Ultrahigh-energy gamma-ray photons up to 1.4 PeV from
twelve stable gamma-ray sources were detected by China’s high-altitude
observatory, LHAASO. Copyright $\copyright$ 2021, LHAASO Collaboration.

The first discovery of cosmic gamma-rays at such high energies with LHAASO is
of great significance. Since cosmic rays are high-energy particles from space,
they may encode important information of the source objects. Precise
measurements of these high-energy photons may provide clues to the
astrophysical origins of cosmic rays and the processes that accelerate
particles to such great speeds from the surrounding extreme regions and can
give valuable hints on the evolution of celestial bodies and the early history
of the Universe. On the other hand, such high-energetic photons in the PeV
range may also enable us to test fundamental principles of physics and probe
new physics that goes beyond the standard model (SM) of particle physics and
general relativity (GR) of gravitation, such as ideas of violation of Lorentz
invariance (LV) [2, 4].

In Einstein’s theory of special relativity, the speed of light and the laws of
physics appear to be the same in all frames. The underlying symmetry, Lorentz
invariance, is also the one that is mandated in the SM of particle physics.
However, some extensions of the SM, especially those quantum gravity (QG)
scenarios that try to reconcile quantum mechanics with GR, motivate the
breaking down of exact Lorentz symmetry in the proximity of the Planckian
regime. Such a symmetry violation can lead to observable consequences that
could serve as QG imprints at lower energies [4]. Due to LV impacts, for
instance, photons would no longer propagate through space at the constant
light speed, but at a velocity that depends on the energy of that photon. In
addition, photons with energies above a certain threshold would decay away.
Such photon self-decay process is generally forbidden in the SM, since it
breaks energy-momentum conservation. The ordinary reaction that two photons
can interact and produce an electron-positron pair with an existing threshold
could also possess an anomaly [5] relative to special relativity if Lorentz
invariance does not hold. These new-physics effects are so miniscule that it
is hard to probe them from most terrestrial experiments, even at the LHC,
whilst the high-energy astrophysical observation, such as those carried out by
the LHAASO instrument, provides us an unprecedented laboratory to explore
Lorentz symmetry violations.

Using the latest LHAASO data, we point out [2] that the first detection of PeV
gamma-rays may offer a very sensitive probe to test new physics of Lorentz
violation. The detailed measurement of these UHE events by LHAASO allows us to
constrain the photon decay motivated by superluminal type of Lorentz violation
to a higher level [2, 6, 7], leading towards strong limits to certain Lorentz-
violating frameworks. Quite remarkably, our result sets the tightest bound on
the first-order Lorentz-violation energy scale ever reported, in excess of
$\sim 2.7\times 10^{24}$ GeV, over about 200,000 times the Planck scale. This
is deduced from the fact that LHAASO should not observe the 1.4 PeV photon
that would undergo a rapid self-decay if superluminal Lorentz violation is
allowed to take place. This result is very similar to that derived recently in
Refs. [6, 7], using the same scenario.

Moreover, we indicate that a different type, i.e., subluminal type, of Lorentz
invariance violation may be viable. Actually, a subluminal light speed
variation at the scale of $\sim 3.6\times 10^{17}$ GeV has been proposed
previously from analyses on gamma-ray burst (GRB) data of multi-GeV photons
[8, 9, 10]. We further suggest that such type of subluminal deformation of
Lorentz invariance is permitted to exist at the Planck-scale level in order to
account for the latest LHAASO observation of PeV photons. This is based upon
the observation that subluminal Lorentz violation could drastically alter the
conventional kinematics of photon-photon annihilation processes [2], as
discussed also in Ref. [5]. Thus these PeV gamma-rays would have chance to
reach the Earth without attenuated by low-energy background lights diffusing
in the sky (e.g., cosmic background radiations) during their cosmological
propagation. We propose [2, 5] that any identification of such high-energy
cosmic photons from extragalactic sources can be considered as support for
Lorentz-violating properties of photons of the kind listed above.

On the theoretical side, we notice [2] also a more fascinating implication for
the LHAASO discovery, that is, certain models could account for all observed
phenomenology on Lorentz violation within the photon sector, available to
date, including the latest LHAASO result. In fact, it has been noted [11, 12]
that specific QG theories, especially the D-brane/string-inspired model for
space-time foam [13, 14, 15] would be capable of explaining the light speed
variation from GRB photons [8, 9, 10], in agreement with many other
astrophysical data currently available. We further point out that this stringy
QG framework may be also compatible with the data on the LHAASO UHE gamma-ray,
since this type of models of (super)string theories naturally prevents the
self decay for photons at high energies, but may predict an astrophysical
signature of photon-photon annihilation threshold anomalies which, as explain
above, exactly is what the new LHAASO experiment might indicate.

In the future, more observations with LHAASO could cast more stringent
constraints to certain types (e.g., superluminal type) of Lorentz invariance
violation, but might also serve as supports for some Lorentz-violating QG
schemes such as the string-inspired quantum space-time foam. These new
incoming observations can thus offer hopefully valuable perspectives on the
ultimate fate of the Planck-scale Lorentz symmetry.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Acknowledgments

This work is supported by National Natural Science Foundation of China (Grant
No. 12075003).

## References

## References

  * [1] Z. Cao, F. A. Aharonian, Q. An, et al. [LHAASO Collaboration], Ultrahigh-energy photons up to 1.4 petaelectronvolts from 12 $\gamma$-ray Galactic sources, Nature 594 (2021) 33, https://doi.org/10.1038/s41586-021-03498-z.
  * [2] C. Li, B.-Q. Ma, Ultrahigh-energy photons from LHAASO as probes of Lorentz symmetry violations, Phys. Rev. D 104 (2021) 063012, https://doi.org/10.1103/PhysRevD.104.063012, arXiv:2105.07967.
  * [3] Z. Cao, et al. [LHAASO Collaboration], Introduction to Large High Altitude Air Shower Observatory (LHAASO), Chin. Astron. Astrophys. 43 (2019) 457, https://doi.org/10.1016/j.chinastron.2019.11.001.
  * [4] L. Shao, B.-Q. Ma, Lorentz violation effects on astrophysical propagation of very high energy photons, Mod. Phys. Lett. A 25 (2010) 3251, https://doi.org/10.1142/S0217732310034572, arXiv:1007.2269.
  * [5] H. Li, B.-Q. Ma, Threshold anomalies of ultra-high energy cosmic photons due to Lorentz invariance violation, J. High Energy Astrop. 32 (2021) 1, https://doi.org/10.1016/j.jheap.2021.07.001, arXiv:2105.06647.
  * [6] Z. Cao, F. A. Aharonian, Q. An, et al. [LHAASO Collaboration], Exploring Lorentz invariance violation from ultra-high-energy gamma rays observed by LHAASO, arXiv:2106.12350.
  * [7] L. Chen, Z. Xiong, C. Li, S. Chen, H. He, Strong constraints on Lorentz violation using new $\gamma$-ray observations around PeV, Chin. Phys. C 45 (2021) 095104, https://doi.org/10.1088/1674-1137/ac1166, arXiv:2105.07927.
  * [8] H. Xu, B.-Q. Ma, Light speed variation from gamma-ray bursts, Astropart. Phys. 82 (2016) 72, https://doi.org/10.1016/j.astropartphys.2016.05.008, arXiv:1607.03203.
  * [9] H. Xu, B.-Q. Ma, Light speed variation from gamma ray burst GRB 160509A, Phys. Lett. B 760 (2016) 602, https://doi.org/10.1016/j.physletb.2016.07.044, arXiv:1607.08043.
  * [10] H. Xu, B.-Q. Ma, Regularity of high energy photon events from gamma ray bursts, J. Cosmol. Astropart. Phys. 01 (2018) 050, https://doi.org/10.1088/1475-7516/2018/01/050, arXiv:1801.08084.
  * [11] C. Li, B.-Q. Ma, Light speed variation in a string theory model for space-time foam, Phys. Lett. B 819 (2021) 136443, https://doi.org/10.1016/j.physletb.2021.136443, arXiv:2105.06151.
  * [12] C. Li, B.-Q. Ma, Light speed variation with brane/string-inspired space-time foam, Results Phys. 26 (2021) 140380, https://doi.org/10.1016/j.rinp.2021.104380, arXiv:2109.07096.
  * [13] J. R. Ellis, N. E. Mavromatos, M. Westmuckett, Supersymmetric D-brane model of space-time foam, Phys. Rev. D 70 (2004) 044036, https://doi.org/10.1103/PhysRevD.70.044036, arXiv:gr-qc/0405066.
  * [14] J. R. Ellis, N. E. Mavromatos, D. V. Nanopoulos, Derivation of a vacuum refractive index in a stringy space-time foam model, Phys. Lett. B 665 (2008) 412, https://doi.org/10.1016/j.physletb.2008.06.029, arXiv:0804.3566.
  * [15] T. Li, N. E. Mavromatos, D. V. Nanopoulos, D. Xie, Time delays of strings in D-particle backgrounds and vacuum refractive indices, Phys. Lett. B 679 (2009) 407, https://doi.org/10.1016/j.physletb.2009.07.062, arXiv:0903.1303.