© ESO 2019

1. Introduction

Sagittarius A* (Sgr A*) is a strong radio source at the center of the Milky Way. Proper motion measurements of Sgr A* with respect to two extragalactic radio sources that are close in angular separation have confirmed that it is the dynamical center of the Galaxy (Reid et al. 1999; Reid & Brunthaler 2004). By monitoring stellar orbits at the Galactic center, it was found that at the position of the radio source there is an object with a mass of 4.3 ± 0.4 × 10⁶ M_⊙ at 8.3 ± 0.4 kpc from Earth (Ghez et al. 2008; Gillessen et al. 2009). From the monitoring data, it has been estimated that only up to ∼4 × 10⁵ M_⊙ within a region of 1 pc could be attributed to an extended mass distribution, leaving a supermassive black hole as the only physically feasible explanation. The radio emission is therefore believed to be produced by accretion onto and outflow from the black hole.

Since its discovery by Balick & Brown (1974), Sgr A* has been monitored frequently at various wavelengths. The broadband radio spectrum is flat-to-inverted up to the “submm bump” at ∼10¹² Hz, which is interpreted as a transition from optically thick to optically thin synchrotron emission (Falcke et al. 1998; Bower et al. 2015). The submm bump was explained as originating from a compact synchrotron-emitting region with a radius comparable to that of the event horizon of Sgr A*, which led to the prediction of the appearance of a roughly circular “shadow” of the event horizon surrounded by gravitationally lensed emission from an accretion flow at submm wavelengths (Falcke et al. 2000).

Very Long Baseline Interferometry (VLBI) observations of Sgr A* at centimeter wavelengths have confirmed that the apparent size of the radio source becomes smaller toward higher frequencies (Bower et al. 2004; Shen et al. 2005) as is expected for a scatter-dominated source. Due to this blurring by interstellar scattering, the source looks like a two-dimensional Gaussian distribution that increases in size with the square of the observing wavelength. The intrinsic source structure starts to contribute to the measured size at λ ≲ 6 cm, and dominates over the scattering effects in the millimeter regime (Bower et al. 2006; Doeleman et al. 2008). Hence, it is only at mm (and shorter) wavelengths that an image of the intrinsic structure of the central object showing the black hole shadow can be obtained.

The source Sgr A* has also been detected in the near-infrared and X-ray regimes (e.g., Yusef-Zadeh et al. 2006; Porquet et al. 2008; Do et al. 2009; Meyer et al. 2009), where it shows simultaneous flares on ∼ hour time scales. In the submillimeter regime, variability seems to occur on slightly longer time scales (Marrone et al. 2008; Dexter et al. 2014). With light crossing time arguments, the size of the flare emission region can be constrained to ≲10 Schwarzschild radii.

The apparent angular size of the black hole shadow of Sgr A* assuming zero spin is $2\sqrt{27}GM/c^{2}D\approx 53$μas, where G is Newton’s gravitational constant, M ≈ 4.3 × 10⁶ M_⊙ is the black hole mass, c is the speed of light, and D ≈ 8.3 kpc is the distance to the observer (Falcke et al. 2000; Johannsen et al. 2012). This makes Sgr A* the black hole with the largest angular size in the sky, and therefore the most promising candidate to image a black hole shadow. Another prime candidate is the black hole in the center of the giant elliptical galaxy M 87. This black hole is about 2000 times further away (Bird et al. 2010), but is also about 1000–1500 times more massive than Sgr A* (Gebhardt et al. 2011; Walsh et al. 2013), so that its angular size is comparable to that of Sgr A*.

At 1.3 mm, the black hole shadow of Sgr A* can be resolved with Earth-size baselines of ∼9 Gλ, yielding an angular resolution of ∼23 μas. Resolving the black hole shadow is the main aim of the Event Horizon Telescope (EHT), a VLBI array consisting of (sub)mm stations across the globe (Fish et al. 2013). Observations with the Combined Array for Research in Millimeter-wave Astronomy (CARMA) in California, the SubMillimeter Telescope (SMT) in Arizona, and the SubMillimeter Array (SMA) in Hawaii have resolved structure of Sgr A* and M 87 on event horizon scales (Doeleman et al. 2008, 2012; Fish et al. 2011; Akiyama et al. 2015; Johnson et al. 2015; Lu et al. 2018). With this number of stations, the uv-coverage is not sufficient to image the source, but the size of Sgr A* is determined to be ∼40 μas, indicating structure on scales smaller than the event horizon. The intrinsic source size is measured to be (120 ± 34)×(100 ± 18) μas at 3.5 mm (Shen et al. 2005; Ortiz-León et al. 2016; Issaoun et al. 2019), and $(354\pm 4)\times ({126}_{-41}^{+55})$μas at 7 mm (Bower et al. 2014).

From the measurement of non-zero closure phases, the sum of interferometric phases on a triangle of baselines, Fish et al. (2016) concluded that the source is asymmetric at 1.3 mm, which can be attributed to either the intrinsic source structure or scattering effects. Asymmetric structure due to scattering or instrinsic source morphology is also found in closure phase and amplitude measurements at 3.5 mm (Ortiz-León et al. 2016; Brinkerink et al. 2016, 2019). In April 2017, 1.3 mm observations of Sgr A* have been performed with eight stations as part of the EHT: the IRAM 30 m telescope on Pico Veleta in Spain, the Large Millimeter Telescope (LMT) in Mexico, the Atacama Large Millemeter Array (ALMA), the Atacama Pathfinder Experiment (APEX) telescope in Chile, the SMT in Arizona, the SMA and James Clerk Maxwell Telescope (JCMT) in Hawaii, and the South Pole Telescope (SPT). With the increased uv-coverage and sensitivity of, for example, ALMA and the LMT, these observations may be suitable for image reconstruction.

Imaging the black hole shadow could provide a strong-field test of general relativity as it predicts its size and shape (e.g., Bambi & Freese 2009; Johannsen & Psaltis 2010; Goddi et al. 2017; Psaltis 2018). Psaltis et al. (2015) show that Sgr A* is the optimal target for a general relativistic null hypothesis test because of strong constraints on the opening angle of one gravitational radius m = GM/Dc², which is known to within ∼4% from stellar monitoring observations (Ghez et al. 2008; Gillessen et al. 2009). Uncertainties in our prior knowledge of the effect of interstellar scattering, which is still severe at 230 GHz, pose limitations to the accuracy of shadow opening angle measurements (see also Sect. 3.2). Psaltis et al. (2015) infer that the general relativistic null hypothesis can in principle be tested down to the ∼10% level at 230 GHz. Mizuno et al. (2018) compare synthetic EHT observations of a general relativistic magnetohydrodynamics (GRMHD, see also Sect. 3) model with accretion onto a Kerr black hole to synthetic observations of a GRMHD model with accretion onto a dilaton black hole, the latter of which is taken as a representative solution of an alternative theory of gravity. They show that, with the observational setup of the 2017 EHT observations, it could be extremely difficult to distinguish between these two cases. Similarly, Olivares et al. (2018) show that it could also be difficult for the EHT to distinguish between a black hole and a boson star, although in that case the differences in source size are slightly larger.

Observations from the EHT could also shed light on the nature of the accretion flow, which may be dominated by emission from an accretion disk or relativistic jet (e.g., Falcke & Markoff 2000; Yuan et al. 2003; Fish et al. 2009; Dexter et al. 2010; Gold et al. 2017). Mościbrodzka et al. (2014) generate ray-traced GRMHD images for different electron temperature prescriptions leading to disk- or jet-dominated models, the latter providing a better fit to the radio spectrum of Sgr A*. Chan et al. (2015) argue that with EHT images, one may be able to distinguish disk-dominated from jet-dominated models, but that additional information would be needed to measure plasma and black hole properties within these models. Broderick et al. (2016) use EHT closure phases to constrain the dimensionless black hole spin parameter a_*, disfavoring values larger than ∼0.5. A value of 1 would correspond to a maximally spinning black hole. They also determine the black hole inclination and position angle to within approximately a few tens of degrees, all within the context of semi-analytic radiatively inefficient disk models. The model fits show consistency over several observation epochs spanning seven years.

Observations from the EHT could thus produce the first image of a black hole shadow and put constraints on different accretion flow models. However, due to interstellar scattering effects for Sgr A* and limited uv-coverage, it will likely be difficult to perform high-precision tests of general relativity and measure black hole and plasma parameters with high accuracy. Observations performed at substantially higher frequencies would be less affected by interstellar scattering and increase the image resolution, allowing for more precise tests of general relativity and accretion models. For example, a resolution of ≲10 μas would start to make it possible to visually distinguish between the Kerr and dilaton black hole shadows in Mizuno et al. (2018), and between a black hole and boson star in Olivares et al. (2018).

VLBI is not only carried out from the ground, but also from space. This allows one to observe with longer baselines and thus obtain a higher angular resolution. There are also no atmospheric corruptions for space-based antennas. The first Space VLBI (SVLBI) observations were done by Levy et al. (1986), who detected fringes for three active galactic nuclei (AGN) at 2.3 GHz on baselines of up to 1.4 Earth diameters between the Tracking and Data Relay Satellite System (TDRSS) and ground-based telescopes in Australia and Japan. The first dedicated SVLBI mission was VSOP (Hirabayashi et al. 1998, 2000), with an eight-meter antenna carried by the satellite HALCA, orbiting between 560 (perigee) and 21 000 km (apogee) above the Earth’s surface. It was operational between 1997 and 2003, imaging bright AGN and masers at 1.6 and 5.0 GHz with a network of ground-based telescopes.

The second and currently only operational SVLBI mission is RadioAstron, which has a ten-meter antenna carried by the Spektr-R spacecraft, operating at wavelengths between 1.3 and 92 cm (Kardashev et al. 2013). It has an orbital perigee altitude of about 10 000 km and apogee of about 350 000 km, making it the largest interferometer to date. At the maximum frequency of 22 GHz, the resolution achievable with RadioAstron (7 μas) is in principle high enough to resolve event horizon-scale structures of Sgr A*. However, at this frequency, the intrinsic structure of Sgr A* is blurred by interstellar scattering too severely to image the black hole shadow (see Sect. 3.2). Also, the high resolution is only achievable in one direction because of the highly elliptical orbit of Spektr-R. Studies of two-element SVLBI setups were performed in the iARISE project (Ulvestad 1999; Murphy et al. 2005) and studies for the Chinese space Millimeter-wavelength VLBI array (Hong et al. 2014; Ji Wu, priv. comm.).

There are some examples of space-based submillimeter telescopes operating at the high frequencies (up to ∼690 GHz) considered in this work. The ESA Herschel satellite had a 3.5 m dish, and its onboard spectrometer HIFI covered wavelengths between ∼0.16 and 0.6 mm (480–1250 and 1410–1910 GHz; de Graauw et al. 2010). The Swedish-led Odin satellite had a 1.1 m dish and operated between 486 and 580 GHz and at 119 GHz (Frisk et al. 2003). The ESA Planck satellite had a 1.6 by 1.9 m primary mirror and instruments sensitive to frequencies between 30 and 857 GHz (Tauber et al. 2010).

In this paper we have investigated the possible imaging capabilities of a new SVLBI system concept consisting of two or three satellites in polar or equatorial circular medium Earth orbits (MEOs). With the individual satellites in slightly different orbits, this configuration has the capability to image Sgr A* and other black holes (e.g., M 87) with a resolution that is an order of magnitude higher than the resolution that can be obtained from Earth. We performed simulated observations of Sgr A* with realistic source and system parameters in order to assess the expected image quality that could be obtained with this setup. The system concept is introduced in Sect. 2. Our source models are described in Sect. 3. In Sect. 4 we describe our simulated observations. The simulation results are presented in Sect. 5, and our conclusions and ideas for future directions are summarized in Sect. 6.

2. System setup

2.1. Antennas, orbits and uv-coverage

The SVLBI setup considered in this paper consists of two or three satellites orbiting Earth in circular MEOs at slightly different radii. We first consider the setup of the initial design study conducted by Martin-Neira et al. (2017) and Kudriashov et al. (2019) for the purpose of the Event Horizon Imager (EHI). They propose launching satellites equipped with approximately three-meter reflectors into MEOs with radii of around ∼14 000 km. The MEOs should be circular and either polar or equatorial for stability purposes. In contrast to VSOP and RadioAstron, the EHI concept is a pure space-space interferometer rather than a ground-space interferometer, observing at frequencies up to ∼690 GHz.

Observing at high frequencies is possible in space because there is no phase corruption or signal attenuation by the atmosphere. There are several reasons for increasing the observing frequency. Firstly, the angular resolution of the array increases with frequency as the baseline length measured in wavelengths increases. Secondly, the effect of interstellar scattering on the observed image (discussed in Sect. 3.2) will be considerably smaller. Also, the emitted radiation will originate from closer to the black hole, tracing the lensed photon ring more closely (e.g., Falcke et al. 1993; Mościbrodzka et al. 2009, see also Sect. 3). The latter effect also causes the image variability to be generally less profound at high frequencies, since it is confined to a smaller region that is more dominated by general relativistic effects.

Since there is a small difference between the orbital radii of the satellites, the inner satellite orbits slightly faster than the outer one, slowly increasing the distance between the two as they move from their initial positions on the line intersecting them and the Earth’s center. As the baseline between them constantly changes orientation as seen from a fixed source, the resulting uv-coverage will have the shape of a spiral, with a dense and isotropic sampling of the uv-plane, allowing for high-fidelity image reconstructions. The angular resolution of the reconstructed image is determined by the maximum baseline length, which is limited not only by the orbital radius, but also by the occultation of the required intersatellite link (ISL, Sect. 2.2) by the Earth. We illustrate the concept with an example in Fig. 1, which also shows a three-satellite configuration. Such a system would allow the use of closure phase, which could relax technical system requirements (Sect. 2.2). Also, a three-satellite system has a faster uv-plane filling rate, and measurements could continue with two baselines when one baseline is occulted by the Earth.

Fig. 1. Example of satellite positions and uv-coverage for circular MEOs at different time stamps (vertical) for a system consisting of two (Cols. 1 and 2) and three (Cols. 3 and 4) satellites. The orbital radii are 13 500 and 13 913 km for the two-satellite system, and for the three-satellite system a satellite is added at an orbital radius of 13 638 km, at one third of the distance between the two. In order to make the spiral structure visually clearer, the inner radii were set smaller than the ones used in our imaging simulations. In Col. 1, the satellites are shown in red and the corresponding baseline in orange. The orbits are shown in blue, and the Earth is represented by the black disk. In Col. 2, the red and orange points show the past uv-track for the two directions along the baselines. The current uv-coordinates are shown as larger dots. In Cols. 3 and 4, three baselines and uv-tracks are shown in corresponding colors. From the initial satellite positions on the line between them and the Earth’s center, the distance between them increases and the uv-spiral spreads outwards (upper panels) until the Earth atmosphere occults the line of sight between two satellites (middle panels). As the inner satellite catches up with an outer satellite, the spiral is traversed inwards (lower panels).

With full understanding that many of the SVLBI system parameters mentioned above are extremely challenging, we take them as given and as input into the analysis presented in this paper. Some engineering aspects are discussed in Sect. 2.2. We emphasize that the aim of this work is not to provide a technical justification for the concept, but merely to investigate imaging capabilities, which could serve as input for future engineering design studies.

2.2. Directions for future engineering system analysis

Martin-Neira et al. (2017) and Kudriashov et al. (2019) conduct a first design study addressing the engineering domain of the two-satellite EHI. Although this paper is not an engineering study, we summarize some of the technical aspects and challenges here. These details should be addressed in greater depth in future engineering studies.

The concept assumes that each satellite generates the local oscillator signal by combining a sufficiently stable reference produced on board with the one received over the ISL from the other satellite. It further assumes that the cross correlation of the data streams from both satellites is conducted onboard in real time within a delay window compatible with the on-the-fly relative positioning of a sufficient accuracy. Whether or not the onboard correlation and processing would not have to be prohibitively complicated in order to reach a sufficiently low rate for the data transfer to the ground, where the actual fringe fitting would take place should also be investigated. The interferometer will need to be phase stable and phase calibrated in such a way that long integrations can be performed without coherence loss. It may be necessary to perform an acceleration search in fringe fitting as is done for RadioAstron.

For real-time cross-correlation, the baseline vector and its time derivatives should be known with an accuracy that scales down with the observing wavelength. This defines the requirements for the state vector knowledge of the EHI satellites. Traditionally, in all Space VLBI systems implemented to date (TDRSS, VSOP, and RadioAstron), the correlator delay model is based on estimates of the spacecraft vector provided by ground-based orbit determination assets. In the EHI concept we consider a different approach in which the baseline vector and its time derivatives are obtained (measured) directly between the EHI spacecraft in real time. As demonstrated by RadioAstron, orbit determination measurements (not to be confused with orbit reconstruction) are provided by means of radio measurements at the level that correspond to ∼20–1500 observing wavelengths at 92–1.3 cm, respectively, and radial velocity (Doppler) measurements at the mm/s accuracy level (Zakhvatkin et al. 2014; Duev et al. 2015). Escalating similar requirements for the EHI submm wavelength range leads to the baseline measurement requirements at a precision level of ∼20 cm for the baseline vector and 0.01 mm s⁻¹ for its time derivative.

While these precision levels are challenging for an onboard real-time system, they are achieved by modern ground-based systems supporting interplanetary missions (e.g., Iess et al. 2009; Duev et al. 2016). One should expect further improvement of the real-time space-borne baseline measurement accuracy based on the relative position determination from the Global Navigation Satellite Systems (GNSS) satellites, which orbit further out at ∼20 000 km altitude. These may be able to determine position at centimeter precision (e.g., Allende-Alba & Montenbruck 2016; Park & Kim 2016; Jäggi et al. 2016). Also, ranging measurements with an accuracy down to ∼30 μm could possibly be performed with the ISL (Zech & Heine 2014; J. Perdigues, priv. comm.). Lessons could also be learned from the technological preparation of the LISA mission (Johann et al. 2008, and references therein). The European Space Agency is currently undertaking a study into the precise relative positioning of MEO satellites.

The idea currently explored in the EHI engineering study is to carry out the fringe-search among the measured correlations on the ground, using for this a finer orbit reconstruction to get coherent phases for visibilities with a low signal-to-noise ratio. While the satellites are close together, Sgr A* could be detected as a strong point source, which could be used as a starting point for fringe finding and orbit determination. As this is a new way of doing space VLBI, its feasibility still needs to be demonstrated.

Furthermore, a stable frequency reference will be needed to perform VLBI at 690 GHz. There are two active hydrogen masers onboard RadioAstron (Kardashev et al. 2013), but these have served as frequency references up to observing frequencies of 22 GHz. There is room for improvement of such instrumentation (Rodrigo et al. 2018).

The three-meter reflectors considered by Martin-Neira et al. (2017) and Kudriashov et al. (2019) would fit in a medium-sized space launcher such as a Soyuz fairing. ESA’s larger space launcher Ariane 6 will have a useable payload diameter of 4.6 m (Arianespace 2016), so monolithic reflectors up to this diameter can in principle be launched. Larger reflectors would have to be deployable, although a large (up to 25 m in diameter) monolithic reflector has been considered with side-mounting on a super-heavy launcher for the ESA International VLBI Satellite (IVS) study (Pilbratt 1991; B. Ye. Chertok 1989, priv. comm.).

The choice of orbits is limited by several factors. The maximum orbital radius, and hence the array resolution and filling speed of the uv-plane, is limited in particular by the visibility of GNSS satellites. Three GNSS satellites should be visible at the same time for a reliable real time position determination. Assuming two navigation antennas with a field of view of ±30°, the maximum orbital radius is then 13 913 km (Kudriashov et al. 2019). The minimum orbital radius is determined by the inner Van Allen radiation belt, which is confined to a radius of ∼12 400 km (e.g., Bakhtiyarov 2014). The maximum baseline length then ranges from 21 252 to 24 726 km (Kudriashov et al. 2019), corresponding to a resolution of to 4.2–3.6 μas at 690 GHz. All orbits in this range would thus allow the imaging of the black hole shadow with a resolution that is about an order of magnitude higher than the resolution that can be obtained from the Earth. The time it takes for the spiral to be completed depends on the radial separation of the satellite orbits. Placing one of the satellites at the maximum radius, the spiral is completed in one month for a radial separation of 20 km and in six months for a radial separation of 3 km.

Note that in practice measurements can most likely not be performed at all points in the uv-spiral due to functional constraints. For example, the Sun or Moon may be in the line of sight to the observed object or perturbing the measurements (a detailed discussion on similar constraints for the RadioAstron mission is presented by Gurvits et al. 1991). Also, the observations may have to be carried out over several months, with possible interruptions due to attitude control, command and communication, orbit determination and correction, and other operational activities. These effects should be given close attention in the project design. There would be no further geometric seasonal effects because polar (and equatorial) circular MEOs have no nodal precession: the line of nodes (i.e., the intersection of the orbital plane and the Earth’s equatorial plane) will always be perpendicular to the line of sight toward the observed object. The orbits are thus fixed against for example, the Galactic center.

2.3. Other SVLBI setups for imaging Sgr A*

Adding a third satellite to the MEO system would triple the number of baselines at each time, so that the uv-plane can be filled much faster. Such a system would also enable the use of closure phases, which are immune to station-based phase corruptions and useful for non-imaging analysis. However, adding a third satellite would increase the complexity and cost of the already challenging mission concept. The data would need to be exchanged and correlated for three baselines instead of one. Since there are no atmospheric corruptions in space and the local oscillator signal is exchanged between the satellites, the advantage of having closure phases and an increase in uv sampling speed may not weigh up against the increase in mission complexity and cost. The necessity for closure quantities will need to be assessed as the concept develops further and the expected satellite-based phase corruptions due to, for example, the orbit determination are better quantified.

The EHI MEO concept is not the only setup one could consider for an SVLBI mission. Instead of two or three satellites forming a space-space interferometer, one could launch one or multiple satellites and observe in conjunction with ground-based stations, similarly to RadioAstron. Using a high-sensitivity station like ALMA, one could significantly reduce the required integration time and track the time evolution of the source using dynamical imaging (Johnson et al. 2017; Bouman et al. 2017) from Low Earth Orbits (LEOs, Palumbo et al. 2018, 2019). The angular resolution could also be increased by launching satellites into MEOs or Geosynchronous Equatorial Orbits (GEOs; Fish et al. 2019). However, with this setup, observing at the high frequencies we consider here would be difficult, as the raw data from the satellites would have to be sent to the ground, and the ground data would still be affected by atmospheric attenuation and phase corruption, which is severe at high frequencies. Also, the uv-coverage would not be as dense and uniform as with two satellites in MEOs. An advantage of this method would be that it has been done successfully at low frequencies for RadioAstron, so that less-advanced technology would have to be developed than for the two- or three-satellite space-space interferometer.

Although we consider a space-ground system focusing on resolving source dynamics valuable for understanding black hole accretion, in this work we have focussed on space-space systems. These are more suitable for high-resolution static imaging.

3. Synthetic image generation

In this section, we describe the generation of the general relativistic hydrodynamics (GRMHD) simulations, and present ray-traced images. These images are used as input for our simulated observations in the follow-up sections, where we present reconstructions of these images under the assumption of different variants and parameters of the space-space interferometer concept outlined in Sect. 2.

3.1. Theoretical emission maps of Sgr A* at 230 and 690 GHz

Sgr A* is a quiescent galactic nucleus. Its bolometric luminosity is low in units of the Eddington luminosity: L_Bol/L_Edd  ∼  10⁻⁹. The mass accretion rate onto Sgr A* is estimated to be 10⁻⁹ ≲ Ṁ ≲ 10⁻⁷ M_⊙ yr⁻¹ (Bower et al. 2005, 2018; Marrone et al. 2007). The spectral energy distribution of Sgr A* can be fit with models of radiatively inefficient accretion flows (RIAFs; Yuan & Narayan 2014) that are coupled to a jet model (Falcke et al. 1993; Falcke & Markoff 2000; Yuan et al. 2002; Mościbrodzka et al. 2014).

Because in RIAFs there are few particle interactions and electron cooling is inefficient, the accretion flow becomes a two-temperature advection dominated accretion flow (Narayan et al. 1998), where ions and electrons are described by different temperatures, T_p and T_e. The behavior of the infalling magnetized plasma can be modeled numerically using GRMHD simulations. Starting from an initial RIAF-type plasma density and magnetic field configuration, these simulations solve the equations of magnetohydrodynamics within a specified spacetime metric (e.g., Gammie et al. 2003; Mościbrodzka et al. 2009; Porth et al. 2017).

The synchrotron emission and absorption are calculated from the physical GRMHD quantities such as the plasma density, magnetic field, and temperatures. The gas pressure in the simulations is dominated by protons, so T_p is computed from the simulations. Additional assumptions must be made for the electron temperature T_e. Radiative transfer equations with source and sink terms are then integrated along the geodesics running from the pixels of a virtual “camera” located far away from the source to each point in the simulation domain. This results in an image of the source as seen by a distant observer.

In this work, we have used models of Sgr A* presented in Mościbrodzka et al. (2014). They generated radiative transfer models based on the 3D GRMHD simulation b0-high from Shiokawa (2013). This GRMHD simulation starts with a Fishbone-Moncrief torus (Fishbone & Moncrief 1976) with inner radius 12 GM/c² and pressure maximum at 24 GM/c² in Keplerian orbit at the equator of a rotating supermassive black hole. The black hole spin parameter a_* is set to 0.94 and its mass to 4.5 × 10⁶ M_⊙. Mościbrodzka et al. (2014) consider various electron temperature models. As examples, we use their models 16, 24, 31, and 39 (see Table 1 in Mościbrodzka et al. 2014) as these models represent different electron heating scenarios within the same physical model of the accretion flow. In models 16 and 31, the electrons are heated primarily in the turbulent accretion disk and hence only the disk around the black hole is visible. In models 24 and 39, the electrons are hot in the magnetized jet outflow while the electrons in the disk are cooler, so that the jet is visible in the images. These different heating scenarios are motivated by more detailed collisionless plasma models (Ressler et al. 2015; Kawazura et al. 2019).

We present total intensity images time-averaged over 810 GM/c³ (about 5 h) in Fig. 2 (left column). All images were generated with the relativistic ray tracing radiative transfer code ibothros (Noble et al. 2007; Mościbrodzka & Gammie 2018). We assumed that the source is at a distance of 8.5 kpc from Earth. The inclination angle between the black hole spin axis and line of sight is 60° for models 16 and 24, and 30° for models 31 and 39. Images were generated with a field of view of 210.44 μas (corresponding to 40 × 40 GM/c² at the distance of the black hole) and a resolution of 256 × 256 pixels at the two frequencies of 230 GHz (EHT frequency) and 690 GHz that EHI aims to use. In this work, we have demonstrated the space VLBI array performance in reconstructing images of the black hole shadow in case of these four, quite distinct models of plasma around the black hole.

Fig. 2. Time-averaged GRMHD source models used for simulated observations of Sgr A*. Images with the note “scattered” were convolved with the scattering kernel from Bower et al. (2006). The total flux of the model is given in the bottom right corner of each image. The field of view is 210 μas for all images. Colors indicate brightness/pixel in mJy (square root scale).

The images show emission in a region close to the event horizon. Gravitational lensing causes emission originating close to the black hole to bend around it, leading to the appearance of a lensed photon ring and the shadow (Falcke et al. 2000) in the center. Doppler boosting of emission from plasma moving toward the observer causes the apparent asymmetry. Near the horizon, the emission pattern of the disk (models 16 and 31) and jet (models 24 and 39) models is similar: at 690 GHz especially, the reconstructed image will be dominated by general relativistic effects that are not strongly dependent on the exact nature of the accretion flow. At 230 GHz, the difference between the models is more pronounced. Due to strong lensing of the emission originating in a small region close to the black hole, the observed variability in the image plane is generally less profound at higher frequencies. At lower frequencies, larger moving structures that are less easily averaged out can be seen further out in the accretion flow, especially for the jet models.

The total flux density of Sgr A* at 690 GHz is variable on intra-day time scales (Dexter et al. 2014). In the simulations, we set the accretion rate such that the total flux density at 230 GHz is within 30% of the 2.4 Jy measured by Doeleman et al. (2008). As the 1.3 mm flux density varies as well (e.g., Fish et al. 2011; Bower et al. 2018), all models shown are considered to be reasonable representations of the expected brightness of Sgr A*.

3.2. Interstellar scattering

Scattering of radio waves due to electron density fluctuations in the interstellar medium between the Earth and the Galactic center causes phase fluctuations of the incoming plane wave. The effect of scattering can usually be treated as a random phase-changing screen described by a spatial structure function D_ϕ(x) ≡ ⟨[ϕ(x₀ + x)−ϕ(x₀)]²⟩, where ϕ is the change in phase and x is a transverse screen coordinate (Narayan & Goodman 1989; Goodman & Narayan 1989; Johnson & Gwinn 2015). Scattering occurs in two main regimes. Diffractive scattering is dominated by fluctuations on the phase coherence length, which is the length scale r₀ on the scattering screen corresponding to a change in phase of one radian, i.e., D_ϕ(r₀) ≡ 1. It is quenched for sources larger than r₀ and has therefore only been relevant in observations of extremely compact sources such as pulsars. Refractive scattering is dominated by fluctuations on the refractive scale r_R, which corresponds to the apparent size of a scattered point source. For Sgr A*, refractive scattering is expected to affect observations at frequencies up to about 2 THz (Johnson & Gwinn 2015).

The scattering screen is generally assumed to be frozen, with a transverse velocity v_⊥ with respect to Sgr A* as the only source of variations in time. Three averaging regimes introduced by Narayan & Goodman (1989) and Goodman & Narayan (1989) are important for interferometric imaging of a scattered source. In the snapshot regime, source and background noise are averaged for a single scattering realization. Observing a source that extends over scales larger than r₀ or time scales longer than t_dif = r₀/v_⊥ brings one into the average image regime, where only the refractive noise is relevant. In the average image regime, the source image contains spurious refractive substructure. This substructure is quenched, but not smoothed for a source size that exceeds the refractive scale (Johnson & Gwinn 2015). The ensemble average is the average over many realizations of the (refractive) scattering screen. Averaging scattering screen realizations over a time t ≫ t_R = r_R/v_⊥ brings one into the ensemble average regime. For Sgr A*, t_R is expected to be about a day at 230 GHz and about three hours at 690 GHz. In the ensemble average regime, the source appears as the unscattered source convolved with a scattering kernel that effectively blurs the image. The size of the kernel, which is Gaussian at least down to centimeter wavelengths, increases with the square of the observing wavelength λ. Bower et al. (2006) measured the full-width-at-half-maximum (FWHM) of the scattering kernel major axis to be 1.309 ± 0.015 mas cm⁻², the FWHM of the minor axis to be $0.{64}_{-0.05}^{+0.04}$ mas cm⁻², and the position angle to be ${78}_{-1.0}^{+0.8}$ degrees east of north. The right column of Fig. 2 shows our model images blurred with this scattering kernel. The blurring effect is significant at 230 GHz, but hardly visible at 690 GHz due to the λ² size law. Using a physically motivated scattering model including refractive effects, Johnson et al. (2018) infer that the scattering kernel at millimeter wavelenghts may be smaller than predicted by extrapolating the kernel from Bower et al. (2006).

4. Simulated observations

In this section, we outline the process of simulating observations of the Sgr A* source models described in Sect. 3.

4.1. uv-sampling

To calculate model visibilities, we used the eht-imaging¹ software (Chael et al. 2016). This software package calculates the observed complex visibilities corresponding to a given source model and uv-coverage. As an example, we adopted the orbital parameters consistent with the two-satellite MEO setup as discussed by Martin-Neira et al. (2017) and Kudriashov et al. (2019). The two satellites are in circular orbits with radii of 13 892 and 13 913 km, the latter of which is the maximum based on the requirement of having simultaneous visibility of at least three GNSS satellites (Sect. 2.2). With this setup, the maximum baseline length is 1.9 × 10¹⁰λ for 230 GHz and 5.7 × 10¹⁰λ for 690 GHz, corresponding to an angular resolution of 11 and 3.6 μas, respectively. We also performed simulations with a three-satellite system. In this case, the third satellite was placed at a radius of 1 899 km, which is at one third of the distance between the inner- and outermost satellites.

The orbital period of the satellites is ∼4.5 h. The completion time of the full spiral, starting with the satellites at their minimum distance and ending when the line of sight between them is occulted by the Earth, is set by the radial difference between the satellite orbits, which is 21 km in this case. The full spiral is then traversed in 29 days. The orbital plane is initially set perpendicular to the line of sight to the observed source in order to keep the simulations free from any preferential directions initially (we study different geometries in Sect. 5.2.3).

The integration time per measurement t_int should be set short enough to avoid image corruption by uv-smearing. If, within an integration time, a displacement is made in the uv-plane that corresponds to an angle on the sky that is smaller than the source size θ_source, the reconstructed image will be affected (Thompson et al. 2017; Palumbo et al. 2019). In the case of our SVLBI system, motions in the uv-plane are dominated by the azimuthal component since the spiral contains many (29 days/4.5 h = 155) loops. The uv-vector rotates fastest at the edge of the spiral, where the satellites are at their maximum separation. Here, the uv-separation per integration time is given by

$$\begin{array}{c}\hfill ds=\frac{2\pi B_{\max }{t}_{\mathrm{int}}}{P},\end{array}$$(1)

where B_max is the maximum baseline length and P is the orbital period. The uv-smearing limit on the integration time is then

$$\begin{array}{c}\hfill t_{\mathrm{int}}<\frac{P}{2\pi B_{\max }{\theta }_{\mathrm{source}}}·\end{array}$$(2)

With B_max = 5.7 × 10¹⁰ Gλ, P = 4.5 h, and θ_source = 150 μas (the important emission features in our model images are within this field of view), we get t_int <  62 s. Since this limit only holds when the satellites are at their maximum separation, the integration time can be made longer when they are closer together. When calculating our uv-spirals, we set a uv-distance-dependent integration time that is equal to half the limit from Eq. (2). The integration time is then well within the uv-smearing limit everywhere while a sufficient signal-to-noise ratio can be accumulated. At the shortest baselines, we set the maximum integration time to 454 s in order to limit the uv-arcs to ten degrees.

The total integration time required for imaging is at least one iteration of the spiral (29 days), which is much longer than both the expected source variability time scale and the expected refractive time scale. We comment on mitigating source and scattering variability in Sects. 5.2.4 and 4.3, respectively. Simulated observations involving a third satellite are presented in Sect. 5.3.

4.2. System noise calculation

Thermal noise can be characterized by a circular Gaussian in the visibility plane with zero mean and standard deviation σ. The value of σ can be calculated with the standard noise equation used in radio interferometry

$$\begin{array}{c}\hfill \sigma =\frac{1}{0.88}\sqrt{\frac{{\mathrm{SEFD}}_1{\mathrm{SEFD}}_2}{2\mathrm{\Delta }\nu t_{\mathrm{int}}}},\end{array}$$(3)

where Δν is the observing bandwidth and t_int is the integration time (Thompson et al. 2017). The factor 1/0.88 results from two-bit quantization losses. The system equivalent flux densities (SEFD_1, 2) of the antennas can be calculated as

$$\begin{array}{c}\hfill \mathrm{SEFD}=\frac{2k_{\mathrm{B}}{T}_{\mathrm{sys}}}{\eta A},\end{array}$$(4)

where k_B is Boltzmann’s constant, T_sys is the system temperature, and A = π(D/2)² is the area of an antenna with diameter D. In our simulations, η = η_apη_corη_clock includes the efficiencies of the aperture, correlator, and clock, respectively. The system parameters adopted in this work are consistent with the EHI setup considered by Martin-Neira et al. (2017) and Kudriashov et al. (2019). The parameters and resulting noise for the most important system setups considered in this paper are shown in Table 1. In addition to these, we have also performed simulations for a system consisting of three 4.0 m satellites (Sect. 5.3).

The system temperature is consistent with single side band (SSB) SIS receivers as installed in the Herschel HIFI instrument (de Graauw et al. 2008), with a 10 K antenna temperature. The 4.4 m antennas would fit in the Ariane 6 spacecraft (Arianespace 2016). In the case of a three-satellite system, three 4.0 m dishes would fit. The aperture, correlator, and clock efficiencies are the current baseline figures for the EHI design study. The estimate for η_clock is based on the coherence of the Space Hydrogen Maser of the Atomic Clock Ensemble in Space (ACES; Goujon et al. 2010), at a time scale of one second. Since the clock connection time over the ISL would be of order 10² ms for the longest EHI baselines, this estimate may be pessimistic. The bandwidth of the intersatellite link is 10 GHz, which is the sum of the bandwidths in two polarizations, assumed to be 5 GHz each.

For all simulations in this paper, we have assumed that the baseline vector is known exactly, and hence no phase corruptions due to uncertainties in the orbital model were introduced. In practice, the baseline vector and its time derivatives are only be known up to a certain precision. Future engineering studies should determine the accuracy that can be reached, so that the phase corruptions and their effect on the images can be modeled properly.

4.3. Deblurring

Interstellar scattering introduces variable refractive substructure as the scattering screen traverses in front of the source (Johnson & Gwinn 2015, see also Sect. 3.2). The average of different realizations of the scattering screen is the scattering kernel described in Sect. 3.2, which is convolved with the background source image and has a blurring effect. The refractive time scale of approximately one day at 230 GHz and three hours for 690 GHz is much shorter than the planned observation time (in the order of months) in which visibilities are averaged, so that the scattering effect in the reconstructed average image can be approximated by the ensemble average kernel. If the size and orientation of the scattering kernel are known (Bower et al. 2006), its effect may be mitigated (Fish et al. 2014). Calculating the Fourier transform of the convolution of two functions I(ξ, η) and G(ξ, η) is equivalent to pointwise multiplication of the Fourier transforms Ĩ(u,v) and $\stackrel{\sim }{G}(u,\mathit{v})$ of those functions:

$$\begin{array}{c}\hfill I(\xi ,\eta )\ast G(\xi ,\eta )\rightleftharpoons \stackrel{\sim }{I}(u,v)\stackrel{\sim }{G}(u,v),\end{array}$$(5)

where denotes convolution and ⇌ denotes a Fourier transform. The visibilities of a scattered source may thus be corrected by dividing them by the Fourier transform of the scattering kernel. Assuming that the scattering kernel is known, the remaining effect of a scattered source on the deblurred visibilities is that its signal-to-noise ratio is smaller than for an unscattered source. As the scattering kernel from Bower et al. (2006) is Gaussian, so is its Fourier transform, and the scattering causes the measured S/N to drop steeply as a function of baseline length. This effect is relevant for the setup discussed here as the array performance is sensitivity-limited and the baselines are long (although recent modeling by Johnson et al. (2018) suggests that the kernel may in fact be non-Gaussian and smaller at millimeter wavelengths, mitigating this effect). At 230 GHz, where the scattering is much stronger than at 690 GHz (Fig. 2), deblurring will become a problem at long baselines as the low-S/N visibilities are divided by small numbers, blowing up the errors. In order for the deblurring to work, we imposed an S/N cutoff of 4.5. In order to remove some remaining outliers, we impose the requirement that the error on the visibility in a grid cell should be smaller than the flux of the source on the shortest baseline. The S/N cutoff of 4.5 was determined empirically by selecting the value that gave the best image reconstruction quality (Sect. 4.4).

4.4. Image quality metric

In order to quantify the quality of our image reconstructions, we calculated the normalized root-mean-squared error (NRMSE). The NRMSE is defined as (Chael et al. 2016)

$$\begin{array}{c}\hfill \mathrm{NRMSE}=\sqrt{\frac{\sum _{i=1}^{n^2}{|I{\prime }_i-I_i|}^2}{\sum _{i=1}^{n^2}{\left|I_i\right|}^2}},\end{array}$$(6)

where I_i is the ith pixel of the n × n pixels model image (Fig. 2), and $I_i^{\prime }$ is the same for the reconstructed image. The NRMSE is thus a pixel-by-pixel comparison of two images: the lower the NRMSE, the more similar the two images are. Since in our case the model and reconstructed images do not have the same number of pixels, the reconstructed image was regridded to the pixel size of the model image and the two images were aligned before calculating the NRMSE. The NRMSE is only a coarse image comparison metric: it does not compare, for example, image smoothness or contrast, or the reconstruction of specific features.

5. Simulation results

In this section, we describe the outcome of the simulations for which the setup is described in the previous sections. We start with a two-satellite system, first reconstructing images using conventional VLBI techniques and then using the assumption that the system can be made to behave like a connected interferometer, which requires sharing the LO signal and obtaining an excellent orbit reconstrucion in post-processing. We then describe imaging techniques and results for a three-satellite system.

5.1. Two-satellites: conventional VLBI techniques

Figure 3 shows the expected signal-to-noise ratio (S/N) of the visibility at the spiral points for the scattered source models (Fig. 2) and system noise (Table 1) described above, assuming two 4.4 or 25 m reflectors. The S/N is highest at short baselines sampling the integrated large-scale source structure. The S/N drops more steeply toward longer baselines at 230 GHz than at 690 GHz due to blurring by interstellar scattering. The contours show that the region in uv-space where an S/N of greater than seven can be reached is confined to ∼10–20 Gλ for a 4.4 m reflector, whereas the proposed orbital and frequency setup allows for baselines up to ∼60 Gλ at 690 GHz. With a 25 m reflector, S/N of greater than seven can be reached near the maximum baseline length for most models.

Fig. 3. Signal-to-noise ratio (colored) for the spiral uv-points calculated with the system parameters in Table 1 (but with a 25 m reflector diameter in the bottom row) and scattered source models in Fig. 2. Contours indicate S/N values of 3, 5, 7, and 20.

We note that in practice the uv-coverage will not be circular, but elliptical due to the declination of the source (−29° for Sgr A*). The decrease of the structural source information that can be obtained as the baselines get shortened in one direction will depend on the shape and orientation of the source on the sky (see Sect. 5.2.3 for additional discussion).

In ground-based VLBI, one can use only the visibilities with sufficient S/N over a single integration time for fringe detection. The S/N-limit for fringe detection is typically set to seven, which is sufficient for a false detection rate of less than 0.01% in a search of 10⁶ values of delay and delay rate (Thompson et al. 2017). Figure 4 shows images for the four source models at 690 GHz that were reconstructed using only visibilities that fulfill this requirement within the set integration time of half the uv-smearing limit over a total observation time of one month, as indicated in Fig. 3. With 4.4 m dishes, the image resolution is considerably lower than with 25 m dishes. With 25 m dishes, conventional VLBI techniques could be used to reconstruct images with a nominal resolution of 4 μas within one month of observing time. However, launching 25 m dishes that have sufficient surface accuracy to observe at 690 GHz is certain to be extremely challenging, if not impossible within the next few decades.

Fig. 4. Image reconstructions of simulated observations of all 690 GHz models (left to right) with a MEO system with two 4.4 (top row) or 25 (bottom row) meter dishes, using only the data points that have S/N of greater than seven over a total observation time of one month. Images were reconstructed with MEM using eht-imaging (Chael et al. 2016). The NRMSE-values when comparing to the input models (Fig. 2) are shown in Table 2. The field of view is 210 μas for all images. Colors indicate brightness/pixel in mJy (square root scale).

5.2. Two satellites: two-stage correlation

The S/N detection threshold of seven may be reduced considerably depending on the system setup. Because there is no atmosphere in space, the system could be made to behave like a connected interferometer using a two-stage correlation scheme. In this setup, the two local oscillator signals should be shared (Sect. 2.2) and an (a posteriori) orbital reconstruction down to the subwavelength level should be obtained, using, for example, the intersatellite link for high-accuracy ranging measurements. After the onboard correlation has been performed and the data has been sent to the ground, the refined orbit reconstruction could then be used to expand the fringe search to longer solution intervals. Visibilities that have S/N of less than seven on short time intervals (set by the uv-smearing limit) may then still be detected because the fringe can be tracked over longer timescales.

If this behavior can be achieved, visibilities with low S/N will be coherent. Once the fringe has been detected and provided that the system is phase-stable over long timescales, multiple low-S/N visibilities can be averaged to obtain high-S/N visibilities that can be used for image reconstruction. The uv-spiral may be traversed for multiple iterations to build up S/N over time. If, for example, the uncertainties on the reconstructed baselines are too large, the low-S/N visibilities will be incoherent and averaging the visibilities will not yield robust higher-S/N data points.

Kudriashov et al. (2019) show that for the orbits considered in this paper, a baseline vector knowledge uncertainty (1-sigma 3D) of 0.1 mm leads to a directivity loss 3 dB. It is not yet clear whether it will be possible to achieve such accuracy within a reasonable budget. Another challenging task is to obtain a sufficiently accurate estimate of the baseline velocity vector and possibly its acceleration. Specific requirements for these parameters depend on the characteristics of the processing system (the size of the delay and delay-rate windows) and will be addressed in other studies.

The dense uv-coverage of the two-satellite system allows to divide the uv-plane into a square grid, and then to average all the spiral points that lie within the same grid cell. The resulting uniform uv-coverage allows for image reconstruction by simply taking the Fourier transform of the gridded visibilities. The size of the grid cells should be kept limited in order to keep a field of view that is large enough to image the entire source (the finer the grid cell spacing, the larger the field of view of the resulting image) and avoid uv-smearing. In our simulations, we set the grid size to be equal to the uv-distance that corresponds to the field of view of the model image. With our source model and observational parameters, this results in a grid of 39 × 39 pixels at 230 GHz, and 116 × 116 pixels at 690 GHz, with a grid cell size of 0.49 Gλ. With our baseline-dependent integration time (Sect. 4.1), each grid cell typically contains one or two measurements per observing epoch.

In the following subsections, we first present the simulation results for model 39 observed with the system described above while varying the total integration time (Sect. 5.2.1). We then compare results for different source models (Sect. 3). In these simulations, we observe the scattered time-averaged source models in Fig. 2, setting the orbital plane perpendicular to the line of sight to the source. We study the effects of source declination and time variability in Sects. 5.2.3 and 5.2.4, respectively.

5.2.1. Total integration time

Figure 5 shows the improvement in signal-to-noise ratio for the gridded visibilities of model 39 observed at 230 GHz with increasing total integration time with a 4.4 m dish. The assumed bandwidth is 10 GHz, which one would reach by combining all polarizations (Sect. 4.2). Figure 6 shows the Fourier transform (dirty image) of these visibilities. The visibilities were deblurred (Sect. 4.3) before performing the FFT. The same plots are shown in Figs. 7 and 8 for an observing frequency of 690 GHz. Deblurring was not applied here as the extrapolated major axis of the scattering kernel at 690 GHz is only 2 μas, which is a factor of two smaller than the maximum angular resolution that can be obtained with the investigated setup. NRMSE values (Sect. 4.4) comparing the reconstructed images to the original model are shown in Table 3. For all images, the pixel values in the reconstructed image were scaled such that the total flux matched the total flux of the input model before the NRMSE was calculated.

Fig. 5. S/N map of the gridded visibilities of model 39 (scattered) at 230 GHz after integrating for 1, 6, and 24 months (left to right), with a reflector diameter of 4.4 meters. Contours indicate the points with an S/N of 3, 5, 7, and 20.

Fig. 6. FFT of the gridded and deblurred visibilities of model 39 (scattered) at 230 GHz after integrating for 1, 6, and 24 months (left to right), with a reflector diameter of 4.4 meters. The field of view is 210 μas for all images. Colors indicate brightness/pixel in mJy (square root scale). The NRMSE-values when comparing to the input model (Fig. 2) are shown in Table 3.

Fig. 7. Same as Fig. 5, but for an observation frequency of 690 GHz.

Fig. 8. Same as Fig. 6, but for an observation frequency of 690 GHz.

Comparing the S/N plots after one month in Figs. 5 and 7 to Fig. 3, the S/N has increased due to the combining of polarizations and averaging visibilities that are in the same grid cell. The latter effect is marginal as the number of measurements in a single grid cell is of order unity in most grid cells. The S/N improvement as a function of total integration time is smaller for 230 GHz than for 690 GHz because blurring by interstellar scattering decreases the S/N on long baselines at this frequency. This is also reflected in the NRMSE values (Table 3), which show a plateau in image quality as the noise decreases at 230 GHz, but keep decreasing at 690 GHz.

At 230 GHz, the resolution of the reconstructed images is comparable to the resolution of the EHT (with visibilities that have S/N greater than seven, on baselines up to ∼8 Gλ after one month), despite the longer baselines available in the SVLBI setup. Due to the dense uv-coverage and long integration, robust image reconstructions can be obtained. Observations at 230 GHz could also be useful for (initial) fringe detection at 690 GHz, orbit reconstruction, and cross-comparison with EHT results.

5.2.2. Different source models

Figure 9 shows the image reconstructions for the different source models (Fig. 2) and frequencies for a six-month observation with 4.4 m reflectors. NRMSE values are shown in Table 4. For the low-inclination models (31 and 39), the black hole shadow can be traced more easily than for the high-inclination models. At 690 GHz, the apparent difference between the disk (models 16 and 31) and jet (models 24 and 39) image reconstructions is small, because the image morphology is dominated by general relativistic effects such as gravitational lensing and Doppler boosting, and the jet feature in the simulation is relatively faint. At 230 GHz however, extended structure associated with the jet can be seen more clearly in the reconstructions of the jet models, especially for model 24. 230 GHz observations may thus be more useful to discriminate between disk and jet models, but 690 GHz observations will allow for significantly sharper reconstructions of the black hole shadow.

Fig. 9. FFT of the gridded visibilities of all models at all frequencies, integrated for six months with a 4.4 m reflector. The NRMSE-values when comparing to the input models (Fig. 2) are shown in Table 4. The field of view is 210 μas for all images. Colors indicate brightness/pixel in mJy (square root scale).

For comparison, Fig. 10 shows image reconstructions of simulated EHT observations of all scattered models made with the eht-imaging software (Chael et al. 2016, 2018). The stations included are the same as for the April 2017 EHT observations (Sect. 1). The integration time per measurement was set to 30 s, with a measurement cadence of 300 s, observing for one day in total. The bandwidth was set to 4 GHz. No atmospheric or instrumental effects were included except for thermal noise. Images were reconstructed from the simulated visibility amplitudes and closure phases using a maximum entropy (MEM) algorithm with the Gull-Skilling entropy function (Chael et al. 2016). Closure phases were used for image reconstruction instead of visibility phases because atmospheric corruptions would severely corrupt visibility phases, while closure phases are immune to these. The images were blurred with a Gaussian kernel with a size of half the beam corresponding to the array resolution, in order to mitigate spurious superresolved structures. The NRMSE values (Table 4) are higher than the NRMSE values of the high-S/N 230 GHz simulations of the SVLBI array, indicating less similarity between the input models and reconstructions. Visually comparing the reconstructions, the SVLBI reconstructions are more robust in that they contain less spurious structure than most of the EHT reconstructions, especially for the low-inclination models 31 and 39. Comparing Fig. 10 to the bottom row of Fig. 9, observing Sgr A* for multiple months with two 4.4 m reflectors in space at 690 GHz could produce an image of the black hole shadow with a quality that significantly surpasses the image quality that can be expected for the EHT.

Fig. 10. Image reconstructions of simulated EHT 2017 observations of all models including ensemble-average scattering at 230 GHz. Images were reconstructed using a MEM algorithm implemented in eht-imaging (Chael et al. 2016, 2018). The NRMSE-values when comparing to the input models (Fig. 2) are shown in Table 4.

5.2.3. Source declination

In previous simulations, the orbital plane of the satellites was set perpendicular to the line of sight to the source, so that the uv-coverage had the shape of a circular spiral. This orientation was chosen in order to keep the simulations free from any preferential directions with respect to the source geometry. In practice the orbits are polar or equatorial, while the source is at a certain angle with respect to the orbital plane. The line of nodes will remain perpendicular to the line of sight, so that in case of a polar orbit the projected east-west baselines can be maximized (although a compromising orientation may be chosen if multiple sources are observed from the same orbit). The declination of Sgr A* is −29°, so that the angle between the line of sight and the orbital plane is 29° in case of an equatorial orbit, and 61° if the orbit is polar. The effect of source declination on the uv-coverage is twofold. Due to the projection of the orbital plane as seen from the source, the baselines will get shortened in the north-south direction by a factor sinα, where α is the angle between the orbital plane and the line of sight. Also, depending on α and the orbital radius, the satellites may traverse the Earth’s shadow during their orbit, so that source visibility is temporarily lost and gaps occur in the uv-plane. For an orbital radius R, this will occur if α is smaller than

$$\begin{array}{c}\hfill {\alpha }_{\text{crit}}=arcsin\left(\frac{R_{\text{E}}}{R}\right),\end{array}$$(7)

where R_E is the radius of the Earth. Using R = 13 892 km (as considered for the simulations above) and R_E = 6378 km, α_crit = 27.3°. Thus, in the case of Sgr A*, the satellites will not be in the Earth’s shadow for the orbital radius considered here, for both the polar and equatorial configurations. The only effect of the source declination on the uv-coverage will be the foreshortening of the baselines as described above.

We simulated observations of model 39 with a 4.4 m reflector and 24 months of integration considering the declination of Sgr A* in either polar or equatorial orbits. The model image was rotated on the sky by 90 degrees to maximize the effect of baseline shortening in the direction where the S/N is highest. Figure 11 shows the resulting beam pattern, S/N map, and dirty images. The beam pattern (left panel) is the FFT of the uv-coverage, assigning a value of one to the grid cells containing data points, and zeroes to the empty grid cells. The beam pattern is indeed more elongated for an equatorial orbit, and the baselines get shortened by a factor sin29° ∼ 0.5 in the v-direction. This factor is cos29° ∼ 0.87 for a polar orbit. Projected baseline shortening leads to a slight increase in S/N (middle panel) in the u-direction as the grid cells contain more points. The reconstructed image for an equatorial orbit (right panel) shows some artifacts due to the beam pattern, which may be taken out by using for example a CLEAN algorithm (Högbom 1974). For both the polar and equatorial orbits, the black hole shadow is still well visible in the dirty images. Because of the low S/N on long baselines and extremely high resolution, a foreshortening factor of two does not severly limit our ability to image the source. With a single orbital setup, it is therefore possible to observe sources in a wide angular range on the sky.

Fig. 11. Beam pattern (left), S/N (middle), and FFT (right) of the gridded visibilities of model 39 at a frequency of 690 GHz, rotated by 90 degrees on the sky and observed for 24 months with a 4.4 m reflector. Upper panels: simulations in case of polar satellite orbits, and lower panels: equatorial satellite orbits, taking the declination of Sgr A* (−29°) into account.

5.2.4. Source and scattering variability

GRMHD simulations of Sgr A* and VLBI data products simulated from these exhibit variability on short time scales (∼ minutes, e.g., Mościbrodzka et al. 2009; Dexter et al. 2010; Shiokawa 2013; Dexter & Fragile 2013; Kim et al. 2016; Medeiros et al. 2017; Roelofs et al. 2017) at mm-wavelengths due to orbital dynamics of the turbulent structure. Integrating observations for multiple months thus strongly violates the static source assumption in standard aperture synthesis imaging.

In case of the EHT, attempting to image a single day of observations of Sgr A* as a variable source with standard imaging methods will indeed lead to unsatisfactory results because of this violation: the measured visibilities correspond to different images at different uv-points (Lu et al. 2016). However, Lu et al. (2016) also show that one can make use of the linearity of the Fourier transform to reconstruct the average image of the time-variable source. Averaging multiple images is equivalent to averaging the corresponding visibilities in uv-space. The important features imprinted on the observed source by general relativistic effects, such as the size of the lensed photon ring and crescent shape caused by Doppler boosting, are continuously present in the image and will therefore remain prominent in the average image of the source. The turbulent substructure can then be averaged out if enough epochs are observed, provided that the variability indeed occurs on small spatial scales. In combination with additional methods such as normalizing the visibilities to the total flux of the source and applying a smoothing algorithm in the uv-domain, averaging eight days of observations before imaging leads to a reconstruction that is almost equally similar to the input model as observing the time-averaged source for one day (Lu et al. 2016). Furthermore, if sufficient uv-coverage is obtained on source variability time scales, dynamical imaging methods (Johnson et al. 2017; Bouman et al. 2017) could be used to reconstruct movies of the source and solve for a time-averaged image simultaneously.

In the case of our space VLBI system, the situation is similar to that of the EHT in that the total integration time is much longer than the variability time scale of the source. Dynamical imaging methods would likely not work here because there are only one or three baselines at each time. It is already necessary to observe multiple epochs and averaging visibilities when using two small dishes, in order to obtain an S/N that is sufficient for imaging on the longest baselines (Fig. 3). Hence, the method from Lu et al. (2016) described above could be used to mitigate source variability.

Simulations of source variability using GRMHD over time scales of months are not available. However, we demonstrate here that this method works, in principle, by simulating SVLBI observations of the 81 GRMHD movie frames of model 39 from Mościbrodzka et al. (2014) that were used to obtain the averaged image in Fig. 2. The frames were spaced by 10 GM/c³, corresponding to 221 seconds for Sgr A*, resulting in a total movie duration of five hours. To include the effect of refractive substructure rather than just the ensemble average scattering kernel, the movie was scattered with a a scattering screen traversing in front of the source using the stochastic optics module in eht-imaging (Johnson 2016; Johnson et al. 2018). The frames were observed with two 4.4 m reflectors, and the movie was repeated every time the last frame was reached. The resulting visibilities were gridded and averaged as described in Sect. 5.2.

Figure 12 shows reconstructed images after integrating for 1, 6, and 24 months. The average source structure showing the size and shape of the black hole shadow can in principle be recovered using this method. After one iteration of the spiral, the source structure is well visible already. Of course, the five-hour GRMHD movie may not be representative of the source variability over multiple months. Due to for example flaring activity, the source may undergo more radical and large-scale changes over this time period, and recovering the average quiescent structure may be more challenging. In future studies, the limitations of imaging a possibly more strongly varying source should be investigated more deeply.

Fig. 12. FFT of the gridded visibilities of model 39 observed as a movie with a 4.4 m reflector for 1, 6, and 24 months (left to right) at 690 GHz. The movie was scattered with a scattering screen traversing in front of the source (Johnson 2016; Johnson et al. 2018). NRMSE values from left to right are 0.75, 0.56 and 0.39.

5.3. Three-satellite system

As an alternative to building a system that behaves similarly to a connected interferometer, a third satellite could be added to the system so that closure phases can be formed. Closure phase is the phase of the bispectrum, which is the product of complex visibilities on a triangle of baselines (Jennison 1958; Rogers et al. 1974). Hence, closure phase is the sum of the individual phases on the triangle baselines. They are immune to station-based phase errors due to, for example, positioning offsets in the reconstructed orbital model, as these are canceled out when the phases are summed.

5.3.1. Static source

Figure 13 shows image reconstructions for a three-satellite system observing the time-averaged models, where the third satellite was put at a radius of 13 899 km, which is at one third of the distance between the inner- and outermost satellite. Apart from the reflector diameter, which was set to 4.0 m so that an Ariane 6 spacecraft could fit three, the noise parameters were kept equal to the two satellite system. The images were reconstructed with eht-imaging using the bispectrum accumulated over 1, 6, and 24 months. For the averaging of 6 or 24 iterations of the uv-spiral, we assumed a system for which the individual phases are corrupted (i.e., no connected interferometer-like behavior as outlined in Sect. 5.2). Hence, we did not average complex visibilities, but we averaged the bispectrum, which still has coherent (closure) phases. Thermal noise on closure phases is Gaussian down to an S/N of approximately three, where it starts to deviate (e.g., Rogers et al. 1995; Roelofs et al. 2017). It is still close to Gaussian (at the level of approximately a few percent) down to an S/N of one, which makes multi-epoch averaging of the bispectrum a viable method to accumulate S/N. This is also reflected in the images, which show an increasing amount of detailed structure as the number of observing epochs increases.

Fig. 13. Image reconstructions of simulated observations of all 690 GHz models (left to right) with a three-satellite MEO system using the bispectrum alone. The bispectrum was accumulated over 1, 6, and 24 months (top to bottom). Images were reconstructed using a MEM algorithm implemented in eht-imaging (Chael et al. 2016, 2018). NRMSE-values are shown in Table 5.

The image quality after one month of integration with three 4.0 m satellites is better than for the two 4.4 m satellites that only use visibilities with S/N >  7 (Fig. 4). NRMSE values for the three-satellite reconstructions (Table 5) are generally in between those for the two 4.4 and 25 m dishes using S/N >  7 visibilities (Table 2) after integrating for one month. Comparing these NRMSE values to the ones for images made with gridded visibilities is not reliable because the images were made in a different way (maximum entropy versus a simple FFT). The latter have a systematically higher NRMSE as the noise on the gridded visibilities is transferred to the image plane, while the MEM algorithm fits a model image to the data, resulting in an artificially high dynamic range. As the total integration time increases, the three-satellite images visually do not become quite as sharp as the images made with complex visibility gridding (Fig. 9). Possible contributing causes are a higher noise level due to smaller dishes, systematics caused by averaging of data points with S/N of less than one, where the error distribution starts to significantly deviate from a Gaussian distribution, and the fact that less information was used for image reconstruction (bispectrum vs. complex visibilities).

5.3.2. Variable source

Figure 14 shows the same as Fig. 13, but for the GRMHD model 39 observed as a movie instead of a time-averaged image. Here, the image quality is significantly worse than for the two-satellite observation of the time-variable source (Fig. 12), with more spurious-substructure. The reason for this difference is the fact that the average of a set bispectra does not correspond to the Fourier transform of the average of the set of images corresponding to those bispectra, which is also noted by Lu et al. (2016). The relation holds for complex visibilities, but not for the triple product of these. So, for a time-variable source, the bispectrum alone is not sufficient to reconstruct a static image with a quality similar to the two-satellite system employing a two-stage correlation scheme. Either a combination of these techniques should be used, or more advanced (dynamical) imaging techniques should be developed for this purpose.

Fig. 14. Same as Fig. 13 but for model 39 at 690 GHz observed as a movie. NRMSE-values are 0.39, 0.52, and 0.45 from left to right.

6. Summary and outlook

In this paper, we have presented imaging simulations of the EHI SVLBI system consisting of two MEO satellites in circular orbits at slightly different radii, as discussed by Martin-Neira et al. (2017) and Kudriashov et al. (2019). The EHI could be used to image the black hole shadow of Sgr A* up to frequencies of about 690 GHz. Such high observing frequencies can be reached in space because of the absence of atmospheric corruptions. The setup allows for long baselines (up to ∼60 Gλ at 690 GHz), resulting in a maximum image resolution of 4 μas, which is a significant improvement on the ∼23 μas resolution that can be obtained with EHT baselines at 230 GHz. The two-element interferometer setup results in a spiral-shaped sampling of the uv-plane with a density that cannot be obtained with Earth-based VLBI, so that high-fidelity images can be reconstructed. Apart from the higher resolution, advantages of observing at higher frequencies are the small interstellar scattering and source variability effects at 690 GHz compared to 230 GHz, and the closer origin of the emission to the event horizon.

Using GRMHD simulations of Sgr A* and model system parameters, we have performed simulated observations in order to assess the image quality that can be expected. The signal-to-noise ratio of the measured visibilities is expected to be less than seven on baselines longer than 10–20 Gλ, preventing robust fringe detection on these baselines using conventional VLBI methods. However, the detection threshold may be decreased by using a system with excellent clock and orbit reconstruction (≲0.1 mm) accuracy. Higher-S/N measurements may then be obtained by averaging visibilities measured in different iterations of the uv-spiral. If such a system cannot be built within a reasonable budget, one would need to launch two 25 m antennas rather than 4.4 m antennas, in order to obtain sufficent S/N for conventional fringe fitting on long baselines.

At 230 GHz, the expected image resolution is comparable to the expected resolution of the images produced by the EHT because of stronger scattering effects on long baselines, although the reconstructed SVLBI images are more robust due to the dense and uniform uv-coverage. At 690 GHz, interstellar scattering has only a small effect on the observed image, and the proposed setup could allow for reconstructed images of Sgr A* with unprecedented angular resolution and fidelity within one or a few months of integration.

We show that source variability can be averaged out to reconstruct an image of the quiescent source structure showing the photon ring and Doppler-boosted emission. We note that the ability to reconstruct an average image from a time-variable source using this method depends on the nature of the variability. If the variability is caused by small-scale turbulent structures while the large-scale features remain prominent, such as in the GRMHD simulations we have considered, variability can indeed be averaged out. If, on the other hand, there are large-scale structural changes in the source, this will become more difficult. Since Sgr A* is a variable source, future studies leading to a full mission proposal should further investigate the ability to reconstruct an image under the assumption of different variability scenarios within the parameter space allowed by existing (EHT) measurements.

If the phase stability and orbit reconstruction accuracy of the two-satellite system are not sufficient to obtain detections on long baselines, three four-meter antennas could be launched so that closure phases could be formed and used for imaging. Since closure phases are immune to station-based phase errors, such a system could relax the orbit reconstruction and stability requirements. However, a system solely relying on measurements of the bispectrum poses challenges for imaging a time-variable source.

There are still significant technical challenges to overcome for the concept to be turned into an actual mission. The main issues to be resolved are the maximum orbit reconstruction accuracy that can be obtained, the complexity of the onboard correlation and processing that would be needed to send reduced data to the ground, and the frequency reference stability for 690 GHz observations. These challenges should be addressed in future engineering studies. More investigations should also be made into the possibilities of reducing the system noise as this is an important determining factor for the image quality.

The EHI concept could be of great astrophysical interest as it allows for precise tests of general relativity and accretion models. A quantitative comparison of the precision of these tests between the EHT and the SVLBI experiments discussed here, with inclusion of all instrumental corruptions, should be done as the project develops further. Furthermore, observations of GRMHD data that are ray-traced in full Stokes (e.g., Gold et al. 2017; Mościbrodzka & Gammie 2018) could be simulated, to infer what could be learned about the magnetic field structure near the event horizon.

Apart from Sgr A*, other sources will be interesting to observe with the SVLBI concept presented here as well. Emission from M 87, the black hole with the second largest apparent size on the sky and also a prime EHT target, is not affected by interstellar scattering. Imaging it at 230 GHz with the long baselines of the MEO SVLBI experiment could thus have a significant advantage over imaging it from the ground at the same frequency. Another advantage of imaging M 87 is that it is variable at ∼10³ times longer time scales than Sgr A*, possibly allowing for static snapshot reconstructions and multi-epoch dynamical reconstructions depending on the satellite separation which sets the radial uv-filling speed. Since GRMHD simulations of M 87 exhibit similar features as GRMHD simulations of Sgr A* at millimeter wavelengths (e.g., Mościbrodzka et al. 2016), the static imaging results presented here for Sgr A* may be largely applicable to M 87 as well, provided that its mass is close to the 6.6 × 10⁹ M_⊙ measured by Gebhardt et al. (2011).

The Sobrero Galaxy M 104 hosts a supermassive black hole of ∼10⁹ M_⊙ (Kormendy et al. 1996) at a distance of 9.55 ± 0.13 ± 0.31 Mpc (McQuinn et al. 2016), which yields an apparent event horizon size of ∼11 μas, which can be resolved by EHI baselines. The black hole at the center of the elliptical galaxy M 84 has a mass of $8.5_{-0.8}^{+0.9}\times {10}^8{M}_{\odot }$ (Walsh et al. 2010). At a distance of 17 Mpc, the size of the event horizon is ∼5 μas on the sky, which is comparable to the EHI resolution at 690 GHz. M 81* has a black hole with mass $7.0_{-1}^{+2}\times {10}^7{M}_{\odot }$ (Devereux et al. 2003) at a distance of 3.63 ± 0.34 Mpc (Freedman et al. 1994), yielding an apparent event horizon diameter of ∼2 μas. Another close active galactic nucleus is Centaurus A, at a distance of 3.8 ± 0.1 Mpc (Harris et al. 2010) and with a black hole mass of (5.5 ± 3)×10⁷ M_⊙ (Neumayer et al. 2010). For this black hole, the event horizon (∼1 μas) may be too small to resolve with the setup discussed here, but at this distance the 4 μas angular resolution at 690 GHz corresponds to a linear scale of only two light hours. This would enable one to image the structure of the relativistic jet on a length scale that is two orders of magnitude shorter than that already achieved (Müller et al. 2014). Similarly, jets of several other AGN could be studied in detail in order to improve our understanding of jet launching and collimation.

Different variations of the presented concept could be explored. Depending on the technical possibilities, one could try to push for even higher frequencies, which would increase the resolution further. A shorter separation of the orbits might enable studying various objects, such as protoplanetary disks (e.g., Hogerheijde et al. 2011), at lower resolution for many orbits before the satellites separate. Further studies should assess whether this could lead to valuable science. Another possibility is investigating a space-space-ground hybrid system that can perform both lower-frequency space-ground observations for dynamical imaging (Palumbo et al. 2019) and higher-frequency space-space observations for high-resolution static imaging.

Acknowledgments

This work is supported by the ERC Synergy Grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (Grant 610058). We thank Andrew Chael and Katie Bouman for making the eht-imaging software publicly available, and for providing the code for regridding and aligning the images and calculating the NRMSE. The design study was the result of a colloquium on space interferometry by H. Falcke at ESTEC on Sep 2, 2015 and an internal ESA study organized by M. Martin-Neira. We are grateful to the anonymous referee for useful and constructive comments. M. Mościbrodzka acknowledges H. Shiokawa for providing the HARM3D GRMHD simulation data used in our ray-tracing simulations. We thank Alan Roy, Michael Bremer, Jason Dexter, Itziar Barat, Thijs de Graauw, Vincent Fish, Andrey Baryshev, André Young, and Daniel Palumbo for useful comments and discussions on this work.

References

All Tables

All Figures

Fig. 1. Example of satellite positions and uv-coverage for circular MEOs at different time stamps (vertical) for a system consisting of two (Cols. 1 and 2) and three (Cols. 3 and 4) satellites. The orbital radii are 13 500 and 13 913 km for the two-satellite system, and for the three-satellite system a satellite is added at an orbital radius of 13 638 km, at one third of the distance between the two. In order to make the spiral structure visually clearer, the inner radii were set smaller than the ones used in our imaging simulations. In Col. 1, the satellites are shown in red and the corresponding baseline in orange. The orbits are shown in blue, and the Earth is represented by the black disk. In Col. 2, the red and orange points show the past uv-track for the two directions along the baselines. The current uv-coordinates are shown as larger dots. In Cols. 3 and 4, three baselines and uv-tracks are shown in corresponding colors. From the initial satellite positions on the line between them and the Earth’s center, the distance between them increases and the uv-spiral spreads outwards (upper panels) until the Earth atmosphere occults the line of sight between two satellites (middle panels). As the inner satellite catches up with an outer satellite, the spiral is traversed inwards (lower panels).

In the text

	Fig. 2. Time-averaged GRMHD source models used for simulated observations of Sgr A*. Images with the note “scattered” were convolved with the scattering kernel from Bower et al. (2006). The total flux of the model is given in the bottom right corner of each image. The field of view is 210 μas for all images. Colors indicate brightness/pixel in mJy (square root scale).
In the text

	Fig. 3. Signal-to-noise ratio (colored) for the spiral uv-points calculated with the system parameters in Table 1 (but with a 25 m reflector diameter in the bottom row) and scattered source models in Fig. 2. Contours indicate S/N values of 3, 5, 7, and 20.
In the text

Fig. 4. Image reconstructions of simulated observations of all 690 GHz models (left to right) with a MEO system with two 4.4 (top row) or 25 (bottom row) meter dishes, using only the data points that have S/N of greater than seven over a total observation time of one month. Images were reconstructed with MEM using eht-imaging (Chael et al. 2016). The NRMSE-values when comparing to the input models (Fig. 2) are shown in Table 2. The field of view is 210 μas for all images. Colors indicate brightness/pixel in mJy (square root scale).

In the text

	Fig. 5. S/N map of the gridded visibilities of model 39 (scattered) at 230 GHz after integrating for 1, 6, and 24 months (left to right), with a reflector diameter of 4.4 meters. Contours indicate the points with an S/N of 3, 5, 7, and 20.
In the text

	Fig. 6. FFT of the gridded and deblurred visibilities of model 39 (scattered) at 230 GHz after integrating for 1, 6, and 24 months (left to right), with a reflector diameter of 4.4 meters. The field of view is 210 μas for all images. Colors indicate brightness/pixel in mJy (square root scale). The NRMSE-values when comparing to the input model (Fig. 2) are shown in Table 3.
In the text

	Fig. 7. Same as Fig. 5, but for an observation frequency of 690 GHz.
In the text

	Fig. 8. Same as Fig. 6, but for an observation frequency of 690 GHz.
In the text

	Fig. 9. FFT of the gridded visibilities of all models at all frequencies, integrated for six months with a 4.4 m reflector. The NRMSE-values when comparing to the input models (Fig. 2) are shown in Table 4. The field of view is 210 μas for all images. Colors indicate brightness/pixel in mJy (square root scale).
In the text

	Fig. 10. Image reconstructions of simulated EHT 2017 observations of all models including ensemble-average scattering at 230 GHz. Images were reconstructed using a MEM algorithm implemented in eht-imaging (Chael et al. 2016, 2018). The NRMSE-values when comparing to the input models (Fig. 2) are shown in Table 4.
In the text

	Fig. 11. Beam pattern (left), S/N (middle), and FFT (right) of the gridded visibilities of model 39 at a frequency of 690 GHz, rotated by 90 degrees on the sky and observed for 24 months with a 4.4 m reflector. Upper panels: simulations in case of polar satellite orbits, and lower panels: equatorial satellite orbits, taking the declination of Sgr A* (−29°) into account.
In the text

	Fig. 12. FFT of the gridded visibilities of model 39 observed as a movie with a 4.4 m reflector for 1, 6, and 24 months (left to right) at 690 GHz. The movie was scattered with a scattering screen traversing in front of the source (Johnson 2016; Johnson et al. 2018). NRMSE values from left to right are 0.75, 0.56 and 0.39.
In the text

	Fig. 13. Image reconstructions of simulated observations of all 690 GHz models (left to right) with a three-satellite MEO system using the bispectrum alone. The bispectrum was accumulated over 1, 6, and 24 months (top to bottom). Images were reconstructed using a MEM algorithm implemented in eht-imaging (Chael et al. 2016, 2018). NRMSE-values are shown in Table 5.
In the text

	Fig. 14. Same as Fig. 13 but for model 39 at 690 GHz observed as a movie. NRMSE-values are 0.39, 0.52, and 0.45 from left to right.
In the text