MeqSilhouette : A mmVLBI observation and signal corruption simulator
Abstract
The Event Horizon Telescope (EHT) aims to spatially resolve the silhouette (or shadow) of the supermassive black holes in the Galactic Centre (Sgr A) and M87. The primary scientific objectives are to test general relativity in the strongfield regime and to probe accretion and jetlaunch physics at eventhorizon scales. This is made possible by the technique of Very Long Baseline Interferometry (VLBI) at (sub)millimetre wavelengths, which can achieve angular resolutions of order arcsec. However, this approach suffers from unique observational challenges, including scattering in the troposphere and interstellar medium; rapidly timevariable source structure in both polarized and total intensity; as well as nonnegligible antenna pointing errors. In this, the first paper in a series, we present the MeqSilhouette software package which is specifically designed to accurately simulate EHT observations. It includes realistic descriptions of a number of signal corruptions that can limit the ability to measure the key science parameters. This can be used to quantify calibration requirements, test parameter estimation and imaging strategies, and investigate systematic uncertainties that may be present. In doing so, a stronger link can be made between observational capabilities and theoretical predictions, with the goal of maximising scientific output from the upcoming order of magnitude increase in EHT sensitivity.
keywords:
instrumentation: interferometers, submillimetre: general, Galaxy: centre, atmospheric effects, techniques: high angular resolution1 Introduction
The principal goal of the Event Horizon Telescope (EHT) is to spatially resolve the gravitationallylensed photon ring (or ‘silhouette’) of a supermassive black hole (Doeleman et al., 2010). The two primary targets are Sgr A and M87, which are expected to have gravitationallylensed photon rings with apparent angular diameters of and arcsec (Falcke & Markoff, 2013; Broderick & Loeb, 2009) respectively. The extreme angular resolution required, blurring effects due to scattering by the interstellar medium (ISM; e.g. Fish et al., 2014), and the transition to an optically thin inner accretion flow at (sub)mmwavelengths (Serabyn et al., 1997; Falcke et al., 1998), necessitates that the EHT be comprised of antennas separated by thousands of kilometres that operate at high radio frequency ( GHz).
Performing what is known as Very Long Baseline Interferometry (VLBI) at mmwavelengths presents unique calibration challenges, including short atmospheric coherence times that are typically 10 s (Doeleman et al., 2009), low calibrator source sky density, complex and variable calibrator source structure, and antenna pointing accuracies that are a nonnegligible fraction of the antenna primary beam. These effects may place significant limitations on the sensitivity, image fidelity, and dynamic range that can be achieved with mmVLBI. Furthermore, unaccounted for systematic and/or nonGaussian uncertainties could preclude robust, accurate Bayesian parameter estimation and model selection analyses of accretion flow (e.g. Broderick et al., 2016) and gravitational physics (e.g. Broderick et al., 2014; Psaltis et al., 2016), two of the EHT’s many objectives.
Over the past decade, significant effort has been placed on advanced radio interferometric calibration and imaging algorithms for centimetre and metrewave facilities in response to the large number of new arrays in construction or design phase, including MeerKAT, Australian Square Kilometre Array Pathfinder (ASKAP), Square Kilometre Array (SKA), LowFrequency Array (LOFAR), and the Hydrogen Epoch of Reionization Array (HERA). A leading software package in this pursuit is MeqTrees^{1}^{1}1https://skasa.github.io/meqtrees/ (Noordam & Smirnov, 2010), which was developed to simulate, understand and address the calibration issues to be faced with the greatly enhanced sensitivity, instantaneous bandwidth, and fieldofview of such facilities. MeqTrees is rooted in the Measurement Equation mathematical formalism (Hamaker et al., 1996), which parametrises the signal path into distinct complex matrices called Jones matrices. This formalism and applications thereof are laid out in (Smirnov, 2011b, c, d) and are arbitrarily generalized to model any (linear) effect, including undesired signal corruptions that often may have subtle, yet systematic effects. MeqTrees has been applied to correct for direction dependent calibration errors to Karl. G. Jansky Very Large Array (VLA) and Westerbork Synthesis Radio Telescope (WSRT) observations, achieving recordbreaking high dynamic range images (Smirnov, 2011d). The effectiveness provided by the Measurement Equation formalism in radio interferometric calibration provides a strong motivation to explore its application to the challenging goal of imaging a supermassive black hole silhouette with mmVLBI.
Recently, there has been an increase in the attention given to simulating EHT observations of Sgr A and M87 (Fish et al., 2014; Lu et al., 2014; Bouman et al., 2015; Lu et al., 2016; Chael et al., 2016). However, these are primarily focused on image reconstruction and assume either negligible or Gaussian distributed gain errors; perfect antenna pointing accuracy; and in most cases only Gaussian convolution to simulate ISM scattering. Clearly, as the EHT array is enhanced (and possibly expanded), so too must the interferometric simulations evolve to provide evermore physical predictions on the confidence levels with which parameters can be extracted and hence exclude theoretical models of gravity and/or accretion flows.
Given the significant, yet surmountable, observational challenges that the EHT faces, we have undertaken a project to leverage metre and cmwavelength simulation and calibration successes and build a MeqTreesbased mmVLBIspecific software package called MeqSilhouette. While MeqTrees has not yet been used in the context of mmwavelength observations, the framework is agnostic to higher frequency implementation as long as the Measurement Equation is appropriately constructed. MeqSilhouette enables realistic interferometric simulations of mmVLBI observations in order to gain deeper understanding of a wide range of signal propagation and calibration effects. In this paper we describe the simulation framework and illustrate some of its key capabilities. These include the ability to simulate tropospheric, ISM scattering, and timevariable antenna pointing error effects. As will be demonstrated in a forthcoming series of papers, this technology will enable deeper understanding of a wide range of mmVLBI signal propagation and calibration systematics, quantify their effect on accretion flow and gravitational theoretical model selection, and hence maximise the scientific utility from EHT observations.
The paper is organized as follows: in section 2 we provide an introductory discussion on scattering theory; in section 3 we describe the implementation of the simulator and provide demonstrations of the most important modules; section 4 summarises our results and outlines our current plan for investigations with and future implementations of MeqSilhouette; and finally we conclude in section 5.
2 Theoretical background
Millimetre wavelength radiation originating at the Galactic Centre is repeatedly scattered along the signal path to the Earthbased observer. The first occurrence is due to electron plasma in the ISM (e.g. Bower et al., 2006; Gwinn et al., 2014), while the second is due to poorlymixed water vapour in the Earth’s troposphere (e.g. Carilli & Holdaway, 1999; Lay, 1997). It is essential that the effects of the scattering phenomena are understood for accurate calibration and robust inference of the intrinsic source properties. To this end, simulation modules approximating scattering in both media are implemented in MeqSilhouette. As an introduction to the separate descriptions of each, we review a simple scattering model.
An electromagnetic wave is scattered when it passes through a medium with refractive index inhomogeneities. Following Narayan (1992), this effect can be modeled as a thin screen, located between source and observer planes and orientated perpendicular to the lineofsight. The screen, indexed by coordinate vector , adds a stochastic phase to the incoming wave at each point on the screen, yielding a corrugated, outgoing wavefront. We define the Fresnel scale as , where is the observerscatterer distance, or the distance where the geometrical path difference rad.
To determine the resultant electric field at a point in the plane of the observer, indexed by coordinate vector , one has to take into account all possible ray paths from the screen to . To illustrate the model, a calculation of the scalar electric field generated by a point source, yields the FresnelKirchoff integral (Born & Wolf, 1980)
(1) 
where C is a numerical constant.
The statistical properties of can be described by a power spectrum or equivalently the phase structure function,
(2) 
where and represent two points on the screen and denotes the ensemble average.
can be reasonably approximated by a power law dependence on the absolute distance between points on the screen
(3) 
where is the phase coherence length scale defined such that rad.
Kolmogorov turbulence, which describes how kinetic energy injected at an outer length scale cascades to increasingly smaller scales until finally dissipated at an inner length scale , predicts in the domain . This scaling has been demonstrated to be a reasonable approximation for the ISM over scales km to AU (Armstrong, Rickett & Spangler, 1995), and also for the troposphere with , where is the thickness of the turbulent layer (Coulman, 1985; Carilli & Holdaway, 1997).
The two length scales, and , define the nature of the scattering which is split into the strong and weak regimes. In weak scattering, and hence by equation (3), . This implies that most of the radiative power measured on a point will originate from a screen area . Whereas in the regime of strong scattering, yielding . This results in coherent signal propagation onto the point from multiple disconnected zones each of area (Narayan, 1992). Scattering at millimetre wavelengths in the troposphere and the ISM in the direction of the Galactic Centre fall into the regimes of weak and strong scattering respectively.
To evolve the screen in time, one typically assumes a frozen screen i.e. that the velocity of the individual turbulent eddies is dominated by the bulk motion of scattering medium (e.g. Lay, 1997). This allows us to treat the screen as frozen but advected over the observer by a constant motion. Hence, time variability can be easily incorporated by the relative motion between source, scattering screen and observer.
3 Central components and layout of the simulator
MeqSilhouette is an observation and signal corruption simulator written in Python and MeqTrees using the measurement set^{2}^{2}2https://casa.nrao.edu/Memos/229.html data format. A flow diagram of the simulator is shown in Fig. 1. Input to the simulator is a sky model and configuration file. The former is typically a timeordered list of fits images, where each image represents the source total intensity^{3}^{3}3Later versions of MeqSilhouette will enable the full Stokes cubes as input. over a time interval , where is the observation length and is the number of source images. The configuration file specifies all parameters needed by the pipeline to determine the particular observation configuration (array, frequency, bandwidth, start time, etc.) and which signal corruption implementation should be employed. The visibilities are calculated through evaluation of the Fourier Transform at each UVW coordinate in the dataset, the time and frequency resolution of which is specified by the user. The primary outputs of the pipeline are an interferometric dataset in measurement set format along with the closure phases and uncertainties and a dirty and/or deconvolved image (or spectral cube if desired). The modular structure of the pipeline allows for multiple imaging and deconvolution algorithms to be employed. The rest of this section is devoted to describing the implementation of each signal corruption module.
3.1 Interstellar medium
Scattering in the ISM at millimetre wavelengths towards the Galactic Centre falls into the strong scattering regime, which can be further subdivided into snapshot, average and ensembleaverage regimes (Narayan & Goodman, 1989; Goodman & Narayan, 1989). Following from section 2, patches on the scattering screen with linear size will emit electromagnetic waves into singleslit diffraction cones of angular size . For a point source, an observer will be illuminated by many patches spanning with projected size on the screen equal to the refractive scale,
(4) 
The diffraction cones from each small region will interfere, resulting in a multislit diffractive scintillation pattern. A single realisation of this pattern falls in the snapshot regime. An extended source will average over many realisations and quench the diffractive scintillation. In the average regime, although diffractive scintillation has been averaged over, there still exists scintillation over scales comparable to the size of the scattered image of a point source , termed refractive scintillation. This scintillation acts to focus/defocus the ensemble of coherent patches of linear size . This weak, largescale scintillation is more difficult to average over, requiring multiepoch observations over weeks to months in order to allow the scattering material to move across the source (assuming transverse ISM velocities of a few 10s of km s). An extended source size will quench refractive fluctuations but only when . In the ensembleaverage regime, all scintillation has been averaged and the scattering is equivalent to Gaussian convolution.
An algorithm which approximates scattering in the average regime, which is relevant to VLBI observations of Sgr A, has been implemented in the Pythonbased Scatterbrane^{4}^{4}4http://krosenfeld.github.io/scatterbrane package, based on Johnson & Gwinn (2015). This approach extends the structure function shown in equation (3) to regimes where the inner and outer turbulent scales as well as the anisotropy of scattering kernel are considered. In this framework the scattered image is approximated by ‘reshuffling’ of the source image through
(5) 
where is the directional derivative. Even though is only coherent to , remains spatially coherent over much larger scales, leading to the presence of refractive substructure (Johnson & Gwinn, 2015).
We include the ScatterBrane software, which has already yielded important context for mmVLBI observations towards Sgr A (e.g. OrtizLeón et al., 2016), within the MeqSilhouette simulation framework.Our ISM module interfaces the Scatterbrane code within an interferometric simulation pipeline. This module enables simultaneous use of timevariable ISM scattering and timevariable intrinsic source structure within a single framework. The user is able to select a range of options relating to the timeresolution and epoch interpolation/averaging of both. By default, if the time resolution chosen to sample the source variability and screen variability are unequal, we set

if

if ,
where rounds the fraction to the nearest integer. This modification to the ISM sampling resolution avoids interpolation between different snapshots of the intrinsic source structure. Note that even though the ISMscattering corruption is applied in the correct causal position in the signal propagation chain, equation (5) is nonlinear and hence can not be written in the Measurement Equation formalism.
To demonstrate the implementation and provide an example of intraday ISM variability, we present the results of a simulated observation of 10 minutes duration at 14:00 UTC on four consecutive days in Fig. 2. To compare to published observations, we use the threestation EHT array consisting of the Submillimeter Telescope (SMT) in Arizona, the Combined Array for Research in Millimeterwave Astronomy (CARMA) in California and the James Clerk Maxwell Telescope (JCMT) on Mauna Kea, Hawaii. The distance to the screen is taken as kpc (Bower et al., 2013). The relative transverse velocity between the observer and scattering screen is set to to be consistent with OrtizLeón et al. (2016). The source is a circular Gaussian with a arcsec, approximately the angular distance that a scattering screen would travel over days. The source size has been chosen such that it is consistent with the latest estimate of the size of Sgr A at GHz (Fish et al., 2011). Closure quantities are model dependent and calculated as specified in Rogers et al. (1995), where the thermal noise was added based on the system equivalent flux density (SEFD) table in (Lu et al., 2014).
3.2 Pointing Errors
All antennas suffer pointing errors to some degree as a result of a variety of factors, including dish flexure due to gravity, wind and thermal loading, as well as imperfect drive mechanics. This corresponds to an offset primary beam, which should only translate to minor amplitude errors if the pointing error is significantly smaller than the primary beam (i.e. ). In the Measurement Equation formalism, this offset can be represented by a modified (shifted) primary beam pattern in the EJones term
(6) 
where correspond to the directional cosine offsets.
We investigate the effect of pointing errors on the 50 m (i.e. fully illuminated) Large Millimeter Array (LMT) dish configured in an eight station VLBI array. The LMT has been measured to have an absolute pointing accuracy of arcsec, where smaller offsets occur when observing sources closer to zenith, and a tracking pointing accuracy arcsec^{5}^{5}5http://www.lmtgtm.org/telescope/telescopedescription/. We investigate the observational effect of these errors through three different pointing error models which explore different instructive and plausible scenarios. The LMT has been singled out as this may well serve as a reference station for the EHT array given its sensitivity and central geographic location. The source used is a circular Gaussian of characteristic size arcsec, located at the phase centre. For this investigation, as long as , the exact structure of the source is unimportant. We approximate the LMT beam profile using an analytic WSRT beam model (Popping & Braun, 2008) with a factor of two increase in the beam factor to take into account the increased dish size
(7) 
where is a constant, with value GHz. Note that the power beam becomes , resulting in a arcsec at 230 GHz. We make use of the RMS fractional visibility amplitude error , where and are the visibility amplitudes with and without pointing errors respectively, and . In Fig. 3, is plotted against pointing error over the range arcsec.
In the first case we assume a constant pointing error. This simulation is meant to be instructive as to the typical amplitude error in the simplest possible scenario.
Also interesting to consider is a slower, continuous timevariable pointing error associated with the tracking error . Physically, this could be attributed to changes in wind, thermal and gravitational loading which all change with telescope pointing direction and over the course of a typical few hour observation. Using the MeqTrees software package, such behaviour has been demonstrated to occur with the WSRT (Smirnov, 2011a, d). This is modeled as sinusoidal variability with period sampled from a uniform distribution between 0.5 and 6 hours, and a peak amplitude , where the factor relates the peak amplitude to the RMS of a sinusoidal, zeromean waveform.
Whilst a stationary phase centre is tracked, the pointing error should evolve slowly and smoothly, however, in mmVLBI observations the phase centre is often shifted to another source/calibrator. This could cause the pointing error to change abruptly, with an absolute pointing error . Source/calibrator change is scheduled every 510 minutes in a typical millimetre observation. An important point is that even though the EHT will be able to determine the pointing offset when observing a calibrator with well known source structure, when the antennas slew back to a source (e.g. Sgr A) with less certain or variable source structure, the pointing error could change significantly. This is exacerbated by the scarcity of mmwavelength calibrators, which are often widely separated from the source. The antenna pointing error would induce scatter in the visibility amplitudes, which may be difficult to decouple from other effects e.g. intrinsic source variability and/or structure as well as time variable ISM scattering. We simulate this stochastic variability by resampling the pointing error every 10 minutes from a Gaussian of characteristic width equal to the quoted pointing error. We perform 50 realisations of the simulation for each pointing offset to generate reasonable uncertainties.
In this simulation, we only consider LMT pointing errors due to its narrow primary beam and potential to be used as a reference station. However, the capability to simulate independent pointing errors for each station is available. In the case of a phased array, a pointing error simulation could be used to investigate the contribution of the pointing error to a variable phasing efficiency, which can be reasonably approximated by a scalar Jones matrix.
3.3 Troposphere
The coherence and intensity of millimetre wavelength electromagnetic waves are most severely deteriorated in the lowest atmospheric layer, the troposphere, which extends up to an altitude of km above sea level and down to a temperature K (Thompson et al., 2001). The troposphere is composed of primary gases and , trace gases (e.g. water vapour and ), as well as particulates of water droplets and dust. The absorption spectrum in the GHz range (e.g. Pardo et al., 2001) is dominated by several transitions of and as well as a pseudocontinuum opacity which increases with frequency. The pseudocontinuum opacity is due to the far wings of a multitude of pressurebroadened water vapour lines above 1 THz (Carilli & Holdaway, 1999).
In contrast to the other atmospheric chemical components, water vapour mixes poorly and its timevariable spatial distribution induces rapid fluctuations in the measured visibility phase at short wavelengths. The water vapour column density is measured as the depth of the column when converted to the liquid phase and is referred to as the precipitable water vapour (PWV). The PWV is, via the real component of the refractive index, directly proportional to phase offset,
(8) 
where is the depth of the PWV column (Carilli & Holdaway, 1999) and an atmospheric temperature K has been assumed. This relationship between phase and water vapour content has been experimentally verified (Hogg et al., 1981). At 230 GHz, the change in PWV needed to offset the phase by 1 rad is mm. In the mmVLBI case, this sensitive dependence of phase coherence on atmospheric stability is aggravated by typically low antenna elevation angles, uncorrelated atmospheric variations between stations, and the sparsity of the array.
Our focus is to model three primary, interrelated observables which are the most relevant to mmVLBI: turbulencedriven fluctuations in the visibility phase ; signal attenuation due to the atmospheric opacity ; and the increase in system temperature due to atmospheric emission at a brightness temperature .
Our approach is to model these observables as being separable into mean and turbulent components which are simulated independently. The mean tropospheric simulation module performs radiative transfer with a detailed model of the electromagnetic spectrum of each atmospheric constituent. The turbulent simulation module takes a scattering approach to account for the decoherence that results from powerlaw turbulence.
3.3.1 Average Troposphere
The problem of radiative transfer through a static atmosphere is well described and implemented by the Atmospheric Transmission at Microwaves (atm) software (Pardo et al., 2001). atm has been incorporated into MeqSilhouette to provide a fast and sophisticated procedure to calculate average opacities, sky brightness temperatures and time delays. Here we provide a brief summary of the theory underpinning the package but refer the reader to Pardo et al. (2001) for more detail. atm is commonly used in the Atacama Large Millimeter Array (ALMA) community (Curtis et al., 2009; Nikolic et al., 2013) and has been tested with atmospheric transmission spectra taken on Mauna Kea (Serabyn et al., 1998).
We start from the unpolarised radiative transfer equation, which is unidirectional in the absence of scattering,
(9) 
where is the coordinate along the signal path through the atmosphere, is the specific intensity, is the macroscopic emission coefficient and is the macroscopic absorption coefficient.
The goal is to integrate this equation over the signal path which requires as a function of altitude and frequency. The integration naturally yields the mean opacity and sky brightness temperature. The mean time delay is calculated from using the KramersKronig relations. In practice, this involves a triple sum over altitude layer, chemical species and rotational energy transition. Atmospheric temperature and pressure profiles are calculated based on several station dependent inputs, namely, ground temperature and pressure and the precipitable water vapour column depth.
A general equation to determine the absorption coefficient for a transition between a lower and upper states is given in the original paper. Here we merely point out that it should be proportional to the energy of the photon, , the transition probability or Einstein coefficient, , the lineshape, and the number densities of electronic populations. Line profiles which describe pressure broadening (perturbations to the Hamiltonian due to the presence of nearby molecules) and Doppler broadening are used. The condition of detailed balance further requires that decays from the upper state are included yielding, , where is the degeneracy of the electronic state. Putting this together we find,
(10) 
where the Einstein coefficients are calculated from the inner product of the initial and final states with the dipole transition operator. The number densities of the two states, and in local thermodynamic equilibrium (LTE) are simply related to the local number density and temperature via Boltzmann statistics.
Typical opacities and sky brightness temperatures for ALMA, the Submillimeter Array (SMA) and the South Pole Telescope (SPT) are shown in Fig. 4. Note that both the opacity and brightness temperature are inversely proportional to the ground temperature and proportional to ground pressure.
3.3.2 Turbulent troposphere
Visibility phase instability due to tropospheric turbulence is a fundamental limitation to producing high fidelity, sciencequality maps with a mmVLBI array (Thompson et al., 2001). The coherence timescale is typically too rapid ( s) for fast switching calibration, so other calibration procedures (e.g. water vapour radiometry, paired antennas, and/or selfcalibration) must be performed. Selfcalibration is the most commonly used but is limited by the integration time needed to obtain adequate SNR to fringe fit. Phase decoherence often leads to the use of closure quantities to perform model fitting (e.g. Doeleman et al., 2001; Bower et al., 2004; Shen et al., 2005).
Following from section 2, we can model the statistics of with a thin, frozen, Kolomogorovturbulent phase screen moving at a bulk transverse velocity, . We set the height of the screen at the water vapour scale height of 2 km above ground. We will show later that the thickness of the atmospheric turbulent layer can be neglected in our implementation. At an observing wavelength of mm, the Fresnel scale is m and experiments show annual variations of m above Mauna Kea (Masson, 1994) and m above Chajnantor (Radford & Holdaway, 1998), where both sites are considered to have excellent atmospheric conditions for (sub)millimetre astronomy. As , this is an example of weak scattering.
The required fieldofview (FoV) of a global mmVLBI array is typically FoV mas or m at a height of 2 km, which is roughly 78 orders of magnitude smaller than the tropospheric coherence length. The tropospheric corruption can therefore be considered constant across the FoV and, from the perspective of the Measurement Equation, modeled as a diagonal Jones matrix per time and frequency interval. As VLBI baselines are much longer than the turbulent outer scale, km km, the phase variations are uncorrelated between sites and can be simulated independently. This assumption only holds for VLBI baselines and the framework needs to be extended to simulate the effects of turbulence on individual phased arrays stations (e.g. SMA) and short ( km) baselines (e.g. JCMT  SMA).
Our aim then is to produce a phase error time sequence for each station which is added to the visibility phase. We invoke the frozen screen assumption and write the structure function as a function of time, . The temporal structure function provides an efficient route to sample the variability of the troposphere at the typical integration time of the dataset, sec.
The temporal variance of the phase is a function of the temporal structure function, and accounting for time integration yields (see Treuhaft & Lanyi, 1987, B3)
(11) 
Assuming powerlaw turbulence and integrating yields,
(12) 
where is the coherence time when observing at zenith and is the approximate airmass which arises as . As , where is the thickness of the turbulent layer, an thin screen exponent of is justified (Treuhaft & Lanyi, 1987). The phase error timeseries takes the form of a Gaussian random walk per antenna. At mmwavelengths, the spectrum of water vapour is nondispersive up to a few percent (Curtis et al., 2009) and so we can assume a simple linear scaling across the bandwidth. Fig. 5 shows an example simulation of the turbulent and total delays at the SMA and ALMA sites.
Phase fluctuations can also be simulated by taking the inverse Fourier transform of the spatial phase power spectrum. However this approach is much more computationally expensive, e.g. for an observation length involving independent antennae with dish radii m, wind speed m s and pixel size equal to , the number of pixels . Additionally, due to fractal nature of ideal Kolmogorov turbulence, the power spectrum becomes unbounded as the wavenumber approaches zero which makes it difficult to determine the sampling interval of the spatial power spectrum (Lane et al., 1992).
3.3.3 Limitations to highfidelity image reconstruction
A primary objective of MeqSilhouette is to understand and constrain systematic errors in mmVLBI observations. In this section the tropospheric module is used to estimate the effect on image quality for various levels of calibration accuracy.
We simulate the simple scenario of a sky model that consists of a 2.4 Jy point source at the phase centre, which is the approximate EHTmeasured flux density of Sgr A at 230 GHz. We assume a zenith phase coherence time of s above each station (however, each stations PWV can be independently simulated). We approximate the effect of imperfect calibration by adding a small fraction of the turbulent phase noise. For this example, we do not include the mean delay component, assuming it to be perfectly corrected for during the calibration.
4 Discussion
In section 3 we have described the layout of MeqSilhouette synthetic data simulation framework. A wide range of signal propagation effects can be implemented using the Measurement Equation formalism, with tropospheric scattering and antenna pointing errors given as illustrative examples. The framework is sufficiently general and flexible so that time variability in all relevant domains (source, array, ISM, troposphere) can be incorporated. The run time for a typical simulation with a realistic instrumental setup is on the order of minutes. Implementation of polarisation effects is intended in the next version.
The ISM scattering software ScatterBrane, based on Johnson & Gwinn (2015), has been incorporated into the pipeline. Fig. 2 provides an example of closure phase and flux variability over a 4 day period using a static source. Accurate simulation of the ISMinduced closure phase variation is essential in order to make accurate inferences regarding asymmetric, eventhorizon scale structure from EHT observations (e.g. Fish et al., 2016; OrtizLeón et al., 2016). This will become even more important as the EHT sensitivity increases by an order of magnitude in the near future. Note that if the source position is time variable as in the case of a hotspot model (Doeleman et al., 2009), this will increase ISM variability as the relative motion between source, screen and observer is increased.
Visibility amplitude errors due to antenna pointing error has been investigated for the m LMT dish operating at GHz. In Fig. 3, we show that pointing errors associated with frequent phase centre switching (stochastic variability) could introduce a RMS fractional amplitude error for an absolute pointing accuracy arcsec. In contrast, tracking errors are less problematic with for a tracking accuracy arcsec. The case of a constant error pointing model is comparable to that of the ‘slow variability’ case. If the gain error is nonseparable from the calibration model used, it could be interpreted as intrinsic variability, substructure and/or increased noise. If unaccounted for, this effect has the potential to limit the dynamic range of mmVLBI images. Further tests to constrain the pointing uncertainties of EHT stations will lead to more accurate interferometric simulations and hence the overall impact on black hole shadow parameter estimation. Here we demonstrate the capability to incorporate a range of plausible pointing error effects into a full simulation pipeline. For future observations at 345 GHz, these effects will be even more pronounced, given the narrower primary beam.
In section 3.3.3 we explore the observational consequences of observing through a turbulent troposphere. In this simulation, we assume a simple point source model and apply increasing levels of turbulenceinduced phase fluctuations before imaging using regular sampling and a two dimensional inverse fast Fourier transform. The simulated residual calibration errors result in a significant attenuation in source flux; slight offsets in the source centroid (black crosshairs) and the presence of spurious imaging artefacts. In an upcoming paper, we perform a systematic exploration of the turbulent tropospheric effects on the accuracy of fringefitting algorithms and strategies, through use of an automated calibration procedure and including the added complexity of a timevariable source.
Significant progress has been made in the theoretical and numerical modeling of the inner accretion flow and jet launch regions near a supermassive black hole event horizon (e.g. Del Zanna et al., 2007; Etienne et al., 2010; Dexter & Fragile, 2013; Mościbrodzka et al., 2014; McKinney et al., 2014). As the sensitivity of the EHT stands to dramatically increase, these theoretical efforts must be complemented by advances in interferometric simulations. With MeqSilhouette, we now have the ability to couple these with sophisticated interferometric and signal propagation simulations. Moreover, detailed interferometric simulations will enable us to quantify systematic effects on the black hole and/or accretion flow parameter estimation.
5 Conclusion
In light of the science objectives of mmVLBI observations and software advances in the broader radio interferometry community, a mmVLBI data simulator has been developed. An important feature is that this simulation pipeline is performed using the Measurement Set format, in line with ALMA and future VLBI data formats. The focus has been placed on simulating realistic data given an arbitrary theoretical sky model. To this end, the simulator includes signal corruptions in the ISM, troposphere and instrumentation. Examples of typical corruptions have been demonstrated, which show that each corruption can significantly affect the inferred scientific parameters. Particular focus has been placed on EHT observations, however, the pipeline is completely general with respect to observation configuration and source structure. Time variability in all domains (source, array, ISM, troposphere) is implemented. Future versions of MeqSilhouette will include polarisation dependent corruptions. The creation of a close interface between sophisticated theoretical and interferometric mmVLBI simulations will enhance the scientific opportunities possible with the EHT.
Acknowledgements
We thank Michael Johnson and Katherine Rosenfeld for making the ScatterBrane code publicly available, and for helpful discussions. Similarly, we thank Bojan Nikolic for atm support. We thank the referee for helpful comments and questions which helped to refine this work. We are grateful to Monika Mościbrodwska for supplying theoretical simulations we used in testing and to Ilse van Bemmel and Lindy Blackburn for useful discussions. We thank Paul Galatis for assistance with the MeqSilhouette logo design. The financial assistance of the South African SKA Project (SKA SA) towards this research is hereby acknowledged. Opinions expressed and conclusions arrived at are those of the author and are not necessarily to be attributed to the SKA SA^{6}^{6}6www.ska.ac.za O. Smirnov’s research is supported by the South African Research Chairs Initiative of the Department of Science and Technology and National Research Foundation.
References
 Armstrong et al. (1995) Armstrong J. W., Rickett B. J., Spangler S. R., 1995, ApJ, 443, 209
 Born & Wolf (1980) Born M., Wolf E., 1980, Pergamon, New York
 Bouman et al. (2015) Bouman K. L., Johnson M. D., Zoran D., Fish V. L., Doeleman S. S., Freeman W. T., 2015, preprint, (arXiv:1512.01413)
 Bower et al. (2004) Bower G. C., Falcke H., Herrnstein R. M., Zhao J.H., Goss W. M., Backer D. C., 2004, Science, 304, 704
 Bower et al. (2006) Bower G. C., Goss W. M., Falcke H., Backer D. C., Lithwick Y., 2006, ApJ, 648, L127
 Bower et al. (2013) Bower G. C., et al., 2013, ApJ, 780, L2
 Broderick & Loeb (2009) Broderick A. E., Loeb A., 2009, ApJ, 697, 1164
 Broderick et al. (2014) Broderick A. E., Johannsen T., Loeb A., Psaltis D., 2014, ApJ, 784, 7
 Broderick et al. (2016) Broderick A. E., et al., 2016, ApJ, 820, 137
 Carilli & Holdaway (1997) Carilli C., Holdaway M., 1997, MMA/ALMA Memorandum Series No 173, NRAO
 Carilli & Holdaway (1999) Carilli C. L., Holdaway M. A., 1999, Radio Science, 34, 817
 Chael et al. (2016) Chael A. A., Johnson M. D., Narayan R., Doeleman S. S., Wardle J. F. C., Bouman K. L., 2016, preprint, (arXiv:1605.06156)
 Coulman (1985) Coulman C. E., 1985, Annual Review of Astronomy and Astrophysics, 23, 19
 Curtis et al. (2009) Curtis E. I., Nikolic B., Richer J. S., Pardo J. R., 2009, arXiv preprint arXiv:0912.2852
 Del Zanna et al. (2007) Del Zanna L., Zanotti O., Bucciantini N., Londrillo P., 2007, \aap, 473, 11
 Dexter & Fragile (2013) Dexter J., Fragile P. C., 2013, \mnras, 432, 2252
 Doeleman et al. (2001) Doeleman S. S., et al., 2001, \aj, 121, 2610
 Doeleman et al. (2009) Doeleman S. S., Fish V. L., Broderick A. E., Loeb A., Rogers A. E. E., 2009, ApJ, 695, 59
 Doeleman et al. (2010) Doeleman S., et al., 2010, in astro2010: The Astronomy and Astrophysics Decadal Survey. (arXiv:0906.3899)
 Etienne et al. (2010) Etienne Z. B., Liu Y. T., Shapiro S. L., 2010, \prd, 82, 084031
 Falcke & Markoff (2013) Falcke H., Markoff S. B., 2013, Class. Quantum Grav., 30, 244003
 Falcke et al. (1998) Falcke H., Goss W. M., Matsuo H., Teuben P., Zhao J.H., Zylka R., 1998, ApJ, 499, 731
 Fish et al. (2011) Fish V. L., et al., 2011, ApJ, 727, L36
 Fish et al. (2014) Fish V. L., et al., 2014, ApJ, 795, 134
 Fish et al. (2016) Fish V. L., et al., 2016, ApJ, 820, 90
 Goodman & Narayan (1989) Goodman J., Narayan R., 1989, \mnras, 238, 995
 Gwinn et al. (2014) Gwinn C. R., Kovalev Y. Y., Johnson M. D., Soglasnov V. A., 2014, ApJ, 794, L14
 Hamaker et al. (1996) Hamaker J. P., Bregman J. D., Sault R. J., 1996, \aaps, 117, 137
 Hogg et al. (1981) Hogg D., Guiraud F., Decker M., 1981, \aap, 95, 304
 Johnson & Gwinn (2015) Johnson M. D., Gwinn C. R., 2015, ApJ, 805, 180
 Lane (1998) Lane A. P., 1998, in Novak G., Landsberg R., eds, Astronomical Society of the Pacific Conference Series Vol. 141, Astrophysics From Antarctica. pp 289–295
 Lane et al. (1992) Lane R. G., Glindemann A., Dainty J. C., 1992, Waves in Random Media, 2, 209
 Lay (1997) Lay O. P., 1997, \aaps, 122, 535
 Lu et al. (2014) Lu R.S., Broderick A. E., Baron F., Monnier J. D., Fish V. L., Doeleman S. S., Pankratius V., 2014, ApJ, 788, 120
 Lu et al. (2016) Lu R.S., et al., 2016, \apj, 817, 173
 Masson (1994) Masson C. R., 1994, in IAU Colloq. 140: Astronomy with Millimeter and Submillimeter Wave Interferometry. pp 87–95
 McKinney et al. (2014) McKinney J. C., Tchekhovskoy A., Sadowski A., Narayan R., 2014, \mnras, 441, 3177
 Mościbrodzka et al. (2014) Mościbrodzka M., Falcke H., Shiokawa H., Gammie C. F., 2014, \aap, 570, A7
 Narayan (1992) Narayan R., 1992, Philosophical Transactions of the Royal Society A: Mathematical Physical and Engineering Sciences, 341, 151
 Narayan & Goodman (1989) Narayan R., Goodman J., 1989, \mnras, 238, 963
 Nikolic et al. (2013) Nikolic B., Bolton R. C., Graves S. F., Hills R. E., Richer J. S., 2013, \aap, 552, A104
 Noordam & Smirnov (2010) Noordam J. E., Smirnov O. M., 2010, \aap, 524, A61
 OrtizLeón et al. (2016) OrtizLeón G. N., et al., 2016, preprint, (arXiv:1601.06571)
 Pardo et al. (2001) Pardo J., Cernicharo J., Serabyn E., 2001, IEEE Trans. Antennas Propagat., 49, 1683
 Popping & Braun (2008) Popping A., Braun R., 2008, \aap, 479, 903
 Psaltis et al. (2016) Psaltis D., Wex N., Kramer M., 2016, ApJ, 818, 121
 Radford & Holdaway (1998) Radford S. J., Holdaway M. A., 1998, in Phillips T. G., ed., \procspieVol. 3357, Advanced Technology MMW, Radio, and Terahertz Telescopes. pp 486–494
 Rogers et al. (1995) Rogers A. E. E., Doeleman S. S., Moran J. M., 1995, \aj, 109, 1391
 Serabyn et al. (1997) Serabyn E., Carlstrom J., Lay O., Lis D. C., Hunter T. R., Lacy J. H., Hills R. E., 1997, ApJ, 490, L77
 Serabyn et al. (1998) Serabyn E., Weisstein E. W., Lis D. C., Pardo J. R., 1998, Appl. Opt., 37, 2185
 Shen et al. (2005) Shen Z.Q., Lo K. Y., Liang M.C., Ho P. T. P., Zhao J.H., 2005, \nat, 438, 62
 Smirnov (2011a) Smirnov O., 2011a, presentation at ‘SKA CALIM workshop : Solving for primary beams, pointing errors, and The Westerbork Wobble’, Manchester, 2529 July 2011., https://indico.skatelescope.org/event/171/session/9/contribution/20
 Smirnov (2011b) Smirnov O. M., 2011b, A&A, 527, A106
 Smirnov (2011c) Smirnov O. M., 2011c, A&A, 527, A107
 Smirnov (2011d) Smirnov O. M., 2011d, A&A, 527, A108
 Thompson et al. (2001) Thompson A. R., Moran J. M., Swenson Jr. G. W., 2001, Interferometry and Synthesis in Radio Astronomy, 2nd Edition
 Treuhaft & Lanyi (1987) Treuhaft R. N., Lanyi G. E., 1987, Radio Science, 22, 251