SYMBA: An end-to-end VLBI synthetic data generation pipeline - Simulating Event Horizon Telescope observations of M 87

A&A 636, A5 (2019)

https:

doi.org

10.1051

0004-6361

201936622

ESO 2020

Astronomy

Astrophysics

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

Simulating Event Horizon Telescope observations of M 87

The Event Horizon Telescope Collaboration: F. Roelofs

, M. Janssen

, I. Natarajan

, R. Deane

3 ,2

, J. Davelaar

, H. Olivares

, O. Porth

5 ,4

S. N. Paine

, K. L. Bouman

7 ,6 ,8

, R. P. J. Tilanus

1 ,9 ,10

, I. M. van Bemmel

, H. Falcke

, K. Akiyama

12 ,13 ,14 ,7

, A. Alberdi

, W. Alef

K. Asada

, R. Azulay

18 ,19 ,16

, A. Baczko

, D. Ball

, M. Balokovi

7 ,6

, J. Barrett

, D. Bintley

, L. Blackburn

7 ,6

, W. Boland

G. C. Bower

, M. Bremer

, C. D. Brinkerink

, R. Brissenden

7 ,6

, S. Britzen

, A. E. Broderick

25 ,26 ,27

, D. Broguiere

, T. Bronzwaer

D. Byun

28 ,29

, J. E. Carlstrom

30 ,31 ,32 ,33

, A. Chael

,??

, C. Chan

20 ,35

, S. Chatterjee

, K. Chatterjee

, M. Chen

, Y. Chen (

陈

永

军

)

37 ,38

I. Cho

28 ,29

, P. Christian

20 ,6

, J. E. Conway

, J. M. Cordes

, G. B. Crew

, Y. Cui

40 ,41

, M. De Laurentis

42 ,4 ,43

, J. Dempsey

G. Desvignes

16 ,79

, J. Dexter

, S. S. Doeleman

7 ,6

, R. P. Eatough

, V. L. Fish

, E. Fomalont

, R. Fraga-Encinas

, P. Friberg

C. M. Fromm

, J. L. Gómez

, P. Galison

7 ,44 ,45

, C. F. Gammie

46 ,47

, R. García

, O. Gentaz

, B. Georgiev

26 ,27

, C. Goddi

1 ,9

, R. Gold

98 ,4 ,25

M. Gu (

顾

敏

峰

)

37 ,48

, M. Gurwell

, K. Hada

40 ,41

, M. H. Hecht

, R. Hesper

, L. C. Ho (

何

子

山

)

50 ,51

, P. Ho

, M. Honma

40 ,41

C. L. Huang

, L. Huang (

黄

磊

)

37 ,48

, D. H. Hughes

, S. Ikeda

14 ,53 ,54 ,55

, M. Inoue

, S. Issaoun

, D. J. James

7 ,6

, B. T. Jannuzi

, B. Jeter

26 ,27

W. Jiang (

江

悟

)

, M. D. Johnson

7 ,6

, S. Jorstad

56 ,57

, T. Jung

28 ,29

, M. Karami

25 ,26

, R. Karuppusamy

, T. Kawashima

, G. K. Keating

M. Kettenis

, J. Kim

, M. Kino

14 ,58

, 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 (

李

志

远

)

62 ,63

, M. Lindqvist

, R. Lico

, K. Liu

, E. Liuzzo

W. Lo

17 ,65

, A. P. Lobanov

, L. Loinard

66 ,67

, C. Lonsdale

, R. Lu (

路

如

森

)

37 ,38 ,16

, N. R. MacDonald

, J. Mao (

毛

基

荣

)

68 ,69 ,70

S. Marko

5 ,71

, D. P. Marrone

, A. P. Marscher

, I. Martí-Vidal

, S. Matsushita

, L. D. Matthews

, L. Medeiros

99 ,20

,???

, K. M. Menten

Y. Mizuno

, I. Mizuno

, J. M. Moran

7 ,6

, K. Moriyama

13 ,40

, M. Moscibrodzka

, C. Müller

16 ,1

, H. Nagai

14 ,41

, N. M. Nagar

, M. Nakamura

R. Narayan

7 ,6

, G. Narayanan

, R. Neri

, C. Ni

26 ,27

, A. Noutsos

, H. Okino

40 ,74

, H. Olivares

, G. N. Ortiz-León

, T. Oyama

, F. Özel

D. C. M. Palumbo

7 ,6

, N. Patel

, U. Pen

25 ,75 ,76 ,77

, D. W. Pesce

7 ,6

, 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

7 ,6

, L. Rezzolla

B. Ripperda

79 ,80

, A. Rogers

, E. Ros

, M. Rose

, A. Roshanineshat

, H. Rottmann

, A. L. Roy

, C. Ruszczyk

, B. R. Ryan

81 ,82

K. L. J. Rygl

, S. Sánchez

, D. Sánchez-Arguelles

53 ,84

, M. Sasada

40 ,85

, T. Savolainen

16 ,86 ,87

, F. P. Schloerb

, K. Schuster

, L. Shao

16 ,51

Z. Shen (

沈

志强

)

37 ,38

, D. Small

, B. Won Sohn

28 ,29 ,88

, J. SooHoo

, F. Tazaki

, P. Tiede

26 ,27

, M. Titus

, K. Toma

89 ,90

, P. Torne

16 ,83

E. Traianou

, T. Trent

, S. Trippe

, S. Tsuda

, H. J. van Langevelde

11 ,92

, D. R. van Rossum

, J. Wagner

, J. Wardle

, J. Weintroub

7 ,6

N. Wex

, R. Wharton

, M. Wielgus

7 ,6

, G. N. Wong

46 ,81

, Q. Wu (

吴

庆

文

)

, A. Young

, K. Young

, Z. Younsi

95 ,4

, F. Yuan (

袁

峰

)

37 ,48 ,96

Y. Yuan (

袁

业

飞

)

, J. A. Zensus

, G. Zhao

, S. Zhao

1 ,62

, and Z. Zhu

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 M 87, 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: interferometric

These authors contributed equally to this work.

NASA Hubble Fellowship Program, Einstein Fellow.

???

NSF Astronomy and Astrophysics Postdoctoral Fellow under award

no. AST-1903847.

Article published by EDP Sciences

A5, page 1 of 19

A&A 636, A5 (2019)

1. Introduction

The giant elliptical galaxy M 87 hosts an active galactic nucleus

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

2000). The radio core of M 87 shifts inwards with increasing fre-

quency as the jet becomes optically thin closer to the central black

hole, resulting in a flat radio spectrum as predicted by analytical

models (Blandford & Königl 1979; Falcke & Biermann 1995).

The radio core of M 87 coincides with the central engine at

43 GHz (Hada et al. 2011). At millimetre wavelengths, emis-

sion 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; Mo

scibrodzka et al. 2016). The shadow is a region exhibit-

ing 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 constant,

is the black hole mass, and

is the speed of light. The di

er-

ence 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). Esti-

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

tre of M 87 have historically ranged between (3

−

)

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

from stellar-dynamical measurements (Gebhardt

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

5) Mpc (Bird et al. 2010),

the mass measurements correspond to an apparent diameter of the

shadow between

∼

as and 42

as.

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

tion of

∼

as, which is certainly su

ffi

cient to resolve the

black hole shadow of M 87 for the high-mass estimate. M 87

is therefore one of the prime targets of the Event Horizon Tele-

scope (EHT), the Earth-sized mm-Very Long Baseline Interfero-

metry (VLBI) array aiming to image a black hole shadow (Event

Horizon Telescope Collaboration 2019a). The other prime candi-

date 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 scattering e

ects and vari-

ability 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; Johnson 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 M 87, 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 2019b). 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 2019c).

In 2017, the EHT consisted of the IRAM 30 m telescope

on Pico Veleta in Spain, the Large Millimeter Telescope (LMT)

in Mexico, the Atacama Large Millemeter Array (ALMA), the

Atacama Pathfinder Experiment (APEX) telescope 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 M 87 black hole

shadow within a 42

as asymmetric emission ring (Event

Horizon Telescope Collaboration 2019d,e). 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 2019e). At the adopted dis-

tance 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 Astro-

nomical 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

several 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

M 87. Fish et al. (2014) demonstrate that for Sgr A*, the blur-

ring e

ect of interstellar scattering could be mitigated if the

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

showed that source variability could also be mitigated by observ-

ing the source for multiple epochs and applying visibility aver-

aging, normalization, and smoothing to reconstruct an image of

the average 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

M 87, 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 scattered using a variable refractive scattering

screen, and the scattering can be mitigated by solving for the

scattering screen and image simultaneously (Johnson 2016).

However, scattering e

ects are only relevant for observations of

Sgr A*.

eht-imaging

can also simulate observations following

a real observing schedule, and copy the

-coverage and thermal

noise directly from existing data sets. It also includes polarimet-

ric leakage corruptions (Event Horizon Telescope Collaboration

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

http://www.aips.nrao.edu

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

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

A5, page 2 of 19

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

delays, phase decoherence due to atmospheric turbulence, and

signal attenuation caused by the atmospheric opacity are per-

fectly calibrated. In

eht-imaging

, atmospheric turbulence can

be included by fully randomizing the phases (with the option

of fixing them within a scan). In real mm-VLBI data, atmo-

spheric turbulence 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 aver-

age complex visibilities coherently on time scales set by the

atmospheric coherence time.

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., in prep.), which includes

a tropospheric module and physically motivated antenna point-

ing o

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

new

CASA

(McMullin et al. 2007) VLBI data calibration pipeline

rPICARD

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

amplitude calibration (Sect. 3). The overall computing workflow

SYMBA

is outlined in Sect. 4. We describe our simulated obser-

vational setup (antenna and weather parameters and observing

schedule) in Sect. 5 and our input source models for the syn-

thetic data generation in Sect. 6. In Sect. 7, we demonstrate the

ects of simulated data corruptions and subsequent calibration.

We illustrate the capabilities of

SYMBA

in Sect. 8 based on three

scientific case studies. In these studies we show (1) how well we

can distinguish between two example general relativistic magne-

tohydrodynamics (GRMHD) models with di

erent descriptions

for the electron temperatures with the current and future EHT

array, (2) how the EHT would perform under di

erent weather

conditions, and (3) how pre-2017 models of M 87 compare to the

observed image of the black hole shadow. In Sect. 9, we summa-

rize our conclusions and discuss future work.

2. Synthetic data generation with

MeqSilhouette

MeqSilhouette

(Blecher et al. 2017; Natarajan et al., in prep.)

is a synthetic data generator designed to simulate high fre-

quency VLBI observations. While visibilities of real radio inter-

ferometric observations are produced by correlating recorded

voltage streams from pairs of telescopes,

MeqSilhouette

pre-

dicts visibilities directly from the Fourier Transform of an

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

of Gaussian components, each with an independent spectral

index),

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 formalism

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

ects on

the visibilities into a chain of complex 2

2 Jones matrices

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

MeqSilhouette

generates

frequency-resolved visibilities, with a bandwidth and number

of channels set by the user. Frequency-resolved visibilities are

required for the calibration of signal path variations intro-

duced 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

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 the introduction of frequency dependent leakage of polar-

ized signals at telescopes’ receivers, the inclusion of wavelength

dependent Faraday rotation and spectral indices in source mod-

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

significant improvements to the

-coverage

. It is also possible

to generate corrupted data sets from time-dependent polarized

emission models in full Stokes and to follow an observed sched-

ule from a VEX file

. A key design driver of

MeqSilhouette

is 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 visibil-

ities can also be exported to UVFITS

. We briefly describe the

added tropospheric and instrumental corruptions below, refer-

ring to Blecher et al. (2017) and Natarajan et al. (in prep.) 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

amplitudes due to absorption of radiation in molecular transi-

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

tion 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

(Carilli & 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 coef-

ficients, respectively. In thermodynamic equilibrium, the latter

are related through Kirchho

’s law,

(

)

(2)

where

(

) is the Planck spectrum at temperature

. In order to

integrate Eq. (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 densities

of the lower and upper states, and the line shape including pres-

sure and Doppler broadening.

is related to the refractive index

of the medium via the Kramers-Kronig relations. The introduced

time delay is then calculated from the refractive index.

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

separated by 18 GHz, which

MeqSilhouette

can replicate.

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.

A5, page 3 of 19

A&A 636, A5 (2019)

As is evident from Kirchho

’s law (Eq. (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.

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

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

assuming a planar rather than a spherical atmosphere, which breaks

down at elevations below

∼

◦

(Paine 2019). For the synthetic obser-

vations in this work, we set the elevation limit to 10

◦

as is typically

done for real VLBI observations. Hence, the csc (elevation) approxima-

tion has a negligible e

ect on our results.

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)

where

int

is the data integration time and

is taken as the

lowest frequency in the data.

MeqSilhouette

uses Eq. (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 Eq. (4). The specified

coherence time

(PWV

) should decrease with increasing

precipitable water vapour content in the atmosphere, although

other factors such as wind speed also a

ect

. No sudden 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 integra-

tion 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 decoherence

ects within individual frequency channels are also not simu-

lated. It is assumed that frequency resolution is su

ffi

ciently high

to make this e

ect negligible, as it is done in modern correlators.

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

Jansky (Jy) follows as

DPFU

(5)

Here, the DPFU is the telescope’s “degree per flux unit” gain,

defined as DPFU

dish

(

)

, with

the aperture e

ffi

ciency (taken to be constant during observations),

dish

the geo-

metric 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 inte-

gration 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

data is a factor

√

smaller.

A5, page 4 of 19

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

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

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

FW H M

, 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







FW H M













(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

observed

corrupted (obs) visibilities from a baseline of stations

and

-terms cause artificial instrumental polarization as a

rotation 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 i

√

−

1. The star denotes complex con-

jugation. More complex and realistic polarimetric e

ects are

available in the forthcoming release of

MeqSilhouette v2

(Natarajan et al., in prep.).

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. 2019a). 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. 2019a).

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

5 on all base-

lines for which the source can be detected, is chosen (Janssen

et al. 2019b). 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

delays and absence of instrumental delays in the synthetic data

warrants a combined fringe fit solution over the whole frequency

band for a maximum S

N. Usually,

rPICARD

would smooth

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

already takes out the bulk of the delay o

sets. For the syn-

thetic 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. 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 s inter-

vals and a “network calibration” (Fish et al. 2011; Johnson &

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

http://astrogeo.org/psolve/

A5, page 5 of 19

A&A 636, A5 (2019)

Integration time [seconds]

S/N

sqrt S/N increase

= 15

= 2

= 1

Fig. 1.

N estimates for

rPICARD

fringe solutions. The plotted

points indicate the estimated average FFT S

N values by the

CASA

fringefit

code for di

erent integration times (solution intervals), seg-

menting a 15 min long scan of a

MeqSilhouette

observation of a 4 Jy

point source on the ALMA-APEX baseline. Di

erent symbols corre-

spond to di

erent coherence times (Eq. (4)) used for the simulation of

atmospheric turbulence. The dashed line shows the expected increase in

N for an infinite coherence time without added noise corruptions.

Collaboration 2019f) is performed with the

eht-imaging

soft-

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

a very short baseline to another nearby station can be calibrated

if the model of the observed source is known at large scales. For

the 2017 EHT observations, ALMA was able to provide accu-

rate large-scale source models, allowing for a network calibra-

tion 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

parameters are also set in ASCII files. The input source model

can be provided as FITS or ASCII file, as a single model or

multiple frames from a time-variable source, and contain only

Stokes

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

out-

puts a FITS file of the final reconstructed source model, the cali-

brated and self-calibrated visibilities in UVFITS and ASCII for-

mat, and diagnostic plots of the calibration process. The pipeline

is fully dockerized

. An overview of the workflow is shown in

Fig. 3.

5. Simulated observation setup

SYMBA

is able to create synthetic observations for any VLBI

array. Here, we outline the antenna and weather parameters and

https://www.docker.com/

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.

observing 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 M 87 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. (2019a). Full width at half maximum 230 GHz

beam sizes (

FW H M

) and dish diameters (

) were taken from the

websites and documentation for each individual site. Pointing

rms o

sets (

rms

) have been based on a priori station informa-

tion and typical inter- and intra-scan amplitude variations seen

in EHT data. All o

sets lie within o

ffi

cial telescope specifi-

cations. Aperture e

ffi

ciencies (

) were estimated with

∼

10%

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

Horizon Telescope Collaboration 2019f). In our synthetic obser-

vations, we have added gain errors (

err

) listed in Table 1 in

accordance with these uncertainties. Additionally, a polarization

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 2019f).

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 EHT

2017 campaign (5

−

11 April) at the individual primary sites, logged

by the VLBI monitor (Event Horizon Telescope Collaboration

2019a). No weather information was available from the VLBI mon-

itor for ALMA. We adopted the values measured at the nearby sta-

tion 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;

A5, page 6 of 19

F. Roelofs 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

(m)

(Jy)

err

D P

rms

(

′′

)

FW H M

(

′′

)

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

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

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 mixing

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, southern

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

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

observing season). We then interpolate these to the pressure level

of each EHT site and calculate the PWV from the measured

using Eq. (11).

Atmospheric coherence times

were estimated based on

the characteristics of the 2017 EHT measurements for the pri-

mary array. Precise station-based coherence times are di

ffi

cult

to obtain 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

As available on

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

cookbook/unix/zonal/

A5, page 7 of 19

A&A 636, A5 (2019)

Table 2.

Weather parameters adopted in our synthetic observations.

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

Table 2. A larger parameter space will be studied in future work

to characterize 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

located 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 h, and is distributed on a grid having 0.625

◦

longitude

by 0.5

◦

latitude with 42 vertical pressure levels between 0.1 and

1000 mbar. From this dataset, we took the 25th percentile (repre-

senting good weather) of the air temperature and specific humid-

ity measured on 11 April in the last two decades (1999

−

2018)

At each pressure level, these quantities were then linearly inter-

polated between the four grid points nearest to the observatory

site. We then performed an integration of the humidity over the

pressure levels using the

am

atmospheric model software (Paine

2019) to obtain the total PWV above the site. The starting point

for the integration over the pressure levels was determined 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 result-

ing 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 integration

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

It should be noted that the current EHT observing strategy is to trig-

ger a few observing days in a March

April observing window, based on

optimal weather conditions across all sites (Event Horizon Telescope

Collaboration 2019a).

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 M 87. Di

erent colors show baselines

within the EHT2017 array, baselines between the EHT2017 array and

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

(labelled as “Intra-new”).

ten 15 m dishes, so the sensitivity was scaled accordingly from

the JCMT, including a phasing e

ffi

ciency of 87%. The currently

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

Submillimetre Telescope (SEST, Booth et al. 1989) telescope in

Chile. With a sideband separating receiver, the current estimate

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

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 M 87 for the EHT2017

array and expansions with future stations. The EHT2017 sched-

ule for 11 April was adopted. To accommodate the eastward

expansion of the array with the AMT and PDB, ten-minute

scans were prepended to the schedule at 30 min intervals start-

ing 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

baseline. These gaps on short

-spacings pose challenges for

image reconstruction with the EHT2017 array (Event Horizon

Telescope Collaboration 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 M 87 for the EHT2017 and

EHT2020 arrays, and the extended schedule described above is

used for EHT2020

AMT array.

6. Source models

This section describes the set of input source models we use to

exercise the various aspects of the pipeline and perform scientific

case studies.

A5, page 8 of 19