MNRAS
000
, 116 (2021)
Preprint 24 February 2022
Compiled using MNRAS L
A
T
E
X style le v3.0
MeqSilhouette v2: Spectrally-resolved polarimetric synthetic data
generation for the Event Horizon Telescope
Iniyan Natarajan
1
2
3
★
, Roger Deane
1
4
, Iván Martí-Vidal
5
6
, Freek Roelofs
7
8
9
,
Michael Janssen
10
9
, Maciek Wielgus
8
7
10
, Lindy Blackburn
8
7
, Tariq Blecher
3
2
,
Simon Perkins
2
, Oleg Smirnov
3
2
, Jordy Davelaar
11
12
9
, Monika Moscibrodzka
9
,
Andrew Chael
13
, Katherine L. Bouman
14
, Jae-Young Kim
10
15
, Gianni Bernardi
16
3
2
,
Ilse van Bemmel
17
, Heino Falcke
9
, Feryal Özel
18
, Dimitrios Psaltis
18
1
Wits Centre for Astrophysics, School of Physics, University of the Witwatersrand, Private Bag 3, 2050, Johannesburg, South Africa
2
South African Radio Astronomy Observatory, Observatory 7925, Cape Town, South Africa
3
Centre for Radio Astronomy Techniques and Technologies, Department of Physics and Electronics, Rhodes University, Makhanda 6140,
South Africa
4
Department of Physics, University of Pretoria, Hat eld, Pretoria, 0028, South Africa
5
Departament d'Astronomia i Astrofísica, Universitat de València, C. Dr. Moliner 50, E-46100 Burjassot, València, Spain
6
Observatori Astronòmic, Universitat de València, C. Catedràtico José Beltrán 2, E-46980 Paterna, València, Spain
7
Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA
8
Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA
9
Department of Astrophysics, Institute for Mathematics, Astrophysics and Particle Physics (IMAPP), Radboud University, PO Box 9010,
6500 GL N¼megen, The Netherlands
10
Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany
11
Department of Astronomy and Columbia Astrophysics Laboratory, Columbia University, 550 W 120th Street, New York, NY 10027, USA
12
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
13
Princeton Center for Theoretical Science, Jadwin Hall, Princeton University, Princeton, NJ 08544, USA
14
California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA
15
Korea Astronomy and Space Science Institute, Daedeok-daero 776, Yuseong-gu, Daejeon 34055, Republic of Korea
16
INAF-Istituto di Radioastronomia, via Gobetti 101, 40129 Bologna, Italy
17
Joint Institute for VLBI ERIC, Oude Hoogeveensed¼k 4, 7991 PD Dwingeloo, Netherlands
18
Department of Astronomy and Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present MeqSilhouette v2.0 (MeqSv2), a fully polarimetric, time-and frequency-resolved
synthetic data generation software for simulating millimetre (mm) wavelength very long base-
line interferometry (VLBI) observations with heterogeneous arrays. Synthetic data are a critical
component in understanding real observations, testing calibration and imaging algorithms, and
predicting performance metrics of existing or proposed sites. MeqSv2 applies physics-based
instrumental and atmospheric signal corruptions constrained by empirically-derived site and
station parameters to the data. The new version is capable of applying instrumental polarization
eects and various other spectrally-resolved eects using the Radio Interferometry Measure-
ment Equation (RIME) formalism and produces synthetic data compatible with calibration
pipelines designed to process real data. We demonstrate the various corruption capabilities
of MeqSv2 using dierent arrays, with a focus on the eect of complex bandpass gains on
closure quantities for the EHT at 230 GHz. We validate the frequency-dependent polarization
leakage implementation by performing polarization self-calibration of synthetic EHT data us-
ing PolSolve. We also note the potential applications for cm-wavelength VLBI array analysis
and design and future directions.
Key words:
techniques: interferometric techniques: high angular resolution
★
E-mail: iniyan.natarajan@wits.ac.za
©
2021 The Authors
arXiv:2202.11478v1 [astro-ph.IM] 23 Feb 2022
2
Natarajan et al.
1 INTRODUCTION
Very long baseline interferometry (VLBI) enables the highest angu-
lar resolution in astronomy, on the order of milli-arcseconds (mas)
to micro-arcseconds (
휇
as), by operating radio antennas separated
by thousands of kilometres synchronously using atomic clocks. The
Event Horizon Telescope (EHT, Event Horizon Telescope Collabo-
ration et al. 2019b) is a global mm-VLBI array whose principal goal
is to spatially resolve the supermassive black holes at the cores of
the Milky Way galaxy (Sgr A
) and M87, and image their
shadows
,
the depression in observed intensity inside the apparent bound-
ary of the black hole (e.g. Falcke et al. 2000; Dexter et al. 2012;
Psaltis et al. 2015; Mo±cibrodzka et al. 2016), together with a bright
crescent-shaped emission ring (e.g. Bromley et al. 2001; Broder-
ick & Loeb 2009; Kamruddin & Dexter 2013; Lu et al. 2014).
The EHT 2017 campaign has yielded total intensity images of the
shadow of the black hole at the centre of M87 at 230 GHz (Event
Horizon Telescope Collaboration et al. 2019a,d,f). Assuming statis-
tically preferred geometric crescent models and general-relativistic
magnetohydrodynamic (GRMHD) models, measurements of phys-
ical properties such as the diameter of the shadow (
42
3
μ
as) and
angular size of one gravitational radius (
3
8
0
4
μ
as) have been
obtained (Event Horizon Telescope Collaboration et al. 2019d,e,f).
These measurements correspond to a mass of
6
5
0
7
10
9
M
at
the estimated distance
16
8
̧
0
8
0
7
Mpc (Event Horizon Telescope Col-
laboration et al. 2019f), consistent with prior mass measurements
based on stellar dynamics (Gebhardt et al. 2011).
Synthetic data play a signi cant role in understanding the char-
acteristics of an instrument, developing new algorithms for data
analysis, and realistically representing the underlying physics that
give rise to the observed data. Feasibility studies and identi ca-
tion of new sites for upgrading existing arrays such as the EHT,
the Karl G Jansky Very Large Array (VLA) (Perley et al. 2011),
European VLBI Network (EVN) (Porcas 2010), East Asian VLBI
Network (EAVN) (An et al. 2018), the Atacama Large Millime-
ter Array (ALMA) (Matthews et al. 2018; Goddi et al. 2019), and
MeerKAT (Jonas 2009), and for building new arrays such as the
next-generation EHT (ngEHT) (Blackburn et al. 2019), the next-
generation VLA (ngVLA) (Selina et al. 2018), and the Square Kilo-
metre Array (SKA) (e.g. Schilizzi et al. 2007) can bene t greatly
from realistic simulations of interferometric observations. These
bene ts include predictive analyses of new hardware enhancements,
as well as testing and optimization of calibration and imaging al-
gorithms and pipelines. Powerful analyses can be performed by the
average user when the synthetic data tools are user-friendly and
can be seamlessly integrated with the calibration and analysis tools.
meqsv2
has been designed with this as a guiding principle.
Roelofs et al. (2020) provide a brief review of the various
synthetic data generation approaches used for simulating EHT ob-
servations over the past decade. Most of these incorporate thermal
noise arising from the characteristics of the instrument and the at-
mosphere as the only data corruption eect. More complex signal
corruptions are introduced by
eht-imaging
(Chael et al. 2016, 2018)
and
meqsilhouette
(Blecher et al. 2017, and this paper), the two
synthetic data generation packages used for generating synthetic
M87 observations in Event Horizon Telescope Collaboration et al.
(2019a,d).
eht-imaging
can introduce randomly varying complex gains
and elevation-dependent atmospheric opacity corruptions, and also
simulates scattering due to the interstellar medium (ISM), which
aects observations of Sgr A
(Johnson 2016; Johnson et al. 2018).
Instrumental polarization is introduced using previous measure-
ments of leakage characteristics of the EHT antennas (Johnson
et al. 2015) and the antenna gains are generated as random osets
sampled from a normal distribution with a standard deviation that is
within a xed percentage of the visibility amplitudes of actual EHT
data.
eht-imaging
also adds randomized station-dependent inter-
scan phases, that are kept coherent within a single scan to mimic
the phases of fringe- tted visibilities.
meqsilhouette
takes a complementary approach to synthetic
data generation by introducing corruptions based on physical mod-
els and tuned to match empirical station measurements. It was de-
signed to adapt the simulation and calibration techniques developed
for metre and cm-wavelength observations to mm-VLBI observa-
tions. The version presented in Blecher et al. (2017) could simulate
simple Gaussian sources and narrow- eld complex sky models and
introduce physically-motivated tropospheric phase and amplitude
corruptions, interstellar scattering using
scatterbrane
1
(Johnson
& Gwinn 2015), and time-variable antenna pointing errors. It has
been signi cantly enhanced since then, rst in step with the pub-
lication of the rst results from 2017 EHT observations, and then
with the development of
symba
(Roelofs et al. 2020) and the rst
polarimetric results of the 2017 M87 observations (Event Horizon
Telescope Collaboration et al. 2021a,b).
meqsilhouette v2
(hereafter
meqsv2
2
) introduces full po-
larimetric, time-variable, spectrally-resolved synthetic data gener-
ation and corruption capabilities and is capable of handling wide-
eld source models with complex substructures.
meqsv2
has been
rewritten and expanded from Blecher et al. (2017). The code has
been refactored to be fully compatible with the same pipelines used
for the analysis of real EHT data. The pointing and atmospheric
models have been rewritten to include more sophisticated eects.
Source and instrumental polarization simulation capabilities have
been introduced.
meqsv2
also accounts for the eects of bandwidth
on various propagation path eects at mm-wavelengths. These fea-
tures facilitate a variety of studies for both the EHT and upcoming
VLBI arrays, such as performing rotation measure (RM) synthesis
studies (Brentjens & de Bruyn 2005) and multi-frequency synthesis
imaging with the increasing fractional bandwith of the EHT, as well
as the envisioned multi-band imaging at 230 GHz and 345 GHz
for the ngEHT. The ability to vary instrumental polarization across
the receiver bandwidth is crucial to take full advantage of ultra-
wideband receivers and high dynamic-range polarimetric imaging.
With the polarimetric primary beam module, full Stokes primary
beam modelling for upcoming arrays such as the ngEHT can be
undertaken. Roelofs et al. (2020) seamlessly combines the func-
tionality of
meqsv2
and
rpicard
(Janssen et al. 2019), the
casa
-
based VLBI pipeline for calibrating data from the EHT and other
VLBI facilities. Synthetic data generated by both
eht-imaging
and
symba
, representing complementary approaches, have been found
to be consistent with each other (Event Horizon Telescope Collab-
oration et al. 2019d).
In this paper, we present the components of
meqsv2
and il-
lustrate its simulation capabilities. In particular, we illustrate the
new polarimetric and spectral resolution capabilities using synthetic
data and validate them. We study the eects of bandpass gains on
closure quantities, which can limit or bias constraints on intrin-
sic source structure asymmetry, by generating synthetic data with
realistic bandpasses. The polarimetric capabilities of
meqsv2
are
1
https://github.com/krosenfeld/scatterbrane
.
2
Measurement EQuation" (see section 2) + Silhouette" (referring to the
black hole shadow").
MNRAS
000
, 116 (2021)
MeqSilhouette v2: Synthetic data generation for the EHT
3
validated using
polsolve
, a
casa
task developed for polarization
calibration of VLBI observations (Martí-Vidal et al. 2021). Un-
like conventional VLBI polarimetric calibration software packages
such as
lpcal
(Leppanen et al. 1995) that use a linear approxi-
mation to model polarization leakage,
polsolve
uses a non-linear
model derived from the full RIME to handle high leakage values
for specialised cases such as the EHT. In addition,
polsolve
uses a
combined multi-source model when the parallactic angle coverage
for individual calibrators is limited and can model the frequency-
dependence of the leakage terms for calibrating large fractional
bandwidths.
This paper is organized as follows: Section 2 provides a brief
overview of the RIME formalism that forms the basis of how
meqsv2
models visibilities to make this a self-contained publi-
cation. Section 3 provides a detailed account of the control ow
of
meqsv2
, its signal corruption capabilities, and their Jones ma-
trix implementations. Section 4 describes the
casa polsolve
tool
for estimating polarization leakage in heterogenous VLBI arrays.
Section 5 demonstrates the synthetic data generation capabilities of
meqsv2
for three dierent mm-wave telescopes. Section 6 quanti-
es the eect of complex bandpass gains on EHT observations at
230 GHz and section 7 demonstrates and validates the polarimetric
simulation capabilities of
meqsv2
using multiple polarized source
models. Section 8 provides a general discussion on the potential
uses of
meqsv2
and section 9 summarises the results and future
outlook.
2 THE RADIO INTERFEROMETER MEASUREMENT
EQUATION
Hamaker et al. (1996) originally developed the mathematical
formalism for describing radio polarimetry using the
Jones
(Jones 1941) and
Mueller
(Mueller 1948) calculi from optics.
Smirnov (2011a) extended this formulation to incorporate direction-
dependent eects (DDEs) in calibration. Here we provide only a
brief summary of the relevant aspects of this formalism and, given
the range of notations in use, establish the notation used in this
paper.
An interferometer produces four pairwise correlations between
the voltage vectors from two stations
푝
and
푞
(each with two feeds
푥
and
푦
), that can be arranged into the 2
2
visibility matrix
3
:
V
푝푞
=
2
h
푣
푝푥
푣
푞푥
i h
푣
푝푥
푣
푞푦
i
h
푣
푝푦
푣
푞푥
i h
푣
푝푦
푣
푞푦
i
(1)
where the angled brackets denote averaging over some small time
and frequency bin, based on considerations of
smearing
and
deco-
herence
, eld of interest, and processing requirements. In terms of
the voltage two-vectors
풗
푝
, equation (1) can be represented as
V
푝푞
=
2
*
푣
푝푥
푣
푝푦
푣
푞푥
푣
푞푦
+
=
2
h
풗
푝
풗
퐻
푞
i
=
2
h
푱
푝
풆
¹
푱
푞
풆
º
퐻
i
=
2
h
푱
푝
¹
풆풆
퐻
º
푱
퐻
푞
i
(2)
where
풆
is the incoming electromagnetic wave,
푱
푝
are the 2
2
Jones
matrices that describe any linear transformation acting on
the incoming wave, and
퐻
is the Hermitian conjugate. The matrix
product
풆풆
퐻
in equation (2) is related to the four Stokes parameters
3
The factor of 2 is introduced in this equation to ensure that the bright-
ness matrix (introduced shortly) becomes
1
for a
1
Jy unpolarized source
(Smirnov 2011a).
퐼
,
푄
,
푈
, and
푉
that describe the polarization state of electromag-
netic radiation (Hamaker et al. 1996; Thompson et al. 2017) by the
following relation
4
:
2
h
푒
푥
푒
푥
i h
푒
푥
푒
푦
i
h
푒
푦
푒
푥
i h
푒
푦
푒
푦
i
=
퐼
̧
푄 푈
̧
푖푉
푈
푖푉 퐼
푄
B
(3)
B
is the
brightness matrix
which describes the intrinsic source
brightness.
푒
푥
and
푒
푦
are the orthogonal polarizations as measured
by the two feeds
푥
and
푦
.
In the ideal case of corruption-free reception, the phase delay
associated with signal propagation, denoted by the scalar
K-Jones
matrix, is always present, giving rise to the
source coherency
,
X
푝푞
:
X
푝푞
=
퐾
푝
B
퐾
퐻
푞
where
퐾
푝
=
푒
2
휋푖휙
푝
1
0
0
1
(4)
in which
휙
푝
denotes the phase delay between the antenna
푝
and the
reference antenna. In the presence of multiple discrete sources in
the sky, taking into account the direction-dependence of the source
coherency and some Jones matrices, the RIME generalizes to
V
푝푞
=
퐺
푝
∑︁
푠
퐸
푠푝
X
푠푝푞
퐸
퐻
푠푞
!
퐺
퐻
푞
(5)
where the summation is carried out over all the sources and
푬
푠푝
and
푮
푝
denote generic direction-dependent eects (DDEs) and
direction-independent eects (DIEs) respectively.
meqsv2
does not
simulate any DDEs for EHT observations, since, aside from scat-
tering, there are no major DDEs that occur along the signal path.
3 THE MEQSILHOUETTE FRAMEWORK
meqsv2
was designed to use the Measurement Set
5
(MS), a database
format designed to store radio astronomical data in next-generation
facilities such as JVLA, ALMA, MeerKAT, and SKA. Fig. 1 shows
the basic layout and the components of a typical
meqsv2
run. Scat-
tering by the ISM is not applied to the input sky models and is
assumed to have been applied externally, simplifying the user in-
terface compared to v1.
meqsv2
uses a
driver
script to set up the
sequence of steps to be executed to generate the synthetic data.
The
framework
module contains the various functions necessary to
create synthetic data, corrupt them, and optionally generate addi-
tional data products. The inputs to the driver script are presented as
attribute-value pairs in a le in the
json
format (Crockford 2006)
containing information on the source, weather, and antenna param-
eters necessary for computing various components of the RIME.
The Jones matrices are applied to the uncorrupted visibilities in
the order in which they occur along the signal path (Noordam &
Smirnov 2010), unless they are scalar, in which case they can be
applied anywhere along the signal chain.
Advanced users may write their own driver scripts to tailor the
basic strategy provided by
meqsv2
for their own needs. For instance,
in
symba
(Roelofs et al. 2020), we use this framework to create
synthetic data that follow real EHT observing schedules using input
VEX
6
les (the scheduling protocol for VLBI experiments) and
4
This equation is valid for linear (XY) feeds. See section 3.1 for circular
(RL) feeds.
5
https://casa.nrao.edu/Memos/229.html
.
6
https://vlbi.org/vlbi-standards/vex
.
MNRAS
000
, 116 (2021)
4
Natarajan et al.
Figure 1.
Flowchart showing the basic components of synthetic data generation with
meqsv2
. The inputs and outputs are shaded orange, the processes are
shaded blue, and the decision boxes (diamonds) are shaded green. Multiple input con guration les are used and various output data products are produced.
The input values are determined from empirical values obtained from the individual observations themselves, as well as
vlbimonitor
. Each component in the
diagram is explained in section 3.
Table 1.
Physical sizes and mount speci cations of EHT2017 stations par-
ticipating in the simulations.
Station
Eective Diameter (m)
Mount type
ALMA
y
70
Alt-az
APEX
12
Alt-az+Nasmyth-Right
LMT
32
Alt-az+Nasmyth-Left
PV
30
Alt-az+Nasmyth-Left
SMT
10
Alt-az+Nasmyth-Right
JCMT
15
Alt-az
SMA
y
16
Alt-az+Nasmyth-Left
y
Single-antenna equivalent of phased arrays.
extend the pointing oset module to compute short and long-term
pointing osets mimicking the behaviour of EHT stations.
symba
also performs additional a priori calibration so that the output more
closely resembles the uncalibrated EHT data.
meqsv2
can just as
easily be used for simulating observations with other VLBI arrays
(see section 5), although other propagation path eects that become
signi cant at longer wavelengths may need to be implemented.
The following subsections explain the various modules of
meqsv2
. The plots shown are obtained using a hypothetical 12-
hour observing run using the EHT2017 array listed in Table 1
at an observing frequency of 230 GHz towards M87 (
훼
J2000
=
12
h
30
m
49
s
42
훿
J2000
=
12
23
0
28
00
04
). The SPT station from
EHT2017 has been excluded since M87 is always below the horizon
from the south pole.
3.1 Input sky models
meqsv2
introduces the capability to generate synthetic visibilities
from wide- eld non-parametric images and retains the ability to
input simple parametric source models in the form of ASCII text
les describing point or Gaussian source models. Since images are
gridded representations of the sky, we use
wsclean
, a fast generic
wide eld imager (Oringa et al. 2014), to Fourier-invert the model
image to the
uv
-plane. Each polarized, time and frequency variable
image frame is Fourier-inverted into the appropriate subset of the
generated MS
7
. Parametric sky models are handled by
meqtrees
(Noordam & Smirnov 2010), which performs a direct Fourier trans-
form of the sky model into the MS.
For the EHT, the visibilities are always computed in the circular
polarization basis (RR, RL, LR, and LL) except in the case of ALMA
which records signals in linear basis. In
meqsv2
, we assume that
the ALMA visibilities have been perfectly converted to circular
polarization (Marti-Vidal et al. 2015) so that the basis is uniform
across all stations. Then, equation (3) takes the form
B
=
2
h
푒
푟
푒
푟
i h
푒
푟
푒
푙
i
h
푒
푙
푒
푟
i h
푒
푙
푒
푙
i
=
퐼
̧
푉
푄
̧
푖푈
푄
푖푈 퐼
푉
(6)
where `
' indicates circular polarization (Smirnov 2011a). Figure 2
shows the Stokes I visibility amplitudes for all baselines for a single
channel of the data set generated using a point source sky model
with intrinsic time variability. The scatter of the visibilities is due
to thermal noise. The capability to simulate time-varying sources is
particularly useful for simulating sources exhibiting variability on
timescales of minutes to hours, such as the radio source associated
with the supermassive black hole Sagittarius A
(or Sgr A
) located
at the Galactic centre (e.g. Lu et al. 2016). In addition, studies on
decoupling time-varying instrumental eects from source evolution
could be undertaken.
Figure 3 shows the Stokes
퐼
visibilities of a point source with
a steep spectral index, across a 2 GHz bandwidth divided into 64
channels, centred at 227 GHz. The data shown correspond to a 5-
minute subset of the entire observation for all channels. As before,
the scatter seen is due to thermal noise.
7
The naming conventions for the image frames are explained in the o cial
documentation.
MNRAS
000
, 116 (2021)
MeqSilhouette v2: Synthetic data generation for the EHT
5
Figure 2.
Stokes
퐼
visibility amplitudes for a point source with intrinsic
time variability, as a function of time for one frequency channel.
Figure 3.
Stokes
퐼
visibility amplitudes of a point source with intrinsic
frequency variability, as a function of frequency for 5 minutes of observation.
3.2 Antenna pointing errors
Several factors cause antennas to mispoint and modify its gain re-
sponse, such as dierential exure of the antenna, wind pressure,
and thermal expansion. Moreover, errors in drive mechanics or the
telescope control software also introduce pointing errors, which
attenuate the measured visibility amplitudes,
j
푉
푝푞
j
. Even small o-
sets in antenna pointing can cause signi cant reduction in visibility
amplitudes at mm-wavelengths.
Formulated in terms of the RIME, this eect is represented
by a time and direction-dependent antenna-based E-Jones matrix
(Smirnov 2011b):
퐸
푝
¹
푙푚
º
=
퐸
¹
푙
̧
훿푙
푝
푚
̧
훿푚
푝
º
(7)
where
훿푙
and
훿푚
are the time-variable osets in the
푙푚
-plane.
meqsv2
provides a WSRT-derived analytic
cos
3
beam (from the
previous version) and a Gaussian primary beam pro le. The Gaus-
sian pro le has the advantage that it can be conveniently described
by a single parameter, the full width at half maximum (FWHM),
and does not deviate from the more complex Bessel function from
the centre up to very close to the rst null (Middelberg et al. 2013).
This assumption is justi ed in the context of EHT, since, in most
cases, the eld of interest is much smaller than the location of the
rst null.
Figure 4.
The true primary beam response pattern for the EHT stations at
230 GHz, as a function of the pointing oset,
휌
. Each point corresponds to
one realisation of the pointing oset.
The osets per pointing epoch (or scan) for station
푝
,
휌
푝
,
are drawn from the normal distribution
N¹
0
P
rms
푝
º
.
P
rms
values
depend on station and weather characteristics and are determined
based on empirical measurements. The full width at half maximum
(FWHM) of the primary beam
P
FWHM
푝
of each antenna at 230
GHz is scaled to the centre frequency of the observation. The beam
model with the pointing errors is given by
퐸
푝
=
exp
1
2
"
휌
푝
¹P
FWHM
푝
2
p
2 ln 2
º
#
2
!
where
휌
푝
=
√︃
훿푙
2
푝
̧
훿푚
2
푝
(8)
The
2
p
2 ln 2
factor arises due to the relationship
FWHM
=
2
p
2 ln 2
휎
, where
휎
is the standard deviation. Updating equation
(8) is all that is required to implement more complex, asymmetric
primary beam patterns. The corresponding E-Jones matrix is given
by
퐸
푝
=
퐸
푝푎
¹
푙푚
º
0
0
퐸
푝푏
¹
푙푚
º
(9)
This is implemented as a scalar term (
퐸
푝푎
=
퐸
푝푏
) that can com-
mute with other components of the RIME.
Figure 4 shows an example simulation of pointing osets,
휌
푝
,
and how they aect the primary beam response of the EHT2017
antennas. SMA is the least aected owing to its large beam size,
while LMT is the most aected due its narrow primary beam.
3.3 Tropospheric eects
The troposphere has a signi cant eect on signal propagation at
mm wavelengths (Carilli & Holdaway 1999).
meqsv2
models three
signi cant tropospheric eects (i) the signal attenuation due to
atmospheric opacity
휏
, (ii) the increased system temperature due to
atmospheric emission at the sky brightness temperature
푇
sky
, and
(iii) the phase uctuations due to atmospheric turbulence. These
eects are separated into mean and turbulent components, with the
mean component further divided into wet (due to water vapour)
and dry components. A detailed treatment of the propagation fun-
damentals can be found in Blecher et al. (2017). Here we provide a
summary along with the updated equations and algorithms imple-
mented in
meqsv2
.
MNRAS
000
, 116 (2021)
6
Natarajan et al.
Figure 5.
Variation in transmission with time at the sites of ALMA and
LMT, during a 7-hour tracking of M87.
3.3.1 Mean troposphere
At mm-wavelengths, the mean component of the troposphere causes
signal attenuation due to absorption and time delays, along with a
pseudo-continuum opacity that increases with frequency. The H
2
O
and O
2
lines cause signi cant absorption above 100 GHz, peaking
in the THz regime (Thompson et al. 2017).
meqsv2
retains the use
of the Atmospheric Transmission at Microwaves (
atm
, Pardo et al.
2001) software from Blecher et al. (2017) to compute the mean
opacity, sky noise, and the mean component of the delays.
atm
implements radiative transfer through a static atmosphere. In the
absence of scattering, the radiative transfer equation for unpolarized
radiation is given by
d
퐼
휈
¹
푠
º
d
푠
=
휖
휈
¹
푠
º
휅
휈
¹
푠
º
퐼
휈
¹
푠
º
(10)
where
푠
is the coordinate along the direction of the signal path
through the atmosphere,
퐼
휈
¹
푠
º
is the speci c intensity at frequency
휈
, and
휖
휈
and
휅
휈
are respectively the emission and absorption coe -
cients of the atmosphere. The absorption coe cient
휅
휈
is calculated
as a function of frequency and altitude, enabling equation (10) to
be integrated along the signal path to obtain mean opacity and sky
brightness temperature.
An example of the elevation-dependent transmission for
ALMA and LMT over a 7-hour observing track towards the di-
rection of M87 is shown in Figure 5. The mean time delay is ob-
tained from the absorption coe cient using the Kramers-Krönig
relations (e.g. de L. Kronig 1926), for which the necessary atmo-
spheric temperature and pressure pro les are calculated based on
site-dependent precipitable water vapour (PWV) depth and ground
temperature and pressure (Pardo et al. 2001).
The a priori correlator model used by the EHT corrects for
various eects including the elevation-dependent variations in the
intra-scan delays (Event Horizon Telescope Collaboration et al.
2019b). Optionally,
meqsv2
simulates this correction by introducing
a constant residual delay per scan, that is equal to the mean delay
averaged over all frequency channels.
Integrating equation (10) also yields sky brightness tempera-
ture, which contributes to an increase in the system temperature.
The resultant increase in the sky noise is implemented as an opacity
and elevation-dependent noise component. This noise, along with
the thermal noise component, is accounted for in the total noise
budget used for visibility weighting in
meqsv2
(see section 3.6).
3.3.2 Turbulent troposphere
Turbulence in the troposphere introduces phase instabilities in the
measured visibilities. The spatial distribution of water vapour in
the troposphere evolves rapidly, reducing the coherence timescale
to
10
s (Thompson et al. 2017). This limits the coherent averaging
time (and hence the S/N) of uncalibrated visibilities, and renders
conventional calibration procedures (with interleaved observations
of calibrator and target source) ineective. Self-calibration is also
limited by the S/N required for fringe- tting.
meqsv2
models these phase variations,
훿휙
¹
푡휈
º
, by assuming
a thin, Kolmogorov-turbulent phase screen moving with constant
velocity (e.g. Johnson & Gwinn 2015). This method is applicable
to any situation in which the troposphere induces delays in elec-
tromagnetic wave propagation (e.g. Treuhaft & Lanyi 1987). The
phase osets introduced can be described by the phase structure
function given by (Carilli & Holdaway 1999)
D
휙
¹
풙
풙
0
º
=
h»
휙
¹
풙
̧
풙
0
º
휙
¹
풙
º¼
2
i
(11)
where
풙
and
풙
0
are points on the screen and the angled brackets
denote ensemble average. This equation can be reasonably approx-
imated by a power law (e.g. Armstrong et al. 1995):
D
휙
¹
푟
º
푟
푟
0
훽
(12)
where
푟
2
=
¹
풙
풙
0
º
2
and
푟
0
is the phase coherence scale at which
D
휙
¹
푟
0
º
=
1 rad
.
Scattering can be classi ed as
strong
or
weak
based on the rela-
tionship between
푟
0
and the Fresnel scale,
푟
퐹
, de ned as
√︁
휆퐷
os
2
휋
,
where
퐷
os
is the distance between the observer and the scattering
screen. The radiative power measured at a given point originates
from a single region of area
퐴
weak
=
휋푟
2
퐹
when
푟
퐹
푟
0
and
from multiple disconnected zones each of area
퐴
strong
=
휋푟
2
0
when
푟
퐹
푟
0
(Narayan 1992). Empirical estimates of
푟
0
fall in the range
of
50
500
m above Mauna Kea (Masson 1994) and
90
700
m
above Chajnantor (Radford & Holdaway 1998). At
휆
=
1
3
mm and
퐷
os
=
2
km (the water vapour scale height),
푟
퐹
0
64
m (
푟
0
),
and hence scattering falls under the weak regime.
Equation (12) may be rewritten with explicit time-dependence
as
퐷
휙
¹
푡
º
=
퐷
휙
¹
푟
ºj
푟
=
풗
푡
, where
풗
is the bulk transverse velocity
of the phase screen (Coulman 1985). Assuming that the coherence
timescale when observing towards the zenith,
푡
c
=
푟
0
j
풗
j
(Treuhaft
& Lanyi 1987; Blecher et al. 2017), we get
D
휙
¹
푡
º
푡
푡
푐
훽
(13)
Both Kolmogorov theory and empirical measurements show that
the exponent
훽
should be equal to
5
3
when
푟
Δ
ℎ
, where
Δ
ℎ
is the thickness of the turbulent layer which is
1
km (Carilli &
Holdaway 1997). The processed Field-of-View (FoV) of the EHT
is about
100
100
휇
as
2
and is much smaller than
푟
0
, allowing us
to represent it as a diagonal Jones matrix in the RIME.
The antenna-based turbulent phase errors manifest as a series
of correlated, normally distributed random variables in time. This
time-series, {
훿휙
0
¹
푡
º
}, can be constructed as follows (e.g. Rasmussen
& Williams 2006). From the structure function
D
we construct the
covariance matrix
Σ
. Since
Σ
is symmetric and positive de nite,
MNRAS
000
, 116 (2021)
MeqSilhouette v2: Synthetic data generation for the EHT
7
the lower triangular matrix,
퐿
, resulting from the Cholesky decom-
position of
Σ
(where
Σ =
퐿퐿
푇
) can be applied to a time-series of
the desired length with zero mean and unit variance to arrive at
correlated random samples.
We assume a simple linear scaling with frequency across the
bandwidth since the wet dispersive path delay is not more than a
few per cent of the non-dispersive component at mm-wavelengths
(Curtis et al. 2009). Also taking into account the airmass towards
the horizon when observing away from the zenith, the phase error
time-series for antenna
푝
becomes
f
훿휙
푝
¹
푡휈
ºg
=
1
√︁
sin
¹
휃
el
¹
푡
ºº
f
훿휙
0
푝
¹
푡
ºg
휈
휈
0
(14)
where
휈
is the list of channel frequencies,
휈
0
is the reference fre-
quency, taken to be the lowest channel frequency, and
휃
el
is the
elevation angle.
Since VLBI stations are typically located hundreds or thou-
sands of kilometres from each other, the tropospheric corruptions
over individual stations are uncorrelated with each other and must
be simulated independently. This is not strictly true for the short-
baseline pairs of ALMA-APEX (2.6 km) and JCMT-SMA (160 m)
in the EHT, for which the turbulence may be correlated since the
baseline lengths are so short as to be comparable to
푟
0
(e.g. Carilli
& Holdaway 1997). Currently, intra-site baseline correlations are
not simulated and the phase errors are generated independently for
these stations.
3.4 Instrumental polarization
The two polarization feeds nominally measure two orthogonal po-
larization states of the incoming wave, in either circular (RL) or
linear (XY) bases. In practice, mechanical imperfections in the feed
or electronic imperfections in the signal path cause signals from
each independent signal path to leak into the other (Hamaker et al.
1996; Sault et al. 1996). Additionally, the possible rotation of the
feeds by the parallactic angle depending on the antenna mount type
must be taken into account.
3.4.1 Parallactic angle rotation
The mount type of the antenna determines the rotation of the feeds
as seen from the sky (e.g. Dodson 2009). For equatorially mounted
antennas, the orientation of the primary beam with respect to the
sky remains static. For alt-az mounted antennas, the beam rotates
as the source is tracked over the sky, and the feeds rotate by the
parallactic angle
(
휒
푝푎
) given by
휒
푝푎
=
arctan
sin
¹
Θ
º
cos
¹
푙
º
cos
¹
훿
º
sin
¹
푙
º
cos
¹
Θ
º
cos
¹
푙
º
sin
¹
훿
º
(15)
where
Θ
is the hour angle,
훿
is the Declination, and
푙
is the lati-
tude. For alt-az mounted antennas with Nasmyth focus, the more
generic
feed angle
(
휒
푝
) is used in place of parallactic angle, which
incorporates the elevation of the antenna. Depending on whether
the tertiary mirror focuses the light to the right or to the left, the
elevation angle (
휒
푒푙
) is added or subtracted respectively,
휒
푝
=
휒
푝푎
휒
푒푙
where
휒
푒푙
=
arcsin
sin
¹
푙
º
sin
¹
훿
º ̧
cos
¹
푙
º
cos
¹
훿
º
cos
¹
Θ
º
(16)
Figure 6.
Evolution of the parallactic angles of EHT 2017 array stations
over 12-hours in the direction of M87.
In terms of Jones calculus, the
P-Jones
matrix for parallactic angle
rotation for circular polarization is given by (e.g. Hales 2017)
푃
푝
=
e
푗휒
푝
0
0
e
푗휒
푝
(17)
Figure 6 shows the evolution of the parallactic angle corresponding
to the EHT stations listed in Table 1 over the course of 12 hours
towards the direction of M87. The parallactic angle values of the
short baselines formed by (i) JCMT-SMA and (ii) ALMA-APEX
evolve similarly due to the proximity of these stations to each other.
3.4.2 Polarization leakage
The two receptors in the feed are designed to be sensitive to orthog-
onal polarization states. Ideally, the
D-Jones
terms (or the
D-terms
)
corresponding to polarization leakage is represented by a unit ma-
trix, expressed in the appropriate coordinate system (i.e. that of the
two polarization states measured by the feed receptors) (Hamaker
et al. 1996). In practice, the non-diagonal terms of this matrix are
non-zero, due to the fact that each receptor is sensitive to the oppo-
site polarization state due to electronic or mechanical imperfections
in the feed. Hence, the
feed matrix
is given by
퐷
0
푝푅
푑
0
푝푅
푑
0
푝퐿
퐷
0
푝퐿
!
(18)
where
푅
and
퐿
denote the two feed receptors (here, in the circular
frame). This feed matrix may be decomposed as
퐷
0
푝푅
0
0
퐷
0
푝퐿
!
1
푑
푝푅
푑
푝퐿
1
=
퐺
0
푝
퐷
푝
where
푑
푝푅
=
푑
0
푝푅
퐷
0
푝푅
(19)
The matrix
퐺
0
푝
may be absorbed in the
G-Jones
terms that represent
the receiver gains (see section 3.5). The matrix
퐷
푝
is the feed error
or the leakage
D-Jones
term, with the complex numbers
푑
푝푅
and
푑
푝퐿
representing the fractional leakage from either feed.
3.4.3 Implementation
meqsv2
can introduce station-based, frequency-dependent,
complex-valued instrumental polarization. The per station complex
D-term values for the two polarization feeds are sampled inde-
pendently for each frequency channel from the normal distribution
MNRAS
000
, 116 (2021)
8
Natarajan et al.
N¹
푑
푝휇
푑
푝휎
º
;
푑
푝휇
denotes the characteristic empirical leakage
value and
푑
푝휎
denotes the scatter for station
푝
. The visibilities may
be written to the MS in either the sky or the antenna coordinate
system (i.e. without or with parallactic angle rotation correction re-
spectively). The visibilities in the two coordinate systems are related
by
V
ant
=
푃
푝
V
sky
푃
퐻
푞
(20)
where
푃
푝
are given by equation (17). This rotation also includes a
constant oset in the feed angle per station. Once this rotation has
been applied, the D-terms can be applied to the visibilities in the
antenna frame,
V
=
퐷
푝
¹
푃
푝
º
V
sky
¹
푃
퐻
푞
º
퐷
퐻
푞
(21)
If the visibilities are required to be written in the sky frame,
then
V
sky
is converted to
V
ant
so that the D-terms may be applied,
before being converted back to the sky frame,
V
=
¹
푃
퐻
푝
º
퐷
푝
¹
푃
푝
º
V
sky
¹
푃
퐻
푞
º
퐷
퐻
푞
¹
푃
푞
º
(22)
The product
푃
퐻
푝
퐷
푝
푃
푝
evaluates to a rotation of the D-terms in the
antenna frame by twice the feed angle:
푃
퐻
푝
퐷
푝
푃
푝
=
1
exp
¹
2
푗휒
푝
º
푑
푅
푝
exp
¹
2
푗휒
푝
º
푑
퐿
푝
1
(23)
The visibilities are usually generated in this frame to correspond to
the EHT data.
3.5 Temporal and spectral variability in receiver gains
The antenna-based gain terms are represented by the complex-
valued
G-Jones
matrices which are a generalization of simple an-
tenna gains to polarimetry (e.g. Smirnov 2011a). The electronic
gains of a pair of circular receptors are given by a diagonal G-Jones
matrix in the circular coordinate frame
퐺
푝
¹
푡
º
=
푔
푝푅
¹
푡
º
0
0
푔
푝퐿
¹
푡
º
(24)
Dierent antenna-based factors such as the
퐺
0
term in equation
(19) may be subsumed into the G-Jones terms. The time-variable
complex gains are generated by sampling from the normal distribu-
tion
N¹
푔
푝휇
푔
푝휎
º
per timestamp, for both the polarization feeds
for each station
푝
, where
푔
푝휇
denotes the characteristic station-
dependent gain and
푔
푝휎
denotes the scatter.
Various components along the signal path and the bandpass
lters used at each station determine the frequency response of
the receiving channels over the observing bandwidth (Taylor et al.
1999). This response is not constant over the entire bandwidth and
usually falls to about half of the maximum towards both edges
of the passband. This results in a frequency-dependent component
in the antenna gains which is captured using complex-valued
B-
Jones
matrices that vary as a function of frequency. This eect
is usually corrected for by observing a bandpass calibrator source
whose behaviour across the observing bandwidth is well-known.
As with the G-Jones terms, the B-Jones matrices for circular
receptors are diagonal in a circular coordinate frame,
퐵
푝
¹
휈
º
=
푏
푝푅
¹
휈
º
0
0
푏
푝퐿
¹
휈
º
(25)
For each station,
meqsv2
accepts nominal gain values at a few
representative frequencies within the bandwidth of the observation,
and performs spline interpolation of the bandpass amplitudes for
each frequency channel. To these amplitudes, random phases are
generated and added to produce the complex quantities. The G-Jones
and B-Jones terms are simulated separately to provide independent,
ne-grained variability of gains along the time and frequency axes.
3.6 Noise considerations
3.6.1 Thermal noise
Noise due to various factors such as receiver electronics, atmo-
spheric emission, background radiation and ground (i.e. spillover)
radiation aect the system sensitivity adversely. The system tem-
perature,
푇
s
푦푠
, equivalent to the power per unit frequency due to the
noise,
푃
푁
, accounts for this noise and is given by
푇
s
푦푠
=
푃
푁
푘
퐵
(26)
where
푘
퐵
is the Boltzmann constant
8
. The system temperature is of-
ten expressed in terms of the system equivalent ux density (SEFD),
which is de ned as the ux density of a source that would deliver
the same amount of power as the system
SEFD
=
2
푘
B
푇
sys
휂퐴
e
(27)
where
휂
is the antenna e ciency and
퐴
e
is the eective area of
the antenna.
meqsv2
accepts the station-dependent SEFD values
as inputs from which the per-baseline rms uncertainty,
휎
푝푞
, on a
visibility in units of Jy is computed (Thompson et al. 2017):
휎
푝푞
=
1
휂
√︂
SEFD
푝
SEFD
푞
2
Δ
휈휏
(28)
where
퐴
e
denotes the eective area of the telescope and
휂
comprises
any relevant e ciency terms, such as the antenna aperture e ciency,
휂
ap
, the correlator e ciency,
휂
corr
,
Δ
휈
is the bandwidth and
휏
the
integration time. For standard 2-bit quantization,
휂
is set to 0.88.
Since the system noise is broadband and mostly stationary, it can
be described using a Gaussian distribution and the uncertainty in
their measurements can be reduced by increasing the number of
independent measurements
푁
=
2
Δ
휈휏
.
In the RIME implementation, this thermal noise becomes an
additive term per polarization, distributed normally with zero mean
and a variance of
휎
2
푝푞
per visibility:
V
푝푞
=
V
0
푝푞
̧N¹
0
휎
2
푝푞
º
(29)
This additive thermal noise matrix has the same dimensionality as
the data, varying with time, baseline, frequency, and polarization.
3.6.2 Visibility weighting
The MS format allows for the estimated rms noise values and visibil-
ity weights to be recorded alongside the data. The
휎
푝푞
values com-
puted above are used to generate per-visibility baseline-dependent
thermal noise. These terms are added with the increase in the sky
brightness temperature due to the troposphere (section 3.3.1) and
are used to ll the SIGMA and SIGMA_SPECTRUM columns in
the MS. Inverse-variance weighting is used to ll in the visibility
weights columns WEIGHT and WEIGHT_SPECTRUM.
8
The tropospheric contribution to the increase in
푇
sys
(section 3.3.1) is
added to the noise budget as mentioned below.
MNRAS
000
, 116 (2021)
MeqSilhouette v2: Synthetic data generation for the EHT
9
4 POLSOLVE: POLARIZATION LEAKAGE
ESTIMATION
We use
polsolve
to validate the instrumental polarization capa-
bilities of
meqsv2
.
polsolve
is part of
poltools
9
, the polarimetry
toolbox developed for
casa
, aimed at the simulation, calibration, and
basic analysis of polarimetric VLBI observations (Martí-Vidal et al.
2021). It uses the full RIME, simpli ed for the case of narrow- eld
observations, to estimate and correct for instrumental polarization
using observations of spatially-resolved polarization calibrators.
polsolve
has many advantages compared to the
aips
10
task
lpcal
(Leppanen et al. 1995), the main algorithm used by the VLBI
community for polarimetric calibration. It uses a non-linear model
of polarization leakage derived from the full RIME for handling high
leakage values.
polsolve
can also model the frequency-dependent
variations in D-terms to enable calibration of wide fractional band-
widths and perform D-term estimates based on cross-polarization
self-calibration. We take advantage of many of these features in
performing the validation tests for this paper.
The
aips
-based
gpcal
11
pipeline also addresses many of the
shortcomings of
lpcal
(Park et al. 2021) and has been shown to be
consistent with
polsolve
(Event Horizon Telescope Collaboration
et al. 2021a). A detailed discussion of the impact of these features
on the quality of VLBI polarimetry are discussed in Martí-Vidal
et al. (2021). Below, we give a brief description of the procedure
and equations used by
polsolve
.
4.1 The POLSOLVE tting parameters
polsolve
uses the Levenberg-Marquardt algorithm (Press 2007) to
minimize the error function
휒
2
¹
풙
º
, where
풙
consists of parame-
ters that model both source and instrumental polarization, divided
into two subsets. It uses two dierent approaches to account for the
polarization structure of a source, both of which assume that the
source brightness distribution can be described as a linear combi-
nation of source components (e.g., point sources, for the case of
clean
deconvolution, Högbom 1974).
The rst approach divides the set of CLEAN components
f
퐼
푖
g
into
푁
푠
subsets, also called source subregions, for which a con-
stant fractional polarization is assumed. The corresponding model
visibility functions for Stokes
푄
and
푈
are given by
V
푄
=
푁
푠
∑︁
푠
푞
푠
푛
푠
∑︁
푘
F¹
퐼
푠
푘
º
!
and V
푈
=
푁
푠
∑︁
푠
푢
푠
푛
푠
∑︁
푘
F¹
퐼
푠
푘
º
!
(30)
where
푁
푠
is the number of subregions with constant fractional
polarization,
푛
푠
is the number of CLEAN components inside the
푠
-th subregion, and
퐼
푠
푘
is
푘
-th CLEAN component of Stokes
퐼
in
the
푠
-th subregion. The vectors
풒
=
f
푞
푠
g
and
풖
=
f
푢
푠
g
are the
parameters that
polsolve
ts for.
The second approach uses independent, and xed (
푁
푠
=
0
),
source models for the brightness distributions of
푄
and
푈
. The
model visibilities for
푄
and
푈
then take the simple form
V
푄
=
푁
푞
∑︁
푖
F¹
푄
푖
º
and V
푈
=
푁
푢
∑︁
푖
F¹
푈
푖
º
(31)
where
f
푄
푖
g
and
f
푈
푖
g
are the (e.g.
clean
) components used to model
9
https://launchpad.com/casa-poltools
.
10
http://www.aips.nrao.edu
.
11
https://github.com/jhparkastro/gpcal
.
the polarization source structure. In this case,
polsolve
performs a
cross-polarization self-calibration estimate of the D-terms.
Instrumental polarization is represented by two complex quan-
tities per antenna, corresponding to each polarization. The leakage
terms obey the equation
V
푖
푅푅
V
푖
푅퐿
V
푖
퐿푅
V
푖
퐿퐿
=
1
퐷
푝
푅
퐷
푝
퐿
1
V
푢푣
퐼
V
푢푣
푃
V
푢푣
푃
V
푢푣
퐼
1
¹
퐷
퐿
º
푞
¹
퐷
푅
º
푞
1
(32)
where
V
푖
푘
is the
푖
-th observed visibility of polarization product
푘
(where
푘
is one of
푅푅
,
푅퐿
,
퐿푅
, or
퐿퐿
) and
푢
and
푣
are the
coordinates in Fourier space. We construct the functions for the
complex polarization vector,
푷
, in the form
V
푃
=
V
푄
̧
푗
V
푈
and V
푃
=
V
푄
푗
V
푈
(33)
These visibilities are assumed to be fully calibrated for atmospheric
eects and electronic antenna gains (computed in the frame of the
antenna mounts).
5 SIMULATING MM-WAVE OBSERVATIONS
meqsv2
was developed primarily for generating synthetic data for
the EHT at mm-wavelengths, but it can equally well be applied to
any array, including proposed arrays such as ngEHT and ngVLA.
For instance,
meqsv2
can be used to perform a more elaborate
exploration of the VLBI capabilities of MeerKAT for performing
extragalactic surveys presented in Deane (2016).
meqsv2
has also
been used to generate 5 GHz synthetic EVN observations for sim-
ulating phase corruptions aecting astrometric uncertainties (van
Langevelde et al. 2019). Roelofs et al. (2020) use
meqsv2
for gener-
ating uncalibrated synthetic data observed using an extended EHT
array with additional stations located at potential future sites.
We generate synthetic data with
meqsv2
for three mm-
wavelength arrays: EHT2017 array, ngVLA in its long baseline
array con guration, and ALMA array in its extended con gura-
tion. The EHT2017 array consists of the stations shown in Table
1 observing at a frequency of 230 GHz. The ngVLA array con-
sists of all 18 stations in the Long Baseline Array (LBA) and the
central core reduced to a single high-sensitivity site, observing at
86 GHz. ALMA consists of 43 12-metre diameter antennas in its
most extended con guration, observing at 230 GHz. An asymmet-
ric crescent (Kamruddin & Dexter 2013) created with
eht-imaging
,
with the inner radius oset by 3
휇
as in the horizontal direction was
used as the source model (Figure 7,
right
).
Various propagation path eects described in the previous sec-
tion such as pointing osets, tropospheric eects, receiver gains,
polarization leakage, and thermal noise are introduced. Nominal
pointing osets based on a priori station information are used (e.g.
Roelofs et al. 2020). The aperture e ciencies were determined
using targeted observations of planets of known brightness temper-
atures as calibrators (Event Horizon Telescope Collaboration et al.
2019c). The individual station SEFDs are determined by extrapo-
lating measured system temperatures to zero airmass (Janssen et al.
2019; Roelofs et al. 2020).
We adopt the median values measured by
vlbimonitor
12
at the
individual sites during the 2017 EHT observing campaign (Event
12
https://bitbucket.org/vlbi
.
MNRAS
000
, 116 (2021)
10
Natarajan et al.
Figure 7.
Symmetric ring (
left
) and asymmetric crescent (
right
) models
from Kamruddin & Dexter (2013) used in generating synthetic data for
sections 5 and 6.
Horizon Telescope Collaboration et al. 2019b) for the weather pa-
rameters. For stations for which this information does not exist,
it was obtained through climatological modelling using the data
sets from Modern-Era Retrospective Analysis for Research and Ap-
plications, version 2 (MERRA-2) from the NASA Goddard Earth
Sciences Data and Information Services Center (GES DISC; Gelaro
et al. 2017) and the
am
atmospheric model software (Paine 2019).
The atmospheric conditions over all the ALMA antennas are the
same, though independently simulated for each antenna.
The top row of Figure 8 shows the
uv
-coverages for the three
arrays in the direction of M87 at their respective observing fre-
quencies, showing how the dierent arrays complement each other.
The bottom row of Figure 8 shows both the corrupted and uncor-
rupted Stokes
퐼
visibility amplitudes for each array as a function of
uv
-distance. The 60
휇
as crescent is fully resolved by EHT2017 as
expected. At 86 GHz, ngVLA resolves it partially, clearly seeing the
source as having some extended structure. ALMA, with its longest
baselines of only 12 km, sees the crescent as a point source at 230
GHz and the model visibilities show up as a straight line at 5 Jy
in the gure. In the following sections, we concentrate only on the
EHT2017 array.
6 EFFECT OF BANDPASS ON CLOSURE QUANTITIES
In this section, we estimate the magnitude of systematic non-closing
errors introduced by complex bandpass gains in EHT observations
at 230 GHz using synthetic data. We introduce station-based band-
passes of varying shapes, spline-interpolated to all frequency chan-
nels as shown in Figure 9. The phases are sampled conservatively
within the uniform range of
30
(Michael Janssen, private comm).
Two dierent sky models from (Kamruddin & Dexter 2013), gen-
erated using
eht-imaging
were used: a symmetric ring with outer
radius 30
휇
as and width 3
휇
as, and an asymmetric crescent with
the inner radius oset by 3
휇
as in the horizontal direction (Figure
7). The total ux of the ring was set to 5 Jy to ensure that the S/N
of the data was high so that systematic non-closing errors could
be unambiguously detected. Thermal noise generated using equa-
tion (28) was applied to all datasets. All synthetic data sets used
in this subsection are generated with a time resolution of 1 s and
a frequency resolution of 31.25 MHz per frequency channel, with
64 channels, spanning a bandwidth of 2 GHz. These data are then
band-averaged to a single channel of width 2 GHz. We generated
four categories of synthetic data, with 10 data sets to each category,
diering in the realisation of thermal noise:
(i) symmetric ring with only thermal noise corruption,
(ii) symmetric ring with thermal noise and bandpass corruption,
(iii) asymmetric crescent with only thermal noise corruption,
and
(iv) asymmetric crescent with thermal noise and bandpass cor-
ruption.
Closure quantities are formed around a closed loop of stations
to eliminate the eects of station-based gain errors (e.g. Blackburn
et al. 2020). Closure phases are formed over a closed triangle
푝푞푟
as
휓
퐶푝푞푟
=
arg
¹
V
푝푞
V
푞푟
V
푟푝
º
(34)
and closure amplitudes are formed over a quadrangle of stations
푝푞푟푠
:
ln
퐴
퐶푝푞푟푠
=
ln
V
푝푞
V
푟푠
V
푝푟
V
푞푠
(35)
where ln is the natural logarithm. The closure phases (CP) and
the log closure amplitudes (LCA) are susceptible to systematic non-
closing errors. The intra-site baselines ALMA-APEX and JCMT-
SMA enable one to form trivial closure quantities, between dier-
ent combinations of these four stations. The trivial CP and LCA so
formed should ideally be zero, but factors such as band-averaging
without correcting for bandpass errors and intrinsic large-scale
asymmetry in source structure can break the assumptions of trivial
closure.
We evaluate the magnitude of systematic non-closing errors
in trivial closure quantities following the EHT data validation pro-
cedures outlined in Event Horizon Telescope Collaboration et al.
(2019c) and in Wielgus et al. (2019). We employ the
MAD
0
(me-
dian absolute deviation from zero) statistic, which, for a normally
distributed variable
푌
with zero mean takes the form
MAD
0
¹
푌
º
=
1
4826 median
¹j
푌
jº
(36)
where the subscript 0 denotes that the raw moment is estimated.
The scaling factor 1.4826 ensures that
MAD
0
acts as an unbiased
estimator of standard deviation for normally distributed data. The
total uncertainty,
휎
, associated with the trivial closure quantities
de ned above is given by
휎
2
=
휎
2
th
̧
푠
2
(37)
where
휎
th
is the known thermal component and
푠
denotes the sys-
tematic non-closing error, modeled as a constant. For a trivial clo-
sure quantity
푋
, we solve for the characteristic value of
푠
by enforc-
ing a following condition:
MAD
0
푋
휎
!
=
MAD
0
푋
√︃
휎
2
th
̧
푠
2
!
=
1
(38)
The
MAD
0
values estimated using equation (36) and the char-
acteristic magnitude of systematic errors
푠
calculated with equation
(38) for the four data sets are given in Table 2. The uncertainties are
estimated by obtaining these values for 10 data sets with dierent
thermal noise realisations.
If the error budget is exactly accounted for, then the
MAD
0
estimator is equal to 1. This is evident from the estimated
MAD
0
values being very close to unity when thermal noise is the only cor-
ruption introduced. The corresponding systematic errors are about
0
08
for trivial closure phases and less than
0
4
per cent for trivial
MNRAS
000
, 116 (2021)