MeqSilhouette : A mm-VLBI simulator

MeqSilhouette : A mm-VLBI observation and signal corruption simulator

Tariq Blecher, Roger Deane, Gianni Bernardi, Oleg Smirnov
Centre for Radio Astronomy Techniques and Technologies, Department of Physics and Electronics, Rhodes University,-
Grahamstown 6140, South Africa
Square Kilometre Array South Africa, Pinelands 7405, Cape Town, South Africa

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 strong-field regime and to probe accretion and jet-launch physics at event-horizon 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 time-variable source structure in both polarized and total intensity; as well as non-negligible 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.

instrumentation: interferometers, submillimetre: general, Galaxy: centre, atmospheric effects, techniques: high angular resolution

1 Introduction

The principal goal of the Event Horizon Telescope (EHT) is to spatially resolve the gravitationally-lensed 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 gravitationally-lensed 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)mm-wavelengths (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 mm-wavelengths 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 non-negligible 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 mm-VLBI. Furthermore, unaccounted for systematic and/or non-Gaussian 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 metre-wave 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), Low-Frequency Array (LOFAR), and the Hydrogen Epoch of Reionization Array (HERA). A leading software package in this pursuit is MeqTrees111 (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 field-of-view 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 record-breaking 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 mm-VLBI.

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 ever-more 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 cm-wavelength simulation and calibration successes and build a MeqTrees-based mm-VLBI-specific software package called MeqSilhouette. While MeqTrees has not yet been used in the context of mm-wavelength observations, the framework is agnostic to higher frequency implementation as long as the Measurement Equation is appropriately constructed. MeqSilhouette enables realistic interferometric simulations of mm-VLBI 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 time-variable 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 mm-VLBI 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.

Figure 1: Flow diagram showing basic sequence of the MeqSilhouette simulation pipeline. The sky model could include (a) a time-ordered list of fits images or (b) parametric source model consisting of Gaussians or point sources. The details of the station information, observation strategy, tropospheric and ISM conditions are specified in a user-defined input configuration file. The pipeline is flexible, allowing any additional, arbitrary Jones matrices to be incorporated. Note that the current ISM-scattering implementation is non-linear and hence can not be incorporated into the Measurement Equation formalism. Further details in text.

2 Theoretical background

Millimetre wavelength radiation originating at the Galactic Centre is repeatedly scattered along the signal path to the Earth-based 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 poorly-mixed 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 electro-magnetic 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 line-of-sight. 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 observer-scatterer 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 Fresnel-Kirchoff integral (Born & Wolf, 1980)


where C is a numerical constant.

The statistical properties of can be described by a power spectrum or equivalently the phase structure function,


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


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 set222 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 time-ordered list of fits images, where each image represents the source total intensity333Later 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 ensemble-average regimes (Narayan & Goodman, 1989; Goodman & Narayan, 1989). Following from section 2, patches on the scattering screen with linear size will emit electromagnetic waves into single-slit 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,


The diffraction cones from each small region will interfere, resulting in a multi-slit 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, large-scale scintillation is more difficult to average over, requiring multi-epoch 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 ensemble-average 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 Python-based Scatterbrane444 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


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 mm-VLBI observations towards Sgr A (e.g. Ortiz-Leó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 time-variable ISM scattering and time-variable intrinsic source structure within a single framework. The user is able to select a range of options relating to the time-resolution 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 ISM-scattering corruption is applied in the correct causal position in the signal propagation chain, equation (5) is non-linear 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 three-station EHT array consisting of the Submillimeter Telescope (SMT) in Arizona, the Combined Array for Research in Millimeter-wave 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 Ortiz-Leó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).

Figure 2: An example simulation of ISM scattering towards Sgr A, observed with SMT-JCMT-CARMA. The top panel, left to right, shows the original  -arcsec Gaussian (top left), the simulated ISM scattered image on the first night (top middle) and last night (top right) of the observation, respectively. The bottom panel, left to right, shows the evolution of the 10 minute-averaged closure phase with epoch (bottom left), uv-tracks for each night (bottom middle) and the RMS fractional visibility amplitude differences as a function of uv-distance (bottom right). , where and — are the simulated average and ensemble average visibility amplitudes respectively. Variations from the ensemble-average flux on the shortest baselines reveal total flux modulation while flux variations on longer baselines and non-zero closure phases track the fluctuations in substructure. Furthermore, ISM scattering simulations can constrain the variability fraction associated with the screen, enabling a more robust estimation of source variability, as demonstrated in Ortiz-León et al. (2016). The time-variability of the ISM is built into the MeqSilhouette framework.

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 E-Jones term


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  arcsec555 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


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 time-variable 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, zero-mean waveform.

Whilst a stationary phase centre is tracked, the pointing error should evolve slowly and smoothly, however, in mm-VLBI 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 5-10 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 mm-wavelength 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 re-sampling 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.

Figure 3: RMS fractional amplitude error induced by pointing error with the 50 m (i.e. fully illuminated) LMT antenna as a function of pointing error offset at  GHz. We assume that these errors are degenerate or non-separable from the self-calibration/fringe-fitting model used. This simulation capability enables constraints on the magnitude of pointing-induced errors given a particular pointing calibration strategy. See text for more details.

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 pseudo-continuum opacity which increases with frequency. The pseudo-continuum opacity is due to the far wings of a multitude of pressure-broadened water vapour lines above 1 THz (Carilli & Holdaway, 1999).

In contrast to the other atmospheric chemical components, water vapour mixes poorly and its time-variable 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,


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 mm-VLBI 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 mm-VLBI: turbulence-driven 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 power-law 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,


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 Kramers-Kronig 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 line-shape, 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,


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.

Figure 4: Simulated mean opacity (black) and sky brightness temperature (red) at  GHz for three typical ground pressures and temperatures over a typical PWV range (Lane, 1998) which approximately represent the sites of SPT (dots), ALMA (squares) and SMA (triangles). The legend shows the estimated input ground (pressure, temperature) parameters for each site.

3.3.2 Turbulent troposphere

Visibility phase instability due to tropospheric turbulence is a fundamental limitation to producing high fidelity, science-quality maps with a mm-VLBI array (Thompson et al., 2001). The coherence time-scale is typically too rapid ( s) for fast switching calibration, so other calibration procedures (e.g. water vapour radiometry, paired antennas, and/or self-calibration) must be performed. Self-calibration 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, Kolomogorov-turbulent 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 field-of-view (FoV) of a global mm-VLBI array is typically FoV  mas or  m at a height of 2 km, which is roughly 7-8 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)


Assuming power-law turbulence and integrating yields,


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 time-series takes the form of a Gaussian random walk per antenna. At mm-wavelengths, the spectrum of water vapour is non-dispersive 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).

Figure 5: Simulation of the total delay (top) and the turbulent atmospheric delay (bottom) for SMA (blue) and ALMA (green) sites towards Sgr A. Ground pressures and temperatures are the same as Fig. 4, PWV above each station is set to  mm, and the zenith coherence time is set to  s for both stations. Note that all tropospheric parameters are, however, independently set for each station. The conversion from time delay to phase at 230 GHz is  rad  ps.

3.3.3 Limitations to high-fidelity image reconstruction

A primary objective of MeqSilhouette is to understand and constrain systematic errors in mm-VLBI 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 EHT-measured 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.

Figure 6: The effect of residual troposphere phase noise on uniformly weighted images of a point source observed for 12 hours at 230 GHz with 4 GHz bandwidth with the following array: SPT, ALMA, SMA, SMT, LMT and JCMT, assuming SEFDs from Lu et al. (2014) and an elevation limit of 15. For simplicity the weather parameters at each station were set to: coherence time  sec; PWV depth  mm; ground pressure  mb; ground temperature  K. Top left: Interferometric map with thermal noise only. Top right: Atmospheric attenuation and sky noise (due to non-zero opacity) with 1% of the turbulent phase noise added. Bottom left: As previous, but with 3% of turbulent phase contribution. Bottom right: As previous, but with 6% turbulent phase contribution. The fractional turbulent phase contributions are illustrative of the effect of fringe-fitting errors. Note the source attenuation and centroid shift that results.

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 ISM-induced closure phase variation is essential in order to make accurate inferences regarding asymmetric, event-horizon scale structure from EHT observations (e.g. Fish et al., 2016; Ortiz-Leó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 non-separable 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 mm-VLBI 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 turbulence-induced 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 cross-hairs) 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 fringe-fitting algorithms and strategies, through use of an automated calibration procedure and including the added complexity of a time-variable 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 mm-VLBI observations and software advances in the broader radio interferometry community, a mm-VLBI 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 mm-VLBI simulations will enhance the scientific opportunities possible with the EHT.


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 O. Smirnov’s research is supported by the South African Research Chairs Initiative of the Department of Science and Technology and National Research Foundation.


  • 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
  • Ortiz-León et al. (2016) Ortiz-Leó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, 25-29 July 2011.,
  • 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
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description