2004.01161.pdf

Astronomy

&

Astrophysics

manuscript no. ms

ESO 2020

April 3, 2020

SYMBA: An end-to-end VLBI synthetic data generation pipeline

Simulating Event Horizon Telescope observations of M87

F. Roelofs

1 ,

, M. Janssen

1 ,

, I. Natarajan

, R. Deane

, J. Davelaar

, H. Olivares

, O. Porth

, S. N. Paine

K. L. Bouman

, R. P. J.

Tilanus

, I. M. van Bemmel

, H. Falcke

, K. Akiyama

, A. Alberdi

, W. Alef

, K. Asada

, R. Azulay

, A. Baczko

, D.

Ball

, M. Balokovi

, J. Barrett

, D. Bintley

, L. Blackburn

, W. Boland

, G. C. Bower

, M. Bremer

, C. D. Brinkerink

, R.

Brissenden

, S. Britzen

, A. E. Broderick

, D. Broguiere

, T. Bronzwaer

, D. Byun

, J. E. Carlstrom

, A. Chael

, C.

Chan

, S. Chatterjee

, K. Chatterjee

, M. Chen

, Y. Chen

, I. Cho

, P. Christian

, J. E. Conway

, J. M. Cordes

, G. B. Crew

, Y.

Cui

, M. De Laurentis

, J. Dempsey

, G. Desvignes

, J. Dexter

, S. S. Doeleman

, R. P. Eatough

, V. L. Fish

, E. Fomalont

, R.

Fraga-Encinas

, P. Friberg

, C. M. Fromm

, J. L. Gómez

, P. Galison

, C. F. Gammie

, R. García

, O. Gentaz

, B. Georgiev

, C.

Goddi

, R. Gold

, M. Gu

, M. Gurwell

, K. Hada

, M. H. Hecht

, R. Hesper

, L. C. Ho

, P. Ho

, M. Honma

, C. L.

Huang

, L. Huang

, D. H. Hughes

, S. Ikeda

, M. Inoue

, S. Issaoun

, D. J. James

, B. T. Jannuzi

, B. Jeter

, W. Jiang

, M.

D. Johnson

, S. Jorstad

, T. Jung

, M. Karami

, R. Karuppusamy

, T. Kawashima

, G. K. Keating

, M. Kettenis

, J. Kim

, J.

Kim

, J. Kim

, M. Kino

, J. Y. Koay

, P. M. Koch

, S. Koyama

, M. Kramer

, C. Kramer

, T. P. Krichbaum

, C. Kuo

, T. R. Lauer

S. Lee

, Y. Li

, Z. Li

, M. Lindqvist

, R. Lico

, K. Liu

, E. Liuzzo

, W. Lo

, A. P. Lobanov

, L. Loinard

, C. Lonsdale

, R.

, N. R. MacDonald

, J. Mao

, S. Marko

, D. P. Marrone

, A. P. Marscher

, I. Martí-Vidal

, S. Matsushita

, L. D.

Matthews

, L. Medeiros

, K. M. Menten

, Y. Mizuno

, I. Mizuno

, J. M. Moran

, K. Moriyama

, M. Moscibrodzka

, C. Müller

H. Nagai

, N. M. Nagar

, M. Nakamura

, R. Narayan

, G. Narayanan

, R. Neri

, C. Ni

, A. Noutsos

, H. Okino

, H. Olivares

, G.

N. Ortiz-León

, T. Oyama

, F. Özel

, D. C. M. Palumbo

, N. Patel

, U. Pen

, D. W. Pesce

, V. Piétu

, R. Plambeck

, A.

PopStefanija

, B. Prather

, J. A. Preciado-López

, D. Psaltis

, H. Pu

, V. Ramakrishnan

, R. Rao

, M. G. Rawlings

, A. W. Raymond

, L.

Rezzolla

, B. Ripperda

, A. Rogers

, E. Ros

, M. Rose

, A. Roshanineshat

, H. Rottmann

, A. L. Roy

, C. Ruszczyk

, B. R. Ryan

K. L. J. Rygl

, S. Sánchez

, D. Sánchez-Arguelles

, M. Sasada

, T. Savolainen

, F. P. Schloerb

, K. Schuster

, L. Shao

, Z.

Shen

, D. Small

, B. Won Sohn

, J. SooHoo

, F. Tazaki

, P. Tiede

, M. Titus

, K. Toma

, P. Torne

, E. Traianou

, T.

Trent

, S. Trippe

, S. Tsuda

, H. J. van Langevelde

, D. R. van Rossum

, J. Wagner

, J. Wardle

, J. Weintroub

, N. Wex

, R. Wharton

M. Wielgus

, G. N. Wong

, Q. Wu

100

, A. Young

, K. Young

, Z. Younsi

101

, F. Yuan

102

, Y. Yuan

103

, J. A. Zensus

, G. Zhao

, S.

Zhao

, and Z. Zhu

(The Event Horizon Telescope Collaboration)

ffi

liations can be found after the references)

Received 3 September 2019

Accepted 18 October 2019

ABSTRACT

Context.

Realistic synthetic observations of theoretical source models are essential for our understanding of real observational data. In using

synthetic data, one can verify the extent to which source parameters can be recovered and evaluate how various data corruption e

ects can be

calibrated. These studies are the most important when proposing observations of new sources, in the characterization of the capabilities of new or

upgraded instruments, and when verifying model-based theoretical predictions in a direct comparison with observational data.

Aims.

We present the SYnthetic Measurement creator for long Baseline Arrays (

SYMBA

), a novel synthetic data generation pipeline for Very Long

Baseline Interferometry (VLBI) observations.

SYMBA

takes into account several realistic atmospheric, instrumental, and calibration e

ects.

Methods.

We used

SYMBA

to create synthetic observations for the Event Horizon Telescope (EHT), a millimetre VLBI array, which has recently

captured the first image of a black hole shadow. After testing

SYMBA

with simple source and corruption models, we study the importance of

including all corruption and calibration e

ects, compared to the addition of thermal noise only. Using synthetic data based on two example general

relativistic magnetohydrodynamics (GRMHD) model images of M87, we performed case studies to assess the image quality that can be obtained

with the current and future EHT array for di

erent weather conditions.

Results.

Our synthetic observations show that the e

ects of atmospheric and instrumental corruptions on the measured visibilities are significant.

Despite these e

ects, we demonstrate how the overall structure of our GRMHD source models can be recovered robustly with the EHT2017 array

after performing calibration steps, which include fringe fitting, a priori amplitude and network calibration, and self-calibration. With the planned

addition of new stations to the EHT array in the coming years, images could be reconstructed with higher angular resolution and dynamic range.

In our case study, these improvements allowed for a distinction between a thermal and a non-thermal GRMHD model based on salient features in

reconstructed images.

Key words.

galaxies: nuclei – black hole physics – telescopes – atmospheric e

ects – techniques: high angular resolution – techniques: interfer-

ometric

1. Introduction

The giant elliptical galaxy M87 hosts an active galactic nucleus

(AGN) with a radio jet extending to kpc scales (e.g. Owen et al.

These authors contributed equally to this work.

2000). The radio core of M87 shifts inwards with increasing

frequency as the jet becomes optically thin closer to the cen-

tral black hole, resulting in a flat radio spectrum as predicted

by analytical models (Blandford & Königl 1979; Falcke & Bier-

mann 1995). The radio core of M87 coincides with the central

engine at 43 GHz (Hada et al. 2011). At millimetre wavelengths,

Article number, page 1 of 20

arXiv:2004.01161v1 [astro-ph.IM] 2 Apr 2020

&

A proofs:

manuscript no. ms

emission near the event horizon becomes optically thin. Due to

strong gravitational lensing, the black hole is predicted to cast a

‘shadow’ on this emission (Falcke et al. 2000; Dexter et al. 2012;

scibrodzka et al. 2016). The shadow is a region exhibiting an

emission deficit produced by the capture of photons by the event

horizon, with a size enhanced by strong gravitational lensing.

For a Schwarzschild (non-spinning) black hole, the apparent

radius of the black hole shadow is

√

, with

G M

the gravitational radius where

is Newton’s gravitational con-

stant,

is the black hole mass, and

is the speed of light. The

erence in shadow size between a rotating black hole (Kerr

1963) and the Schwarzschild solution is marginal (

.

4%), since

the apparent size is nearly independent of the black hole spin

(Bardeen 1973; Takahashi 2004; Johannsen & Psaltis 2010). Es-

timates for the mass of the supermassive black hole at the cen-

tre of M87 have historically ranged between (3

−

)

from gas-dynamical measurements (Walsh et al. 2013), and

from stellar-dynamical measurements (Geb-

hardt et al. 2011). At a distance of (16

5) Mpc (Bird et al.

2010), the mass measurements correspond to an apparent diam-

eter of the shadow between

∼

as and 42

as.

At 230 GHz, Earth-sized baselines give a nominal reso-

lution of

∼

as, which is certainly su

ffi

cient to resolve

the black hole shadow of M87 for the high-mass estimate.

M87 is therefore one of the prime targets of the Event Hori-

zon Telescope (EHT), the Earth-sized mm-Very Long Base-

line Interferometry (VLBI) array aiming to image a black hole

shadow (Event Horizon Telescope Collaboration et al. 2019b).

The other prime candidate is Sagittarius A* (Sgr A*). With a

better constrained shadow size of

∼

as, this is the black hole

with the largest predicted angular size in the sky. Interstellar scat-

tering e

ects and variability on short time scales (minutes) may

make reconstructing the black hole shadow challenging for this

source. On the other hand, it provides us with opportunities to

study scattering e

ects (Johnson 2016; Dexter et al. 2017; John-

son et al. 2018) and real-time dynamics of the accretion flow

(e.g. Doeleman et al. 2009; Fish et al. 2009; Dexter et al. 2010;

Medeiros et al. 2017; Roelofs et al. 2017; Johnson et al. 2017;

Bouman et al. 2017). In this paper, we focus on synthetic EHT

observations of M87, where orbital timescales are much larger

than those of the observations.

With the EHT data sets and images, it is possible to test gen-

eral relativity in a unique environment (e.g. Bambi & Freese

2009; Johannsen & Psaltis 2010; Psaltis et al. 2015; Goddi

et al. 2017; Event Horizon Telescope Collaboration et al. 2019a).

Also, constraints can be put on models of the accretion flow

around supermassive black holes (e.g. Falcke & Marko

2000;

Yuan et al. 2003; Dexter et al. 2010; 2012; Mo

scibrodzka et al.

2014; 2016; Chan et al. 2015; Broderick et al. 2016; Gold et al.

2017; Event Horizon Telescope Collaboration et al. 2019e).

In 2017, the EHT consisted of the IRAM 30-metre tele-

scope on Pico Veleta in Spain, the Large Millimeter Tele-

scope (LMT) in Mexico, the Atacama Large Millemeter Ar-

ray (ALMA), the Atacama Pathfinder Experiment (APEX) tele-

scope in Chile, the Sub-Millimeter Telescope (SMT) in Arizona,

the Sub-Millimeter Array and James Clerk Maxwell Telescope

(JCMT) in Hawaii, and the South Pole Telescope (SPT). In the

April 2017 observing run (hereafter EHT2017) and a subsequent

two-year analysis period, the EHT imaged the M87 black hole

shadow within a 42

as asymmetric emission ring (Event

Horizon Telescope Collaboration et al. 2019d;f). The measured

ring size, when associated with a black hole shadow, leads to an

angular size of one gravitational radius of 3

as (Event

Horizon Telescope Collaboration et al. 2019f). At the adopted

distance of 16

−

Mpc that was calculated from multiple mea-

surements (Bird et al. 2010; Blakeslee et al. 2009; Cantiello et al.

2018), this angular size corresponds to a black hole mass of

stat

sys

)

, which is consistent with the

stellar-dynamical mass measurement by Gebhardt et al. (2011).

Over the years, synthetic data have proven to be of impor-

tance for demonstrating the capabilities of the EHT. They were

also essential for developing new strategies to increase the sci-

entific output of the rich, yet challenging, observations.

Doeleman et al. (2009) and Fish et al. (2009) used the As-

tronomical Image Processing System (AIPS)

task UVCON to

calculate model visibilities for the EHT array, showing that sig-

natures of source variability could be detected in Sgr A* by using

interferometric closure quantities and polarimetric ratios. The

MIT Array Performance Simulator (MAPS)

was used in sev-

eral EHT synthetic imaging studies. Lu et al. (2014) used it to

test the ability of the EHT to reconstruct images of the black

hole shadow for several models of the accretion flow of M87.

Fish et al. (2014) demonstrate that for Sgr A*, the blurring ef-

fect of interstellar scattering could be mitigated if the proper-

ties of the scattering kernel are known. Lu et al. (2016) showed

that source variability could also be mitigated by observing the

source for multiple epochs and applying visibility averaging,

normalization, and smoothing to reconstruct an image of the av-

erage source structure.

Typically, the only data corruption included in these syn-

thetic data sets is thermal noise, although Fish et al. (2009)

also included instrumental polarization. More corruptions can be

added with the

eht-imaging

library

. Chael et al. (2016; 2018)

simulated polarimetric EHT images of Sgr A* and M87, and

included randomly varying complex station gains and elevation-

dependent atmospheric opacity terms. With the stochastic optics

module in

eht-imaging

, the input model images can be scat-

tered using a variable refractive scattering screen, and the scat-

tering can be mitigated by solving for the scattering screen and

image simultaneously (Johnson 2016). However, scattering ef-

fects are only relevant for observations of Sgr A*.

eht-imaging

can also simulate observations following a real observing sched-

ule, and copy the

-coverage and thermal noise directly from

existing data sets. It also includes polarimetric leakage corrup-

tions (Event Horizon Telescope Collaboration et al. 2019d).

Despite these recent advances in synthetic data genera-

tion, there are still di

erences between synthetic and real mm-

VLBI data sets. So far, synthetic EHT data sets have not been

frequency-resolved, and gain o

sets have only been included as

random relative o

sets drawn from a Gaussian with a fixed stan-

dard deviation, rather than being based on a physical model.

Moreover, no calibration e

ects are taken into account in the

synthetic data products. It is essentially assumed that residual de-

lays, phase decoherence due to atmospheric turbulence, and sig-

nal attenuation caused by the atmospheric opacity are perfectly

calibrated. In

eht-imaging

, atmospheric turbulence can be in-

cluded by fully randomizing the phases (with the option of fixing

them within a scan). In real mm-VLBI data, atmospheric tur-

bulence results in rapid phase wraps. The correlated phases are

not fully randomized, but evolve continuously over frequency

and time, allowing to perform fringe fitting and average com-

plex visibilities coherently on time scales set by the atmospheric

coherence time.

http://www.aips.nrao.edu

https://www.haystack.mit.edu/ast/arrays/maps

https://github.com/achael/eht-imaging

Article number, page 2 of 20

F. Roelofs, M. Janssen et al.: SYMBA: An end-to-end VLBI synthetic data generation pipeline

In this paper, we present the SYnthetic Measurement creator

for long Baseline Arrays (

SYMBA

) – a new synthetic VLBI data

generation and calibration pipeline.

We generate raw synthetic data with

MeqSilhouette

(Blecher et al. 2017; Natarajan et al. 2019), which includes a

tropospheric module and physically motivated antenna point-

ing o

sets (Section 2). We then calibrate the raw data using the

new

CASA

(McMullin et al. 2007) VLBI data calibration pipeline

rPICARD

(Janssen et al. 2019b), applying a fringe fit and a priori

amplitude calibration (Section 3). The overall computing work-

flow of

SYMBA

is outlined in Section 4. We describe our simu-

lated observational setup (antenna and weather parameters and

observing schedule) in Section 5 and our input source models

for the synthetic data generation in Section 6. In Section 7, we

demonstrate the e

ects of simulated data corruptions and subse-

quent calibration. We illustrate the capabilities of

SYMBA

in Sec-

tion 8 based on three scientific case studies. In these studies we

show 1) how well we can distinguish between two example gen-

eral relativistic magnetohydrodynamics (GRMHD) models with

erent descriptions for the electron temperatures with the cur-

rent and future EHT array, 2) how the EHT would perform un-

der di

erent weather conditions, and 3) how pre-2017 models of

M87 compare to the observed image of the black hole shadow.

In Section 9, we summarize our conclusions and discuss future

work.

2. Synthetic data generation with

MeqSilhouette

MeqSilhouette

(Blecher et al. 2017; Natarajan et al. 2019)

is a synthetic data generator designed to simulate high fre-

quency VLBI observations. While visibilities of real radio in-

terferometric observations are produced by correlating recorded

voltage streams from pairs of telescopes,

MeqSilhouette

pre-

dicts visibilities directly from the Fourier Transform of an in-

put sky model. For simple ASCII input models (e.g. a set of

Gaussian components, each with an independent spectral in-

dex),

MeqTrees

(Noordam & Smirnov 2010) is used for the

visibility prediction. FITS-based

sky models are converted

with the

wsclean

ringa et al. 2014) algorithm. The sig-

nal path is described by the Measurement Equation formal-

ism (Hamaker et al. 1996), breaking down the various e

ects

on the visibilities into a chain of complex 2

2 Jones matri-

ces (Jones 1941; Smirnov 2011a;b;c).

MeqSilhouette

gener-

ates frequency-resolved visibilities, with a bandwidth and num-

ber of channels set by the user. Frequency-resolved visibili-

ties are required for the calibration of signal path variations in-

troduced by the troposphere. In particular, synthetic data from

MeqSilhouette

has been used to validate the CASA-based data

reduction path of the EHT. Moreover, channelized data allows

for the introduction of frequency dependent leakage of polarized

signals at telescopes’ receivers, the inclusion of wavelength de-

pendent Faraday rotation and spectral indices in source models,

and multi-frequency aperture synthesis, which can yield signifi-

cant improvements to the

-coverage.

It is also possible to gen-

erate corrupted data sets from time-dependent polarized emis-

sion models in full Stokes and to follow an observed schedule

https://bitbucket.org/M_Janssen/symba

https://github.com/rdeane/MeqSilhouette_public_v0.1

https://bitbucket.org/M_Janssen/picard

See

https://fits.gsfc.nasa.gov/fits_documentation.

html

for a definition of the FITS standard.

For example, the EHT is currently able to observe with two sidebands

separated by 18 GHz, which

MeqSilhouette

can replicate.

from a VEX file.

A key design driver of

MeqSilhouette

to generate synthetic data (and associated meta-data) in a for-

mat that is seamlessly ingested by the

CASA

software package.

The native format is the MeasurementSet (MS)

, but the visi-

bilities can also be exported to UVFITS.

We briefly describe

the added tropospheric and instrumental corruptions below, re-

ferring to Blecher et al. (2017) and Natarajan et al. (2019) for

more details.

2.1. Tropospheric corruptions

The e

ects of the troposphere on the measured visibilities can be

separated into those resulting from a mean atmospheric profile,

and those resulting from atmospheric turbulence.

2.1.1. Mean troposphere

The mean troposphere causes time delays, resulting in phase

slopes versus frequency and an attenuation of the visibility am-

plitudes due to absorption of radiation in molecular transitions

(Thompson et al. 2017). In the mm-wave regime, absorption

lines are mostly caused by rotational transitions of H

O and

. Apart from the individual lines, there is a general increase

of the opacity with frequency due to the cumulative e

ect of

pressure-broadened H

O lines peaking in the THz-regime (Car-

illi & Holdaway 1999).

MeqSilhouette

calculates the attenuation and time delays

using the Atmospheric Transmission at Microwaves (

ATM

) soft-

ware (Pardo et al. 2001). It integrates the radiative transfer equa-

tion

(

)

(

)

−

(

)

(

)

(1)

where

(

) is the specific intensity at frequency

at path length

coordinate

, and

and

are the emission and absorption co-

ffi

cients, respectively. In thermodynamic equilibrium, the latter

are related through Kirchho

’s law,

(

)

(2)

where

(

) is the Planck spectrum at temperature

. In order to

integrate Equation 1,

ATM

calculates

as a function of altitude.

For a specific transition,

is proportional to the photon energy,

the transition probability (Einstein coe

ffi

cient), molecular densi-

ties of the lower and upper states, and the line shape including

pressure and Doppler broadening.

is related to the refractive

index of the medium via the Kramers-Kronig relations. The in-

troduced time delay is then calculated from the refractive index.

As is evident from Kirchho

’s law (Equation 2), the atmo-

sphere not only absorbs, but also emits radiation. This process

leads to an increase in system temperature (sky noise), which

also follows from the integration in

ATM

and is included in the

noise budget with an elevation and therefore time-dependent

contribution.

See

https://vlbi.org/vlbi-standards/vex/

for a definition

of the VEX file format.

See

https://casa.nrao.edu/Memos/229.html

for the definition

of the MeasurementSet format.

See

ftp://ftp.aoc.nrao.edu/pub/software/aips/TEXT/

PUBL/AIPSMEM117.PS

for a description of the UVFITS file data

format.

Article number, page 3 of 20

&

A proofs:

manuscript no. ms

2.1.2. Turbulent troposphere

Apart from the mean troposphere induced amplitude attenuation,

signal delay, and sky noise, a major source of data corruptions

in the mm regime is tropospheric turbulence. Rapid evolution of

the spatial distribution of tropospheric water vapour causes the

signal path delay to vary on short (

∼

10 s) time scales. This then

leads to rapid and unpredictable rotations of the visibility phase,

posing challenges for fringe fitting. Because of atmospheric tur-

bulence, uncalibrated visibilities can not be coherently averaged

beyond the atmospheric coherence time. Absolute phase infor-

mation can only be recovered with phase-referencing (Beasley

& Conway 1995). For imaging mm-VLBI data, one often needs

to rely on closure phases (e.g. Chael et al. 2018). Closure phase

is the sum of visibility phases on a triangle of baselines, in which

many station-based instrumental and atmospheric corruptions

cancel out (Jennison 1958; Rogers et al. 1974).

MeqSilhouette

, turbulent phase errors are added to the

visibilities assuming that the atmospheric turbulence can be rep-

resented by a thin phase-changing scattering screen. Similar

to simulations of interstellar scattering (e.g. Johnson & Gwinn

2015), the turbulent substructure of the screen is assumed to be

constant in time while the screen itself is moving with a constant

transverse velocity

. The screen velocity sets the atmospheric

coherence time together with the spatial phase turbulence scale

on the screen. The introduced phase o

sets are described by a

phase structure function that takes a power law form,

(

x

x

′

)

〈

[

(

x

x

′

)

−

(

x

)]

〉≈

(

)

(3)

where

x

and

x

′

are spatial coordinates on the screen,

(

x

−

x

′

)

is the phase coherence length such that

(

)

rad,

csc (elevation) is the airmass towards the horizon

, and

3 if one assumes Kolmogorov turbulence, which is sup-

ported by Carilli & Holdaway (1999). The nature of the scatter-

ing is set by the ratio of

and the Fresnel scale

√

where

is the distance between the observer and the scatter-

ing screen. With

measured to be

∼

−

700 m (Masson 1994;

Radford & Holdaway 1998) and a water vapour scale height of

2 km, we have

≈

45 m

and are in the weak scatter-

ing regime. This means that most of the received power on the

ground originates from a screen area

weak

≈

, rather than

from disjoint patches, as is the case for interstellar scattering.

At a distance of 2 km, 1 mas corresponds to

∼

m, and the

Field of View (FoV) of the array is much smaller than

. The

phase error is therefore assumed to be constant across the FoV,

and the structure function can be written as

(

)

(

)

where

is the bulk transverse velocity of the phase screen. From

this, a phase error time sequence can be computed directly. Due

to the long baselines, atmospheric corruptions can be modelled

independently at each station (Carilli & Holdaway 1999). For

a given coherence time

(Treuhaft & Lanyi 1987) at

a reference frequency

, Blecher et al. (2017) showed that the

temporal variance of the phase for a power-law turbulence as a

function of frequency

can be modelled as

(

,ν

)

[

] (

int

)

(

)

rad

(4)

The csc (elevation) dependence of the airmass is an approximation

assuming a planar rather than a spherical atmosphere, which breaks

down at elevations below

∼

10 degrees (Paine 2019). For the synthetic

observations in this work, we set the elevation limit to 10 degrees as

is typically done for real VLBI observations. Hence, the csc (elevation)

approximation has a negligible e

ect on our results.

where

int

is the data integration time and

is taken as the low-

est frequency in the data.

MeqSilhouette

uses Equation 4 to

compute the tropospheric phase turbulence using

3. A

constant amount of precipitable water vapour at zenith (PVW

)

is assumed, mixed evenly into the atmosphere. An increase in

the phase variance due to the PWV therefore enters through the

amount of airmass towards the horizon in Equation 4. The spec-

ified coherence time

(PWV

) should decrease with in-

creasing precipitable water vapour content in the atmosphere,

although other factors such as wind speed also a

ect

. No sud-

den phase jumps due to inhomogeneities in the atmosphere (e.g.

clouds or airmass boundary kinks) along the line of sight are

simulated. Phase turbulence and resulting decorrelation within

an integration time

int

are not simulated by

MeqSilhouette

For realistic results,

int

should therefore preferably be set to well

within

, as is the case for real observations. Delay-related deco-

herence e

ects within individual frequency channels are also not

simulated. It is assumed that frequency resolution is su

ffi

ciently

high to make this e

ect negligible, as it is done in modern corre-

lators.

2.2. Receiver noise

The System Equivalent Flux Density (SEFD) of a station is

a measure for its overall noise contribution.

MeqSilhouette

reads

, the contribution from the receiver noise to the SEFD,

from input files. Receiver temperatures

are typically deter-

mined from real data by extrapolating system temperatures to

zero airmass and the receiver noise contribution in units of Jan-

sky (Jy) follows as

DPFU

(5)

Here, the DPFU is the telescope’s ‘degree per flux unit’ gain, de-

fined as DPFU

dish

(

)

, with

the aperture e

ffi

ciency

(taken to be constant during observations),

dish

the geometric

area of the dish, and

the Boltzmann constant.

2.3. The full noise budget

Visibilities on all baselines are corrupted by the addition of noise

as a complex Gaussian variable with standard deviation

√

SEFD

∆

int

(6)

where SEFD

is the system equivalent flux density from sta-

tion

with combined contributions from the atmosphere and

receiver,

∆

is the channel bandwidth,

int

is the correlator in-

tegration time, and

is a quantization e

ffi

ciency factor, set to

0.88 for standard 2-bit quantization. We assume perfect quan-

tization thresholds when simulating the cross-correlation data.

Therefore, we do not need to simulate the auto-correlations to

correct for erroneous sampler thresholds. All noise sources along

the signal chain (sky noise, turbulence, and thermal noise from

the instrument) enter into

MeqSilhouette

produces visi-

bilities in a circular polarization basis, that is LL, RR, LR, and

RL. The noise on, for example, the Stokes I data is a factor

√

smaller.

2.4. Antenna pointing errors

Pointing o

sets of individual antennas manifest as a time and

station dependent amplitude error. They cause a drop of the vis-

ibility amplitudes

on a

baseline as the maximum of the

Article number, page 4 of 20

F. Roelofs, M. Janssen et al.: SYMBA: An end-to-end VLBI synthetic data generation pipeline

antenna primary beam is not pointed on the source. The primary

beam profile of a station

is modelled as a Gaussian with a

full width at half maximum

FWHM

, which is related to the

Gaussian’s standard deviation by a factor of 2

√

2 ln 2

≈

35. A

Gaussian beam is justified since the pointing o

sets are not large

enough that a Gaussian and Bessel function deviate (i.e. near

the first null), see Middelberg et al. (2013). No further system-

atic point e

ects, such as refraction, are considered here. Point-

ing o

sets

are drawn from a normal distribution

centred

around zero, with a standard deviation given by a specified rms

pointing o

set

rms

. The resulting visibility amplitude loss

∆

exp







−

8 ln







FWHM













(7)

(

,σ

rms

)

describes a data corruption e

ect caused by an erroneous source

tracking of the telescopes.

SYMBA

, we employ two types of pointing o

sets, which

occur on short and long timescales, respectively. The short

timescale variations are caused by the atmospheric seeing and

wind shaking the telescope, resulting in a displacement of the

sky source with respect to an otherwise perfectly pointed tele-

scope beam. Here,

SYMBA

draws values of

from

rms

timescales set by the atmospheric coherence time. The long

timescale variations are caused by sub-optimal pointing solu-

tions adopted by a telescope.

SYMBA

simulates these by adopting

a new value of

every

∼

5 scans and letting these pointing

sets deteriorate by

∼

1 in every scan until a new o

set

is determined. For simplicity, the

are drawn from the same

rms

, multiplied by a factor

∼

5. For a scan number

, the

ect of an incorrect pointing model is thus given as

(

)

mod

(

,σ

rms

)

(8)

2.5. Leakage and gain errors

Complex gain errors

, that would translate to errors in

the DPFUs and phase gains in real observations, and com-

plex leakage e

ects (

-terms) can be added as well. For ob-

served

corrupted (obs) visibilities from a baseline of stations

and

-terms cause artificial instrumental polarization as a ro-

tation of the cross-hand visibilities in the complex plane by twice

the station’s feed rotation angles

(Conway & Kronberg 1969):

obs

true

[

(

)

∗

]

(9)

obs

true

[

−

(

)

∗

−

]

(10)

Here,

are the leakage terms, with a superscript indicating the

polarization, and

√

−

1. The star denotes complex conjuga-

tion. More complex and realistic polarimetric e

ects are avail-

able in the forthcoming release of

MeqSilhouette v2

(Natara-

jan et al. 2019).

3. Synthetic data calibration with

rPICARD

The goal of

SYMBA

is to create synthetic observations which

match real data as closely as possible. After the simulation

of physically motivated data corruptions by

MeqSilhouette

the synthetic data are passed through the

rPICARD

calibration

pipeline (Janssen et al. 2019b). The data are treated in the same

way as actual correlated visibilities and a model-agnostic cali-

bration (Smirnov 2011a) of phases and amplitudes is performed

based on information typically available for real observations.

The

atmospheric

signal

attenuation

introduced

MeqSilhouette

is corrected by recording opacity values

for each station at the start of each scan. This is the equivalent

of measuring opacity-corrected system temperatures with a hot-

load calibration scan in real VLBI observations (Ulich & Haas

1976), which leaves intra-scan opacity variations unaccounted

for. As

MeqSilhouette

does not simulate the digitization

when radio telescopes record data, nor the correlation process,

the simulated visibilities are already scaled to units of flux

density, as derived from the input source model. Therefore,

unity amplitude gains are used and the system temperatures are

set to exp (

) for the amplitude calibration, with

describing the

atmospheric opacity (see Sect. 4.2 in Janssen et al. (2019b)).

Amplitude losses due to pointing o

sets can not be corrected

with this standard VLBI amplitude calibration method.

The phases are calibrated with the

CASA

Schwab-Cotton

(Schwab & Cotton 1983) fringe fitter implementation. With this

method, station gains for phases, rates, and delays are solved

with respect to a chosen reference station.

rPICARD

uses a pri-

oritized list of reference stations (based on availability). For the

EHT, these are ALMA

→

LMT

→

APEX

→

SMT

→

PV. All

solutions are re-referenced to a single common station in the

end. Optimal fringe fit solution intervals are found based on the

signal-to-noise ratio (S

N) of the data in each scan. The search

intervals range from twice the data integration time (typically

∼

0.5-1 s) to 60 s. Within this interval, the smallest timescale

which yields fringe detections with S

5.5 on all baselines

for which the source can be detected, is chosen (Janssen et al.

2019a). Figure 1 shows estimated S

N values for a range of

fringe fit solution intervals and di

erent simulated coherence

times. The presence of (frequency independent) atmospheric de-

lays and absence of instrumental delays in the synthetic data war-

rants a combined fringe fit solution over the whole frequency

band for a maximum S

N. Usually,

rPICARD

would smooth

Integration time [seconds]

S/N

sqrt S/N increase

= 15

= 2

= 1

Fig. 1:

N estimates for

rPICARD

fringe solutions. The plot-

ted points indicate the estimated average FFT S

N values by the

CASA

fringefit

code for di

erent integration times (solution inter-

vals), segmenting a 15 minute long scan of a

MeqSilhouette

observation of a 4 Jy point source on the ALMA-APEX base-

line. Di

erent symbols correspond to di

erent coherence times

(Equation 4) used for the simulation of atmospheric turbulence.

The dashed line shows the expected increase in S

N for an infi-

nite coherence time without added noise corruptions.

Article number, page 5 of 20

&

A proofs:

manuscript no. ms

Fig. 2:

Delay between ALMA and LMT. The delay is solved a

function of time by the fringe fitting calibration step. The input

source model is a 4 Jy point source.

solved delays within scans to remove potential outliers. This

is done under the assumption that an a priori delay model like

Calc

Solve

has been applied at the correlation stage, which al-

ready takes out the bulk of the delay o

sets. For the synthetic

data generation, no atmospheric delay model is applied and

rPICARD

has to solve for steep residual delay gradients caused

by the wet and dry atmospheric components within scans (Fig-

ure 2). Smoothing of solved delays is therefore disabled here.

The last step of the calibration pipeline is the application of

the amplitude and phase calibration tables, and averaging of the

data in frequency within each spectral window. The calibrated

and averaged data are then exported in the UVFITS file format.

Optionally, an additional UVFITS file can be provided as input.

SYMBA

then uses

eht-imaging

to reproduce the

-coverage

from that file. For a UVFITS file from a real observation, this

means taking into account time periods where telescopes drop

out of the observed schedule and all non-detections. Thereby,

a comparison of synthetic and real data is una

ected by

coverage.

Finally, the synthetic UVFITS data are averaged in 10 second

intervals and a ‘network calibration’ (Fish et al. 2011; Johnson

& Gwinn 2015; Blackburn et al. 2019; Event Horizon Telescope

Collaboration et al. 2019c) is performed with the

eht-imaging

software. The gains of non-isolated (redundant) stations, which

have a very short baseline to another nearby station can be cal-

ibrated if the model of the observed source is known at large

scales. For the 2017 EHT observations, ALMA was able to pro-

vide accurate large-scale source models, allowing for a network

calibration of the co-located ALMA

APEX and SMA

JCMT

sites. For our synthetic observations, we use the known total flux

density of the input model.

4. Computing workflow

SYMBA

is controlled by a single input ASCII file. The observed

schedule can either follow a VEX file or explicitly set start time,

duration, number of scans, and gaps between scans. If the VEX

file has been used for a real observation, a UVFITS file can be

provided to match the

-coverage. All antenna and weather pa-

rameters are also set in ASCII files. The input source model can

be provided as FITS or ASCII file, as a single model or multi-

ple frames from a time-variable source, and contain only Stokes

http://astrogeo.org/psolve/

I or full polarization information. The input model is Fourier

Transformed and corrupted by

MeqSilhouette

. The resultant

visibilities are calibrated by

rPICARD

, and optionally network

calibrated and imaged by

eht-imaging

SYMBA

outputs a FITS

file of the final reconstructed source model, the calibrated and

self-calibrated visibilities in UVFITS and ASCII format, and di-

agnostic plots of the calibration process. The pipeline is fully

dockerized.

. An overview of the workflow is shown in Fig-

ure 3.

5. Simulated observation setup

SYMBA

is able to create synthetic observations for any VLBI ar-

ray. Here, we outline the antenna and weather parameters and ob-

serving schedules adopted for the creation of our synthetic data

sets.

5.1. EHT2017 array

Our primary array consists of the 2017 EHT stations, exclud-

ing the SPT station for which M87 is always below the horizon.

The antenna parameters are summarized in Table 1. The receiver

SEFDs of the primary array have been estimated by extrapolat-

ing system temperature measurements to zero airmass, follow-

ing Janssen et al. (2019b). Full width at half maximum 230 GHz

beam sizes (

FWHM

) and dish diameters (

) were taken from

the websites and documentation for each individual site. Point-

ing rms o

sets (

rms

) have been based on a priori station in-

formation and typical inter- and intra-scan amplitude variations

seen in EHT data. All o

sets lie within o

ffi

cial telescope speci-

fications. Aperture e

ffi

ciencies (

) were estimated with

∼

10%

accuracy from planet observations (Janssen et al. 2019a; Event

Horizon Telescope Collaboration et al. 2019c). In our synthetic

observations, we have added gain errors (

err

) listed in Table 1

in accordance with these uncertainties. Additionally, a polariza-

tion leakage corruption has been added at a

5% level for all

stations. This corruption has been left uncalibrated by

rPICARD

to mimic the current capabilities of the EHT, which did not per-

form a polarization calibration for the first scientific data release

(Event Horizon Telescope Collaboration et al. 2019c).

The weather parameters are summarized in Table 2. For

the ground temperature

, pressure

, and precipitable water

vapour PWV, we used the median values measured during the

EHT2017 campaign (5-11 April) at the individual primary sites,

logged by the VLBI monitor (Event Horizon Telescope Collabo-

ration et al. 2019b). No weather information was available from

the VLBI monitor for ALMA. We adopted the values measured

at the nearby station APEX.

The radiometers at the sites measure the atmospheric opacity

, while

MeqSilhouette

takes the PWV as input. The 225 GHz

opacity can be converted to PWV in mm using

PWV

−

dry

−

air

(11)

where

dry

−

air

is the dry air opacity and the slope

is in

mmH

−

and

dry

−

air

have been measured at some sites

and both tend to decrease with site altitude, but the errors on

these measurements are not well known (Thompson et al. 2017;

Thomas-Osip et al. 2007; and references therein): the calibration

needs an accurate independent measure of the water vapour

column density at the same site as the radiometer, which is only

https://www.docker.com/

Article number, page 6 of 20

F. Roelofs, M. Janssen et al.: SYMBA: An end-to-end VLBI synthetic data generation pipeline

Source

model

Antenna

information

Vex

schedule

Real

observation

Master

input file

Observation with

MeqSilhouette

Corruption

and calibration

information

Generate

ANTAB

ANTAB table

Corrupted

observation

Calibration

with rPICARD

Calibrated

observation

Flag unob-

served scans

Flagged

calibrated

observation

Network cal

with eht-imaging

Network

calibrated

observation

Imaging

with eht-

imaging

Reconstructed

image

Frequency

setup,

Observation

schedule,

Source infor-

mation,

Requested

corruptions

Real observation

Total flux, gain tolerance

Field of view,

Gain tolerance

Reference antennas, fringefit search range

Fig. 3:

Computing workflow flowchart of

SYMBA

. Red borders and arrows indicate the main data path. Dashed borders and arrows

indicate optional steps that may be skipped (for example, imaging could be done without network calibration). Yellow boxes are

auxiliary input files; the master input file is indicated by the red box. Green ellipses are actions, and blue boxes are data products.

Text next to arrows lists the information from the master input file that is used for a specific action.

Table 1:

Antenna parameters adopted in our synthetic observations.

Year

Antenna

X (m)

Y (m)

Z (m)

(m)

(Jy)

err

D P

rms

(")

FWHM

(")

2017

ALMA

2225061

-5440057

-2481681

0.73

1.02

0.05

1.0

APEX

2225040

-5441198

-2479303

0.63

3300

0.97

0.05

1.0

JCMT

-5464585

-2493001

2150654

0.52

6500

1.05

0.05

1.0

LMT

-768716

-5988507

2063355

0.31

2400

0.85

0.05

1.0

5088968

-301681

3825012

0.43

1000

1.03

0.05

0.5

SMA

-5464555

-2492928

2150797

0.73

3300

0.96

0.05

1.5

SMT

-1828796

-5054407

3427865

0.57

7700

0.93

0.05

1.0

2018

GLT

541647

-1388536

6180829

0.63

3300

1.08

0.05

1.0

2020

-1994314

-5037909

3357619

0.63

3300

0.96

0.05

1.0

PDB

4523951

468037

4460264

0.52

750

0.95

0.05

1.0

2020

AMT

5627890

1637767

-2512493

0.52

1990

1.03

0.05

1.0

Table 2:

Weather parameters adopted in our synthetic observa-

tions.

Antenna

PWV (mm)

(mb)

(K)

(s)

ALMA

1.5

555

271

APEX

1.5

555

271

JCMT

1.5

626

278

LMT

5.7

604

275

2.9

723

270

SMA

1.5

626

278

SMT

4.4

695

276

GLT

1.7

1000

254

2.5

793

282

PDB

3.0

747

270

AMT

6.2

772

287

available for a few EHT sites. Also,

dry

−

air

is generally small

(order 10

−

), making it challenging to measure.

For these reasons, climatological modelling likely provides

better estimates than empirical measurements here. To estimate

and

dry

−

air

, we use the Modern-Era Retrospective Analysis

for Research and Applications, version 2 (

MERRA-2

) from the

NASA Goddard Earth Sciences Data and Information Services

Center (GES DISC) (Gelaro et al. 2017). In a reanalysis model

MERRA-2

, variables such as the air temperature and mix-

ing ratios of di

erent molecules are computed based on ground-

and space-based measurements. They depend on time, atmo-

spheric pressure level, and latitude and longitude coordinates.

We use 2006-2016

MERRA-2

data averaged over seasons (per

three months) and latitude zones (antarctic and arctic, south-

ern and northern mid-latitudes, and tropical)

. For each pressure

layer and latitude zone, we then perform radiative transfer at 225

GHz with the

am

atmospheric model software (Paine 2019) with

and without water vapour included to calculate

and

dry

−

air

the March-April-May season (which is the usual EHT observ-

As available on

https://www.cfa.harvard.edu/~spaine/am/

cookbook/unix/zonal/

Article number, page 7 of 20

&

A proofs:

manuscript no. ms

ing season). We then interpolate these to the pressure level of

each EHT site and calculate the PWV from the measured

us-

ing equation 11.

Atmospheric coherence times

were estimated based on the

characteristics of the 2017 EHT measurements for the primary

array. Precise station-based coherence times are di

ffi

cult to ob-

tain and will vary from day to day due to changes in the weather

conditions. For this paper, estimates are taken that match well to

decent to poor weather. The values are summarized in Table 2.

A larger parameter space will be studied in future work to char-

acterize the e

ect of varying weather conditions.

5.2. Enhanced EHT array

Apart from simulated observations with the stations that joined

the 2017 EHT campaign, we also simulate observations with

an enhanced EHT array including four additional stations. The

Greenland Telescope (GLT, Ra

ffi

n et al. 2014) is currently lo-

cated at Thule air base (it will be relocated to Summit Station

near the peak of the Greenland ice sheet) and joined the EHT in

2018. The 12-m telescope on Kitt Peak (KP, Freund et al. 2014)

in Arizona and the IRAM NOEMA interferometer on Plateau

de Bure (PDB, Guilloteau et al. 1992) in France were to join

in the cancelled 2020 observations and will join in future cam-

paigns. Finally, the Africa Millimetre Telescope (AMT, Backes

et al. 2016), is planned to be built on the Gamsberg in Namibia.

For these sites, we estimated weather parameters using the

MERRA-2

inst3_3d_asm_Np data product, which has a time reso-

lution of 3 hours, and is distributed on a grid having 0.625 degree

longitude by 0.5 degree latitude with 42 vertical pressure levels

between 0.1 and 1000 mbar. From this dataset, we took the 25th

percentile (representing good weather) of the air temperature and

specific humidity measured on 11 April in the last two decades

(1999-2018).

At each pressure level, these quantities were then

linearly interpolated between the four grid points nearest to the

observatory site. We then performed an integration of the humid-

ity over the pressure levels using the

am

atmospheric model soft-

ware (Paine 2019) to obtain the total PWV above the site. The

starting point for the integration over the pressure levels was de-

termined by interpolating the geopotential height (pressure as a

function of altitude) to the altitude of the site. The geopotential

height data were downloaded through NASA’s Giovanni portal.

The resulting weather parameters are listed in Table 2. The GLT

site is close to sea level, but the closest

MERRA-2

grid points are

further inland at higher altitudes. Hence, the air temperature and

specific humidity were extrapolated from a pressure level of 925

mbar to the GLT site pressure level of 1000 mbar before the in-

tegration was done in

am

The receiver temperatures and aperture e

ffi

ciencies for the

future stations were estimated from existing stations. The GLT

and KP antennas are ALMA prototypes like APEX, so the values

for APEX were adopted here. The NOEMA interferometer has

ten 15-metre dishes, so the sensitivity was scaled accordingly

from the JCMT, including a phasing e

ffi

ciency of 87%. The cur-

rently envisioned dish for the AMT is the now defunct Swedish-

ESO Submillimetre Telescope (SEST, Booth et al. 1989) tele-

scope in Chile. With a sideband separating receiver, the current

estimate for the SEFD of the AMT is 1990 Jy (A. Young, priv.

comm.).

It should be noted that the current EHT observing strategy

is to trigger a few observing days in a March

April observ-

ing window, based on optimal weather conditions across all sites

(Event Horizon Telescope Collaboration et al. 2019b).

10.0

7.5

5.0

2.5

0.0

2.5

5.0

7.5

10.0

)

EHT2017

GLT

PDB

AMT

Intra-new

Fig. 4:

-coverage towards M87. Di

erent colors show base-

lines within the EHT2017 array, baselines between the EHT2017

array and four (potential) new stations, and baselines between

the new stations (labelled as ‘Intra-new’).

Hereafter, the EHT2017 array plus GLT, KP, and PDB are

referred to as EHT2020. When the AMT is also included, it is

referred to as EHT2020

AMT.

5.3.

-coverage

Figure 4 shows the

-coverage towards M87 for the EHT2017

array and expansions with future stations. The EHT2017 sched-

ule for 11 April was adopted. To accommodate the eastward ex-

pansion of the array with the AMT and PDB, ten-minute scans

were prepended to the schedule at 30-minute intervals starting

when the source is at an elevation of more than ten degrees at

both the AMT and PDB. The GLT, strategically located between

the European and American mainland, adds north-south base-

lines to all stations, significantly increasing the north-south res-

olution due to long baselines to ALMA

APEX. KP and PDB

add short baselines to the SMT and PV, respectively, filling the

-gaps between the intrasite baselines and the SMT-LMT base-

line. These gaps on short

-spacings pose challenges for image

reconstruction with the EHT2017 array (Event Horizon Tele-

scope Collaboration et al. 2019d). Finally, the AMT adds north-

south baselines to the European stations, east-west baselines to

ALMA

APEX, and increases the north-east to south-west reso-

lution by adding baselines to the LMT and SMT

KP. The AMT

has a larger impact for observations of more southern sources

like Sgr A*. Unless noted otherwise, all synthetic data sets in

this work are generated based on the 11 April observing sched-

ule for a source in the direction of M87 for the EHT2017 and

EHT2020 arrays, and the extended schedule described above is

used for EHT2020

AMT array.

Article number, page 8 of 20