11institutetext: Department of Astrophysics, IMAPP, Radboud University, PO Box
9010, 6500 GL Nijmegen, The Netherlands 22institutetext: Instituto de
Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile,
Casilla 306, Santiago 22, Chile 33institutetext: Center for Astrophysics $|$
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA
44institutetext: Black Hole Initiative, Harvard University, 20 Garden Street,
Cambridge, MA 02138, USA 55institutetext: NASA Hubble Fellowship Program,
Einstein Fellow 66institutetext: Center for Computational Astrophysics,
Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA 77institutetext:
Department of Astronomy and Columbia Astrophysics Laboratory, Columbia
University, 550 W 120th St, New York, NY 10027, USA 88institutetext: Max-
Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany
99institutetext: ASTRON, The Netherlands Institute for Radio Astronomy,
Postbus 2, NL-7990 AA Dwingeloo, The Netherlands

# Imaging the event horizon of M87* from space on different timescales

A. Shlentsova 1122 F. Roelofs 334411 S. Issaoun 335511 J. Davelaar 667711 H.
Falcke 118899

(Received June 16, 2023; accepted February 4, 2024)

###### Abstract

Context. The concept of a new space very long baseline interferometry (SVLBI)
system named the Event Horizon Imager (EHI) has been proposed to dramatically
improve black hole imaging and provide precise tests of the theory of general
relativity.

Aims. This paper presents imaging simulations for the EHI. We investigate the
ability to make high-resolution movies of the black hole shadow and jet
launching region around the supermassive black hole M87* and other black hole
jets with a three-satellite EHI configuration. We aim to identify orbital
configurations to optimize the $uv$-coverage to image variable sources.

Methods. Observations of general relativistic magnetohydrodynamics (GRMHD)
models were simulated for the configuration, consisting of three satellites in
circular medium earth orbits with an orbital plane perpendicular to the line
of sight. The expected noise was based on preliminary system parameters. Movie
frames, for which a part of the $uv$-coverage may be excessively sparse, were
reconstructed with algorithms that recover missing information from other
frames. Averaging visibilities accumulated over multiple epochs of
observations with an appropriate orbital configuration then improves the image
quality. With an enhanced signal-to-noise ratio, timescales of observed
variability were decreased.

Results. Our simulations show that the EHI with standard system parameters is
capable of imaging the variability in the M87* environment on event horizon
scales with approximately a month-long temporal resolution. The EHI with more
optimistic noise parameters (enhancing the signal-to-noise ratio about
100-fold) would allow for imaging of the variability on gravitational
timescales. Observations with an EHI setup at lower frequencies are capable of
imaging the variability in extended jets.

Conclusions. Our study shows that the EHI concept can be used to image the
variability in a black hole environment and extended jets, allowing for
stronger tests of gravity theories and models of black hole accretion, plasma
dynamics, and jet launching.

###### Key Words.:

Galaxy: center – Techniques: interferometric – Techniques: high angular
resolution – Methods: data analysis

##  1 Introduction

The supermassive black hole in the nucleus of the elliptical galaxy M87 (M87*)
has been an object of great interest for intensive studies with imaging
observations since its discovery in 1978 (Young et al., 1978; Sargent et al.,
1978). Measurements of the black hole mass $M_{BH}=(6.5\pm
0.2{|}_{\mathrm{stat}}\pm 0.7{|}_{\mathrm{sys}})\times 10^{9}M_{\odot}$
extracted from the direct measurements of the angular diameter of the shadow
(Event Horizon Telescope Collaboration et al., 2019a, b, c, d, e, f) are
consistent with the presence of a central Kerr black hole. The estimated mass
from the Event Horizon Telescope (EHT) observations agrees with stellar-
dynamic observations for which the most recently obtained result is
$M_{BH}=(6.6\pm 0.4)\times 10^{9}M_{\odot}$ (Gebhardt et al., 2011). Falcke et
al. (2000) introduced the black hole ‘shadow’ — the region on the sky plane of
a black hole in which there is a noticeable deficit of the observed intensity
due to gravitational lensing effects. They proposed that this shadow is
resolvable with very long baseline interferometry (VLBI) at sub-millimetre
wavelengths. The shadow is caused by two effects: the strong gravitational
redshift and a shorter total path length of photons escaping to the observer
from geodesics intersecting the horizon, while photons on geodesics missing
the horizon can orbit the black hole near the circular radius several times,
which leads to a higher integrated emissivity (Bronzwaer & Falcke, 2021). This
region is a circle of radius $\sqrt{27}R_{g}$ (where $R_{g}=GM_{BH}/c^{2}$,
$G$ is the gravitational constant, $M_{BH}$ is the mass of the black hole and
$c$ is the velocity of light) in the Schwarzschild case (the dimensionless
black hole spin parameter $a_{*}=0$) and has a more flattened shape of a
similar size for a Kerr black hole. Numerical studies on the observability of
the M87* shadow were either done using semi-analytical models (Broderick &
Loeb, 2009) or with general relativistic magnetohydrodynamics (GRMHD)
simulations (Dexter et al., 2012; Mościbrodzka et al., 2016, 2017; Ryan et
al., 2018; Chael et al., 2019b; Davelaar et al., 2019). The shadow is directly
connected to the time-like component of the space-time metric that is also
probed by gravitational wave measurements (Psaltis et al., 2020, 2021).
M87* has a black hole shadow with an angular size of $\sim
42\leavevmode\nobreak\ \mu$as (Event Horizon Telescope Collaboration et al.,
2019a). It is the second-largest black hole shadow on the sky as seen from
Earth, after the Galactic Centre black hole Sagittarius A* (Sgr A*), the
supermassive black hole at the centre of the Milky Way. Sgr A* has a black
hole mass of $M_{BH}=(4.152\pm 0.014)\times 10^{6}{M}_{\odot}$ and an angular
shadow size of $\sim 52\leavevmode\nobreak\ \mu$as (Doeleman et al., 2008;
Johnson et al., 2015; Lu et al., 2018; Abuter et al., 2019; Event Horizon
Telescope Collaboration et al., 2022a, b, c, d, e, f). M87* has several
advantages for observations in comparison with Sgr A*, namely the absence of
scattering by the interstellar medium between us and the central radio source.
The timescales of structural variations of a black hole are characterized by
the gravitational time $t_{\mathrm{G}}=GM_{BH}/c^{3}$ and the innermost stable
circular orbital time $t_{\mathrm{ISCO}}\approx 24.25t_{\mathrm{G}}$ in the
case of the rapidly spinning Kerr black hole ($a_{*}=0.94$). The timescales of
structural variations of M87* are, therefore, $\sim 10^{3}$ longer than those
of Sgr A*, allowing us to assume a static source during single-day
observations and make use of classical aperture synthesis techniques. Another
striking difference is that M87* hosts a powerful relativistic jet observed in
the radio, optical, and X-ray bands (e.g. Sparks et al., 1996; Marshall et
al., 2002). Although Sgr A* may also host a jet, it has not been clearly
observed yet. Altogether this makes M87* an auspicious source for observations
with VLBI at millimetre and sub-millimetre wavelengths.
Multiple VLBI imaging studies of M87* have been carried out at 3.5 mm (a
frequency of 86 GHz), 7 mm (43 GHz), 1.3 cm (24 GHz), and longer wavelengths
(e.g. Hada et al., 2016; Kravchenko et al., 2020; Chang et al., 2010). The
development of high-bandwidth VLBI systems made it possible to decrease
observational wavelength down to 1.3 mm (230 GHz) for the EHT and reach an
instrument angular resolution of about 25 $\mu$as (Event Horizon Telescope
Collaboration et al., 2019b). Therefore, the EHT array can resolve the shadow
of M87* (Event Horizon Telescope Collaboration et al., 2019a, b, c, d, e, f).
The EHT data provided the first image of the black hole shadow with the
surrounding ring of emission. This allowed for the exclusion of some
categories of black hole models (Event Horizon Telescope Collaboration et al.,
2019e). With the addition of linear- and circular-polarimetric images of the
ring, magnetically arrested disk (MAD) models (Narayan et al., 2003) with
$a_{*}$ roughly around -0.5 and 0.0 were found to be favoured (Event Horizon
Telescope Collaboration et al., 2021a, b; Event Horizon Telescope
Collaboration et al., 2023). Nevertheless, conclusions are still strongly
dependent on the imaging assumptions. A further reduction of acceptable models
requires more accurate spin measurements as well as an improved estimation of
the magnetization and electron temperature distribution. Quantitative spin
measurements require higher-resolution imaging of the source (Van der Gucht et
al., 2020; Roelofs et al., 2021). Measurements of the photon subring
parameters of M87* could offer precise measurements of the black hole mass and
spin (Johnson et al., 2020). In addition, more accurate measurements of the
static photon ring allow for deeper tests of the Kerr metric and consequently
the theory of general relativity (GR) (Vagnozzi et al., 2022). Moreover, the
plasma behaviour in a black hole environment is currently not fully
understood. Imaging the variability of the emitting plasma would provide
constraints on the plasma parameters. Hence imaging variability in the M87*
environment with a higher resolution and fidelity would deliver data for more
accurate reconstructions of the plasma surrounding the black hole, as well as
GR tests. Furthermore, such measurements could be used to test different
theories of gravity (Mizuno et al., 2018) and other non-Kerr alternatives such
as boson stars (Olivares et al., 2020), axion models (Chen et al., 2020), or
fuzzballs (Bacchini et al., 2021). However, this requires improvements to be
made to the EHT.
Earth-based extensions of the EHT array, such as the next-generation EHT
(ngEHT, Doeleman et al., 2019; Raymond et al., 2021; Johnson et al., 2023;
Doeleman et al., 2023), will substantially improve imaging with increased
sensitivity and the additional VLBI baseline coverage from new sites. The
ngEHT will also have an expanded frequency coverage up to 345 GHz.
Nevertheless, angular resolution improvements for ground-based EHT extensions
are limited to a maximum physical baseline length of one Earth diameter.
Besides, frequencies higher than 345 GHz are accessible only for very few
sites with excellent atmospheric conditions, which makes VLBI image
reconstruction at these frequencies difficult, if not impossible.
Further enhancement of the EHT can be achieved with the deployment of EHT
antennas into space. The first observations which included space-ground
baselines were performed in 1986 by the Tracking and Data Relay Satellite
System (TDRSS) and ground-based telescopes in Australia and Japan (Levy et
al., 1986). One of the first space VLBI (SVLBI) missions to include space-
ground baselines was the VLBI Space Observatory Programme (VSOP) carried by
the Highly Advanced Laboratory for Communications and Astronomy (HALCA)
satellite that operated in orbit in the period from 1997 to 2003 (Hirabayashi
et al., 1998, 2000). Another SVLBI mission of this kind was RadioAstron, which
was operational between 2011 and 2019 (Kardashev et al., 2013).
In addition to the EHT array, VLBI stations could be located in a nearly
circular Low Earth Orbit (Palumbo et al., 2019) or a Highly Elliptical Orbit
(Andrianov et al., 2021). Such a setup would provide fast baseline coverage,
allowing for dynamical imaging of rapidly varying sources, as Sgr A*. VLBI
stations located in an elliptical medium Earth orbit (MEO) would increase the
angular resolution and the imaging capabilities of the array to distinguish
different space-time metrics around Sgr A* (Fromm et al., 2021). Inclusion of
one telescope in a high-inclination MEO or Geosynchronous Orbit would increase
the angular resolution and the number of available sources for which the
shadow could be resolved (Fish et al., 2019). For example, the central black
holes in IC 1459 and M84 (NGC 4374) with predicted shadow diameters of $\sim
9.2\leavevmode\nobreak\ \mu$as and $\sim 9.1\leavevmode\nobreak\ \mu$as,
respectively, as well as the Sombrero Galaxy (NGC 4594, M104) with an
estimated shadow diameter of $\sim 5.7\leavevmode\nobreak\ \mu$as and several
other sources (Johannsen et al., 2012). The Black Hole Explorer
(BHEX111https://www.blackholeexplorer.org/) mission concept aims to launch a
telescope into a near-geosynchronous orbit to make precision measurements of a
black hole’s photon ring properties using ground-space baselines (Kurczynski
et al., 2022, Johnson et al., Marrone et al., in prep., SPIE), which would
complement high-frequency space-space baseline imaging with the Event Horizon
Imager concept presented in this work. A station on the Moon or in the Earth-
Sun or Earth-Moon Lagrange points, added to the EHT array, would sharpen the
angular resolution sufficiently to resolve signatures of the substructures in
the photon rings of M87* and Sgr A* (Johnson et al., 2020). It was proposed
that the Origins Space Telescope, equipped with a modified version of its
heterodyne instrument, could be used for these purposes (Pesce et al., 2019).
An alternative approach is to consider only space-space baselines. Such a
system has the advantage that atmospheric data corruptions can be avoided,
thus allowing for observation at higher frequencies in addition to the
possibility of reaching baselines longer than an Earth diameter. Therefore,
the resolution can be increased even further.
The Event Horizon Imager (EHI) is a concept for a new SVLBI system that
consists of two or three satellites in polar or equatorial circular MEOs
observing with only space-space baselines at frequencies up to 690 GHz
(Martin-Neira et al., 2017; Kudriashov et al., 2019; Roelofs et al., 2019;
Kudriashov et al., 2021b). The concept envisages the precise localization of
the interferometer baseline based on the real-time relative positioning of the
satellites using the Global Navigation Satellite System (GNSS). In addition,
the concept suggests the on-board cross-correlation for the reduction of the
difficulties of the data downlink to the ground. The accurate positioning of
the satellites and interchange of local oscillator signals permit the phase
calibration of EHI observations and, thus, the use of complex visibilities
(Kudriashov et al., 2021b). Both the proposed on-the-fly positioning and the
on-board data processing require a working inter-satellite link (ISL) to
perform highly accurate ranging measurements, as well as to exchange observed
signals and local oscillator components for coherent operation (Martin-Neira
et al., 2019; Kudriashov et al., 2021a). This implies a direct line of sight
between the satellites. Hence the maximum baseline length, limiting the
angular resolution of the EHI, is set by the radius of the satellite orbits
and the occultation of the ISL by the Earth. The proposed resolution is $\sim
5$ $\mu$as (Kudriashov et al., 2021b), which is an order of magnitude better
than what can be obtained with the EHT from Earth. Due to higher observation
frequencies, the EHI can detect emission that originates from regions that are
closer to the black hole (Mościbrodzka et al., 2009). Those regions are more
dominated by general relativistic effects, leading to the reduced variability
of the images. Additionally, it causes the emission to trace the photon ring
more closely. Moreover, the proposed system can allow one to avoid a
significant part of the interstellar scattering for the observations of Sgr A*
since the scattering kernel size decreases with the square of the observing
frequency (Roelofs et al., 2019).
Simulated observations of Sgr A* with the EHI demonstrate the excellent
imaging capability with this concept (Roelofs et al., 2019). These simulations
have shown that the EHI could image Sgr A* with an order of magnitude higher
angular resolution and accuracy than the EHT within just a few months of
observation, assuming standard system noise parameters.
The overarching EHI science goal is to test theories of gravity. Precise
imaging of a black hole photon ring brings essential information to
distinguish between GR and alternative theories. The imaging of the plasma
variability would provide constraints on models of plasma dynamics and jet
formation. Therefore, the motivation for further testing of the EHI
capabilities is to study aspects limiting the fidelity of imaging on different
timescales of observations in order to determine how accurately GR, as well as
models of black hole accretion, plasma dynamics and jet launching, can be
tested with the EHI. In addition, understanding EHI imaging constraints can
help to optimize the design of the system.
This paper supplements previous EHI studies by examining several
configurations of the system. We study their influence on the possibility of
resolving structural variations of the M87* environment and extended jets of
other black holes (e.g. NGC 1052) at millimetre and sub-millimetre
wavelengths. The research is focused on M87* since its timescales of
structural variations are significantly longer than those of Sgr A*, as
mentioned above. The considered system setups are described in Section 2.
Section 3 introduces source models and image generation. The simulation
results are presented in Section 4, and conclusions are summarized in Section
5.

Figure 1: Time to reach the longest baseline as a function of orbital
separation for the considered EHI system setup. Blue line corresponds to the
most separate pair of satellites; red line corresponds to the pair of the
innermost and the middle satellites; green line corresponds to the pair of the
middle and the outermost satellites. Figure 2: Parts of the $uv$-coverage of
8.9-hour duration for 30, 400, and 1000 km separations of the orbits at 230
GHz; first, 81st, and 161st snapshots are shown (from top to bottom). Blue
points correspond to the most separate pair of satellites; red points
correspond to the pair of the innermost and the middle satellites; green
points correspond to the pair of the middle and the outermost satellites.

##  2 EHI system setup

The EHI concept implies two or three satellites in polar or equatorial
circular MEOs, as mentioned above. The particular configuration is still under
discussion. The EHI system considered in this paper consisted of three
satellites in circular MEOs with slightly different radii. The largest radius
of the satellite orbits was 13 913 km. The selected radius provides the
longest possible baselines taking into account the required simultaneous
accessibility of at least three GNSS satellites from the EHI orbits. At the
same time, this maximal radius ensures the orbits of all three EHI satellites
to be between the Van Allen belts (Kudriashov et al., 2019). The second
satellite was placed on an orbit with a radius smaller by a distance referred
to hereafter as the orbital separation. The third satellite was added between
the first two. Its orbit radius was one-third of the orbital separation bigger
than for the inner satellite. The orbital plane was set perpendicular to the
line of sight to the observed source, taking into account the declination. In
reality, this condition without additional reorientation of the orbital plane
can be fulfilled only for a very limited list of sources. Due to the different
orbital radii and velocities of the satellites, the complete $uv$-coverage
will have the shape of a spiral. We investigated systems with constant orbital
separations of 30, 50, 60, 100, 200, 300, 350, 400, 500, and 1000 km.
The imaging of the horizon-scale variability requires satisfactory coverage of
the whole $uv$-plane over time intervals comparable to the timescale of the
investigated variability. The longest available baseline is limited by the
radius of the satellite orbits and the occultation of the ISL by the Earth.
Let us assume that all three satellites start at the same orbital phase. When
the first pair of satellites reaches the longest baseline, the other two pairs
are also at the relatively long baselines. If the system continues observation
without configuration changes, the $uv$-coverage provided by the remaining
pairs of satellites lacks short baselines. This situation continues until the
occultation of the ISL between the first pair of satellites is finished and
they converge again. During this time, the other two pairs of satellites
sequentially find themselves in a similar situation, which creates a prolonged
period of increased $uv$-coverage sparseness. This sparsity period can in
principle be avoided by interchanging the satellite orbits each time when the
longest baseline length is reached. Figure 1 shows that the time to reach the
longest baseline depends on the orbital separation and rapidly decreases from
approximately a month to several days. Regular changes in the satellite orbits
on such short timescales are exceedingly fuel-consuming and, therefore,
unrealistic. The exclusion of satellite orbit changes will save fuel,
resulting in a prolongation of the overall mission duration and an increase in
both the amount of collected data and the number of objects that can be
observed. We considered satellites to stay in their orbits as long as required
to get a long enough series of $uv$-data points. The full data series can be
used for the reconstruction of a highly detailed time-averaged image. In this
paper, we divided the series of data into parts to obtain the time resolution
of interest. Figure 2 illustrates this concept with an example (see Sec. 3.2
for details). Therefore, the temporal resolution did not depend on orbital
separation and was selected during data processing.
In this work, we considered the standard EHI system composed of 4.0-metre
antennas, fitting in the Ariane 6 spacecraft (Arianespace, 2018), and a
hypothetical EHI+ system consisting of 15.0-metre antennas (see also Gurvits
et al., 2022). System parameters for the noise addition were taken from
Roelofs et al. (2019) for the EHI system. The EHI+ system was assumed to have
a system temperature close to the minimum allowable for a coherent heterodyne
system and a wider observing bandwidth to increase the signal-to-noise ratio
(Sec. 3.3). In addition, EHI+ antennas were considered to observe at higher
frequencies (see Sec. 3.1). Technical specifications assumed in this paper for
EHI+ are foreseen by the authors as technically very challenging and not yet
demonstrated, but not impossible. Phase corruptions due to uncertainties in
the orbital model have not been studied in this work; it was assumed that the
baseline vector is known exactly thanks to the GNSS-assisted orbit
determination, inter-satellite link ranging measurements, a local oscillator
sharing setup, and other measures that are currently under investigation
(Kudriashov et al., 2021b). Image reconstructions were therefore produced
using complex visibilities; the extent of the necessity of using closure
products will depend on the eventual feasibility and performance of the
connected-interferometer concept laid out in Kudriashov et al. (2021b) (see
also Sec. 3.2). An investigation of the influence of the detailed system noise
and calibration parameters on the reconstructed images presented here is left
for future work.

##  3 Simulated observations

In this section, the generation of simulated observations is described.
Theoretical emission maps were produced from GRMHD simulations and used as
input. These source models were used to calculate complex visibilities for the
simulated observations together with the coverage of the $uv$-plane, produced
by the system setup outlined in Section 2. The calculated complex visibilities
were supplemented by data corruption effects (e.g. thermal noise) and then
reconstructed with a Regularized maximum likelihood (RML) algorithm. The
quality of the image reconstructions was quantified by several image quality
metrics.

Table 1: Parameters of the theoretical emission maps used.

source | time interval (hours) | movie duration (months) |  $\nu$ (GHz) | total flux (Jy) | FoV ($\mu$as)
---|---|---|---|---|---
|  |  | 230 | 0.84 - 1.31 | 190.9
M87* | 89 (10 $t_{\mathrm{G}}$) | 10 (800 $t_{\mathrm{G}}$) | 560 | 0.47 - 0.84 | 190.9
|  |  | 5000 | 0.10 - 0.18 | 190.9
|  |  | 43 | 0.60 - 0.82 | 763.7
NGC 1052 | 89 | 10 | 43 | 7.0 - 7.6 | $9.6\times 10^{3}$
|  |  | 230 | 0.93 - 1.04 | 190.9
M87* | 8.9 (1 $t_{\mathrm{G}}$) | 2.5 (200 $t_{\mathrm{G}}$) | 560 | 0.48 - 0.56 | 190.9
|  |  | 5000 | 0.10 - 0.12 | 190.9
|  |  | 43 | 0.69 - 0.72 | 763.7
NGC 1052 | 8.9 | 2.5 | 43 | 7.39 - 7.46 | $9.6\times 10^{3}$

222Duration of movies was calculated as time intervals between frames
multiplied by the number of frames in sets. For $uv$-coverage calculation, the
same time intervals between frames in movies were used as for M87* since the
black hole mass was not properly scaled for NGC 1052 models when changing the
angular size of the source.

###  3.1 Theoretical emission maps

Modelling of accretion flows around black holes is typically done by
performing GRMHD simulations (see Porth et al., 2019, and references therein).
The majority of these simulations only solve for the dynamically important
protons, although more recent works also include information on the electron
population via either including electron thermodynamics in GRMHD (Ryan et al.,
2018; Chael et al., 2019b) or by using Particle-in-Cell simulations instead
(Parfrey et al., 2019; Crinquand et al., 2021; Bransgrove et al., 2021). The
GRMHD models used in this work do not contain information on the electron
population, and therefore use a parametrization for their properties, such as
the shape of the distribution function, in the subsequent ray tracing. In this
work, we used the $\kappa$-jet model of M87* from Davelaar et al. (2019),
first developed for Sgr A* in Davelaar et al. (2018). The $\kappa$-jet model
for M87* is capable of recovering the observed spectral energy distribution
from radio to near-infrared, the jet core shift relation and shows consistent
image morphologies at 230 GHz. The model assumes that the electron
distribution function is a $\kappa$-distribution function, which is a
combination of a thermal core and a power-law. The slope of the power-law is
set by a parametrization from Ball et al. (2018), who studied trans-
relativistic reconnection. The parameterization effectively adds a power-law
population in the jet sheath, while electrons in the disk are in a thermal
distribution.
The dynamics of the accretion flow onto the black hole were simulated using
the Black Hole Accretion Code (BHAC, Porth et al., 2017; Olivares et al.,
2019) that solves the GRMHD equations. The performed simulation assumed a
standard and normal evolution (SANE) accretion disk, the dimensionless black
hole spin parameter was set to be $a_{*}=15/16$. The general relativistic ray-
tracing code RAPTOR (Bronzwaer et al., 2018, 2020) was used to generate
synthetic synchrotron maps (for further details, see Davelaar et al., 2019).
The resulting collections of synthetic synchrotron maps of the jet-launching
region in M87* were used in this paper as input models to simulate
observations. These maps were calculated at an inclination of $i=160^{\circ}$,
which ensured that the orientation of the emitting region corresponds to the
results of the EHT observations (Event Horizon Telescope Collaboration et al.,
2019a, e). The collection of synthetic synchrotron maps represents a black
hole environment with a certain time interval between frames, defined by the
gravitational timescale of a black hole $t_{\mathrm{G}}$. In this work, the
models had 1 and 10 $t_{\mathrm{G}}$ (8.9 and 89 hours, respectively)
intervals between frames.
The analysis for the EHI system was carried out for three frequencies, namely
43 GHz, which is a standard frequency for many VLBI observations, 230 GHz,
which is the operating frequency of the EHT, and 560 GHz, which is an expected
operating frequency of the EHI. 560 GHz was selected, following one of the
secondary EHI science goals of imaging the water line in protoplanetary disks.
If future technical studies find that ground telescope support will be a
system requirement, this observing frequency may be adjusted to, for example,
690 GHz without large consequences for the observable black hole shadow
features.
In the case of the EHI+ system, the analysis was performed for two
frequencies, which are 560 GHz and 5 THz, to test the imaging limitations of
the EHI concept. The frequency of a few THz was selected since this
additionally increases the angular resolution and, therefore, the number of
sources for which a black hole shadow is resolvable, for instance, one in
Centaurus A (Janssen et al., 2021). It should be noted that 5 THz corresponds
to a wavelength of 60 microns, which translates to very stringent requirements
for the phase calibration and antenna surface accuracies of the EHI+ system.
To get models at 5 THz, synthetic synchrotron maps initially calculated at 560
GHz were scaled to the flux at 5 THz according to the flux density spectrum.
The flux density spectrum in the $\kappa$-jet model is described by a power-
law $F_{\nu}\propto\nu^{\alpha}$ with index $\alpha\approx-0.7$ at high
frequencies ($\nu>230$ GHz), as declared in Davelaar et al. (2019).
Another potential secondary EHI science goal is imaging the structural
variations of the extended jet of black holes at 43 GHz. Many active galactic
nuclei (AGNs) are located further away or have smaller black hole masses than
M87*, thus, shadows of these black holes can not be resolved even with the
EHI. Nevertheless, they are suitable for imaging relativistic jets on larger
scales. The considered M87* model provides information about the jet on scales
much shorter than in the observations. Hence additional synthetic synchrotron
maps with $i=90^{\circ}$ inclination of the source were scaled so that the
field of view and the total flux correspond to the parameters of the jet of
another AGN. We note that this scaling is not physical because the black hole
mass and hence the variability timescales are not properly scaled when
changing the angular size of the source. Nevertheless, these coarse maps allow
one to get insight into the ability to resolve AGN jet structures with the
EHI. As a generic example for an AGN jet, NGC 1052 was chosen (Kadler et al.,
2004; Ros & Kadler, 2008; Baczko et al., 2019; Nakahara et al., 2019; Gurvits
et al., 2021). Table 1 summarizes parameters of collections of synthetic
synchrotron maps that are used in this work.

###  3.2 Coverage of the $(u,v)$ plane

Calculation of the complex visibilities was performed with the `eht-imaging`
software (Chael et al., 2018, 2019a). This requires, besides the source model,
a $uv$-coverage as input. As described in Section 2, we fixed the satellite
orbits throughout the observations. The obtained $uv$-data was split in the
time domain into snapshots of a duration corresponding to the time interval
between frames in the theoretical emission maps, namely 1 and 10
$t_{\mathrm{G}}$ (8.9 and 89 hours). Therefore, each part of the $uv$-coverage
was related to the specific frame of the source model when calculating the
complex visibilities.
Following Roelofs et al. (2019), we set an integration time per measurement
that is $uv$-distance-dependent,

$t_{\mathrm{int}}=\frac{P}{4\pi D_{\lambda}\Theta},$ (1)

where $P$ is the orbital period of the innermost satellite, $D_{\lambda}$ is
the length of the corresponding baseline, and $\Theta$ is the field of view
(FoV) of the EHI system. Thereby $t_{\mathrm{int}}$ is within the
$uv$-smearing limit, so we can avoid the corruption of the reconstructed image
due to the displacement of the $uv$-vector during an integration time
(Thompson et al., 2017; Palumbo et al., 2019). For the shortest baselines, the
integration time was chosen so that the $uv$-arcs are limited to 10 degrees.
The FoV of the EHI system was chosen corresponding to the FoV, provided by the
model.
The time it takes to reach the longest baselines (Figure 1) and the spiral
density depend on the orbital separation. Hence the described $uv$-coverage
after splitting into snapshots demonstrates three distinctive possibilities of
points distribution on the $uv$-plane for each pair of satellites. Figure 2
illustrates these possibilities with an example of several snapshots for three
different orbital separations. The first option is a ‘spiral-like’
distribution. For this option, it is typical that the points are allocated
sparsely but homogeneously in the shape of a part of a spiral (e.g. Figure 2,
middle right panel). The second option is a ‘ring-shaped’ distribution. The
points are allocated tightly to each other, however, they are grouped into a
narrow ring (e.g. Figure 2, middle left panel). The third option is an
intermediate one between the previous two. For this option, it is typical that
a part of a spiral forms a ring of a medium width (Figure 2, bottom middle
panel). With a small orbital separation, a spiral-like distribution is present
only in the first snapshot (Figure 2, top left panel). In the other snapshots,
all three pairs of satellites demonstrate a ring-shaped distribution of points
on the $uv$-plane (Figure 2, middle left and bottom left panels). With an
increasing orbital separation, the $uv$-coverage from each pair of satellites
expands through an intermediate phase to a spiral-like distribution (Figure 2,
bottom panels).
The coverage of the $uv$-plane formed by the whole system is one, two or three
concentric narrow rings, depending on how many pairs of satellites are blocked
by the Earth, in most of the snapshots for the small orbital separations.
Therefore, the resulting $uv$-coverage provides information only on a very
limited range of baselines in each snapshot. When the orbital separation
increases, the range of baselines represented in each snapshot also grows due
to wider rings, providing more homogeneous coverage. As a trade-off, the
shortest baselines become longer and the coverage becomes less dense. Hence
optimizing the movie reconstruction quality requires finding a balance between
density and isotropy of the $uv$-coverage in each snapshot, taking the source
evolution timescales into account. Nevertheless, with a suitable orbital
separation, sufficient coverage of the $uv$-plane can be obtained without
additional changes in the satellite orbits during the observations.
Similar to the time it takes to reach the longest baselines, the time it takes
until the temporally blocked satellites catch up with the other satellites,
restoring the ISL as its path is no longer occulted by Earth, is inversely
proportional to the orbital separation. Therefore, over multiple iterations,
an increase of the orbital separation leads to an increase of the fractional
time period with the co-visibility of three satellites, which is necessary for
the calculation of closure phases (i.e. the phase of the bispectrum, the sum
of visibility phases on a closed triangle of baselines). For example, with an
orbital separation of 50 km, only $35\%$ frames of the 10 $t_{\mathrm{G}}$
duration are suitable for reconstructions with amplitudes and closure phases
instead of complex visibilities on the investigated observational period. 100
and 200 km separation provides $45\%$ appropriate frames, while 300 and 400 km
orbital separations grant $75\%$ and $90\%$ suitable frames, respectively.
Since closure phases are immune to station-based phase errors, and the
bispectrum imposes less strict requirements on the phase stability of the
interferometer, the use of these robust quantities relaxes technical system
requirements for the EHI. Therefore, the dependence of the possibility to
calculate closure phases on the orbital separations should be considered in
the design of the system.

Table 2: System parameters and resulting noise.

| EHI system | EHI+ system
---|---|---
| M87* | NGC 1052 | M87*
$\nu$ (GHz) | 230 | 560 | 43 | 43 | 560 | 5000
$D$ (m) | 4.0 | 15.0
$T_{\mathrm{sys}}$ (K) | 150 | 50 | 300
$\eta_{\mathrm{ap}}$ | 0.58 | 0.58
$\eta_{\mathrm{cor}}$ | 0.97 | 0.97
$\eta_{\mathrm{clock}}$ | 0.87 | 0.87
$\Delta\nu$ (GHz) | 5 | 25 | 450
$t_{\mathrm{int,centre}}$ (s) | 452.2 | 452.2
$t_{\mathrm{int,edge}}$ (s) | 65.67 | 27.12 | 87.82 | 6.99 | 27.12 | 3.02
$\mathrm{SEFD}$ (Jy) | $6.7\times 10^{4}$ | $1.6\times 10^{3}$
$\sigma_{\mathrm{centre}}$ (Jy) | 0.036 | 0.00038
$\sigma_{\mathrm{edge}}$ (Jy) | 0.094 | 0.15 | 0.082 | 0.29 | 0.0016 | 0.0047
333$\sigma$-values are calculated with equations 2 and 3 for 30 km orbital
separation at the centre (long integration time) and edge (short integration
time) of the $uv$-spiral. Integration time at the edge of the $uv$-spiral
$t_{\mathrm{int,edge}}$ at 43 GHz is different for NGC 1052 and M87* due to
different fields of view (FoVs).

###  3.3 System noise calculation

Thermal noise in the system can be described as a Gaussian noise in the
complex visibility plane with zero mean and standard deviation $\sigma$. In
radio interferometry, $\sigma$ can be calculated from the integration time
$t_{\mathrm{int}}$, the observing bandwidth $\Delta\nu$ and System Equivalent
Flux Densities of the antennas $\mathrm{SEFD}_{1,2}$ (Thompson et al., 2017).
For the two-bit sampling of the signal with four-level quantization,

$\sigma=\frac{1}{0.88}\sqrt{\frac{\mathrm{SEFD}_{1}\mathrm{SEFD}_{2}}{2\Delta\nu
t_{\mathrm{int}}}}.$ (2)

SEFD is standardly defined as

$\mathrm{SEFD}=\frac{2k_{\mathrm{B}}T_{\mathrm{sys}}}{\eta A},$ (3)

where $k_{\mathrm{B}}$ is Boltzmann’s constant, $T_{\mathrm{sys}}$ is the
system temperature, $\eta$ is the efficiency and $A=\pi(D/2)^{2}$ is the area
of an antenna with diameter D. The efficiency
$\eta=\eta_{\mathrm{ap}}\eta_{\mathrm{cor}}\eta_{\mathrm{clock}}$ includes
efficiencies of the aperture, correlator, and clock, respectively. All
antennas were assumed to have constant efficiency independently of
observational frequency. In practice, the considered efficiency can be
implemented for only one frequency, while at higher frequencies the efficiency
will be lower. The bandwidth of the ISL for all EHI configurations and EHI+,
observing at 560 GHz, was the sum of bandwidths in two polarizations, hence
equals $2\Delta\nu$. For EHI+, observing at 5 THz, the bandwidth was the sum
of bandwidths in two polarizations and two sidebands, so it equals
$4\Delta\nu$. The parameters and resulting noise for the considered system
setups are shown in Table 2. The noise increases towards the edge of the
spiral due to the $uv$-distance-dependent integration times (see Sec. 3.2).

###  3.4 Image reconstruction

With a finite $uv$-sampling, image reconstruction requires imposing
assumptions due to the non-uniqueness of images reproducing the observed data.
The CLEAN algorithm (Högbom, 1974) represents the sky image in terms of a
finite number of point sources to constrain the image. Regularized maximum
likelihood (RML) algorithms (Gull & Daniell, 1978) minimize the weighted sum
of the goodness-of-fit test statistic $\chi^{2}$ and a set of regularizers
that favour certain image characteristics such as smoothness, sparsity, or
similarity to a prior (e.g. Chael et al., 2016, 2018). In this work, two
methods of applying the RML algorithm were used: snapshot imaging and
dynamical imaging, following Johnson et al. (2017).
Imaging of the source dynamics implies long-lasting observations, which can be
used later for the reconstruction of a movie. This movie should display
changes in the surrounding environment of a black hole. In snapshot imaging, a
set of images, making up the movie, is reconstructed from a corresponding set
of observations, and each reconstruction is performed independently. We
selected a combination of the Gull-Skilling entropy function (Skilling & Gull,
1991) and the Total Squared Variation regularizer (Kuramochi et al., 2018)
from several regularization functions that are implemented into the RML
algorithm. The selection of parameters for snapshot imaging is explained in
Appendix A.1.
Dynamical imaging assumes that the images in the set are temporally connected,
each being a perturbation of the previous frame. This allows one to get
information that is lacking for high-quality reconstruction from previous and
subsequent frames. In dynamical imaging, we used a combination of the Total
Squared Variation regularizer, the Second Moment regularizer (Issaoun et al.,
2019) and a generic regularizer $\mathcal{R}_{\Delta t}$, enforcing continuity
from frame to frame in the reconstructed movie (Johnson et al., 2017). The
selection of parameters for dynamical imaging is explained in Appendix A.2.
We reconstructed movies with three temporal resolutions: 1, 10, and 100
$t_{\mathrm{G}}$ (8.9 hours, 3.7 days, and 37 days). The complex visibilities
were calculated from the source model frames and corresponding parts of the
$uv$-coverage frame by frame. Then data files with complex visibilities were
united according to the assumed temporal resolution of the reconstructed
movie. Movies with 10 $t_{\mathrm{G}}$ and 100 $t_{\mathrm{G}}$ temporal
resolutions were simulated based on source models with 1 $t_{\mathrm{G}}$ and
10 $t_{\mathrm{G}}$ time intervals between frames, respectively. Thus, the
variability of the source throughout the observation was included in the
reconstructed movies. For the comparison with the theoretical emission maps,
the latter was averaged over the corresponding number of frames.

###  3.5 Image quality metrics

Quantitative comparison of the reconstructed image with the true model image
was performed via two metrics, namely the Normalized Root-Mean-Square Error
(NRMSE) and the Normalized cross-Correlation (NXCORR). The NRMSE evaluates
images based on pixel-to-pixel similarities and is given by

$\mathrm{NRMSE}=\sqrt{\frac{\sum_{i=1}^{n^{2}}(I_{i}-I_{i}^{\prime})^{2}}{\sum_{i=1}^{n^{2}}I_{i}^{2}}},$
(4)

where $I_{i}$ is the intensity of the $i$th pixel of the $n\times n$ pixels
model image and $I_{i}^{\prime}$ is that of the reconstructed image (Chael et
al., 2016, 2018). For the best fit, the NRMSE should be minimized since it is
zero when the reconstructed image is identical to the true image.
The NXCORR (a sliding inner-product of two normalized functions) is determined
as the cross-correlation of the Fourier transforms of the normalized intensity
patterns of two images at different relative shifts (e.g. Event Horizon
Telescope Collaboration et al., 2019a; Issaoun et al., 2019). Since the NXCORR
compares the bulks of each intensity pattern, it is less sensitive to
individual features in the reconstructed image. For a given relative shift
$\boldsymbol{\delta}$,

$\mathrm{NXCORR}(\boldsymbol{\delta})=\mid\mathcal{F}^{-1}\\{\mathcal{F}\\{I_{\mathrm{norm}}^{*}(\boldsymbol{x})\\}\boldsymbol{\cdot}\mathcal{F}\\{I_{\mathrm{norm}}^{\prime}(\boldsymbol{x+\delta})\\}\\}\mid,$
(5)

where $\mathcal{F}$ is the Fourier transform operator, $I_{\mathrm{norm}}$ is
the normalized intensity pattern of the true image and
$I_{\mathrm{norm}}^{\prime}$ is the same for the reconstructed image. The
normalized intensity for each pixel $i$ in the image can be calculated as

$I_{\mathrm{norm},i}=\frac{I_{i}-\mu_{\mathrm{I}}}{\sigma_{\mathrm{I}}},$ (6)

where $\mu_{\mathrm{I}}$ and $\sigma_{\mathrm{I}}$ are the mean and standard
deviation of the intensity distribution in the image. One of all the different
shifts across the extent of the images resulting in the maximum cross-
correlation is then used to output the final NXCORR value for the two images.
Thus, the maximized NXCORR gives the best fit.
Despite the differences in comparison methods, the NRMSE and the NXCORR
demonstrate the same qualitative dependencies (see Appendix A). Therefore, for
comparison of movies obtained under different conditions, only NXCORR is
shown. To evaluate the quality of a set of images making up the movie, metrics
were calculated for each frame independently and then averaged over the
duration of the movie.

##  4 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 observations of the
M87* shadow with the EHI at 230 and 560 GHz. Further, we describe observations
with the EHI+ at 560 GHz and 5 THz. Finally, observations of AGN jets at 43
GHz are described.

###  4.1 Imaging of the M87* shadow

Simulations of M87* shadow observations were produced from models with 1 and
10 $t_{\mathrm{G}}$ time intervals between frames and then reconstructed with
temporal resolutions of 10 and 100 $t_{\mathrm{G}}$ respectively. These
simulations were reconstructed using both snapshot and dynamical imaging
methods. Additionally, observations with 1 $t_{\mathrm{G}}$ temporal
resolution were reconstructed using dynamical imaging, except the observations
at 5 THz due to the heaviness of the computations at this frequency. In this
subsection, we compare the quality of movie reconstructions and investigate
the optimal orbital separation.

Figure 3: Effect of the reconstruction method on the reconstruction quality.
From left to right: the time-averaged over 3.7 days theoretical emission map
of M87* at 230 GHz (the model with 8.9 hours between frames); frames of the
movies simulated for the EHI with 30 km orbital separation, each lasting 3.7
days, reconstructed independently (snapshot imaging method) and using
dynamical imaging (Johnson et al., 2017). The source is varying during the
simulated observation. Colours indicate brightness/pixel in mJy (square root
scale).

Figure 4: Effect of the temporal resolution on the reconstruction quality.
From left to right: the time-averaged over 3.7 days theoretical emission map
of M87* at 230 GHz (the model with 8.9 hours between frames); frames of the
movies simulated for the EHI with 30 km orbital separation, each lasting 3.7
and 37 days, reconstructed using snapshot imaging (to highlight the effect).
The source is varying during the simulated observation. Colours indicate
brightness/pixel in mJy (square root scale). Figure 5: Effect of the orbital
separation on the reconstruction quality. From left to right: the time-
averaged over 3.7 days theoretical emission map of M87* at 230 GHz (the model
with 8.9 hours between frames); frames of the movies simulated for the EHI
with 30, 400, and 1000 km orbital separations, each lasting 3.7 days,
reconstructed using snapshot imaging (to highlight the effect). The source is
varying during the simulated observation. Colours indicate brightness/pixel in
mJy (square root scale).

Figure 6: Quality of the M87* shadow movies obtained with different orbital
separations of the EHI system setup. The movie quality is shown with the
averaged normalized cross-correlation against the true image, or NXCORR, at
two frequencies: (1) 230 GHz, shown in the top panel; (2) 560 GHz, shown in
the bottom panel. Red, green and blue lines correspond to reconstructed movies
with temporal resolutions of 8.9 hours, 3.7 days and 37 days, respectively.
Dashed lines correspond to snapshot imaging; solid lines correspond to
dynamical imaging. The red semitransparent area indicates orbital separations
with the best quality of reconstructed movies, based on the quality of
individual frames in the movies. Figure 7: Reconstruction of a simulated M87*
shadow observation with the EHI at 230 GHz for 400 km orbital separation. From
left to right: the middle theoretical emission map in the period of 37 days
(the model with 3.7 days between frames); frames of the simulated movies, each
lasting 37 days, reconstructed using snapshot imaging and dynamical imaging
methods. For the EHI observing at 230 GHz, these movies demonstrate the best
quality, according to the NXCORR, among all reconstructions that image source
dynamics. The source is varying during the simulated observation. Colours
indicate brightness/pixel in mJy (square root scale). The full reconstruction
is available as an online movie. Figure 8: Same as Fig. 7, but for an
observation frequency of 560 GHz. The full reconstruction is available as an
online movie.

####  4.1.1 The EHI: 230 and 560 GHz

In Figure 3, we illustrate the quality of individual frames of movies
reconstructed with snapshot and dynamical imaging methods. Snapshot imaging
reconstructs each frame of the movie independently. Hence the absence of
baselines results in an inability to reconstruct some frames with satisfactory
quality, especially for the short orbital separation used. At the same time,
dynamical imaging gets information for image reconstruction from other frames
and, therefore, provides a visually noticeable enhancement of the image
quality throughout the entire movie. Some improvement can be achieved with
lower temporal resolution (Figure 4). The longer duration of movie frames
corresponds to bigger parts of the $uv$-spiral in each snapshot. Nevertheless,
apart from a stronger violation of the static source assumption, a decrease in
temporal resolution implies a loss of information about the source
variability. An increase of the orbital separation leads to a more uniform
$uv$-coverage and a simultaneous decrease in its density (discussed in Sec.
3.2). As demonstrated in Figure 5, the quality of frames in movies obtained
with large orbital separations is significantly higher than with a small
separation. In regions of high intensity, the accuracy and fidelity of movies
are similar for large orbital separations. However, some distortions and
artefacts corresponding to the $uv$-coverage sparsity appear at regions with
low intensity, when the separation is too wide. The visually observable
difference is confirmed by the image quality metrics. Since the discussed
issues of the reconstruction quality depend on the uniformity and density of
the $uv$-plane coverage, they are relevant for all considered frequencies.
Figure 6 shows the averaged quality of the movie reconstructions depending on
the orbital separation for observations with the EHI system at 230 GHz on the
top panel and 560 GHz on the bottom panel. As discussed earlier, small orbital
separations provide coverage of the $uv$-plane per frame in a limited range of
baselines. When the separation is wide enough to ensure a comparatively
uniform $uv$-coverage per each frame, the average quality of the movie is
plateauing. Outside of the best quality range of separations, the averaged
quality of movies stays on a plateau due to individual frames with high image
quality, while the majority of frames demonstrate slightly reduced quality,
compared to movies within the best quality range. The best results for all
simulations displayed in Figure 6 are obtained for orbital separations between
200 and 500 km, based on the quality of individual frames in the movies,
although the source variability can be reconstructed reasonably well for a
wider range of separations. The NXCORR metric indicates that simulations of
observations at 230 GHz with 100 $t_{\mathrm{G}}$ temporal resolution provide
a better quality of movies than ones with higher temporal resolutions
reconstructed with the same method. It also shows a slight advantage in favour
of the snapshot imaging reconstruction method compared to dynamical imaging.
Since the spatial frequency is proportional to the observational frequency
considering the same baselines, the $uv$-coverage is less dense at 560 GHz
compared to 230 GHz with the same system setup. Therefore, dynamical imaging
provides better quality of the reconstructed movies at this frequency for all
tested temporal resolutions with a small difference in quality between the
latter.
Figures 7 and 8 show several frames of the movies from the simulated
observation at 230 GHz and 560 GHz respectively with 400 km orbital separation
and 100 $t_{\mathrm{G}}$ temporal resolution. These movies demonstrate the
best quality among all reconstructions that image source dynamics at given
frequencies according to the NXCORR metric. Although some temporal model
brightness distribution changes are blurred out due to the temporal
resolution, the characteristic shape of the features matches the model across
the corresponding observational period. The right part of the images is
dominated by the emission from the jet, while structures on the left part of
the images display the movement of the plasma in the accretion disk.
Therefore, the demonstrated detailed imaging of the changes around the shadow
and at the beginning of the jet provides high-quality information about M87*
dynamics. The total flux of M87* at 560 GHz is lower than at 230 GHz (Table 1)
which makes the signal-to-noise ratio significantly lower and affects the
quality of the images. However, the spatial resolution is higher and the
environment closer to the event horizon can be probed. An improvement in the
quality of movies may be reached with a different selection of reconstruction
parameters, which were chosen based on simulations at 230 GHz in this work
(see Appendix A). Besides, observations at 560 GHz can be improved
significantly by increasing the system sensitivity (Sec. 4.1.2).

Figure 9: Same as Fig. 5, but for observations with the EHI+ at 560 GHz.

Figure 10: Same as Fig. 6, but for the EHI+ system setup observing at 560 GHz
(the top panel) and 5 THz (the bottom panel). Figure 11: Same as Fig. 8, but
for the EHI+. Simulated movies with frames lasting 8.9 hours are reconstructed
using the dynamical imaging method; ones with frames lasting 3.7 days are
reconstructed using snapshot imaging and dynamical imaging methods. The full
reconstruction is available as an online movie. Figure 12: Same as Fig. 7 and
8, but for the EHI+ with 200 km orbital separation, observing at 5 THz. The
full reconstruction is available as an online movie.

####  4.1.2 The EHI+: 560 GHz and 5 THz

The main trends in the image quality identified for the EHI system and
associated with the $uv$-coverage (Figures 3, 4 and 5) remain relevant for the
EHI+ system. Moreover, Figure 9 shows that reconstruction artefacts and
distortions observed at the largest orbital separations become clearer with an
increasing observation frequency due to a general decrease in $uv$-spiral
density. At the same time, the effect of the lower flux density at higher
frequencies is compensated by the lower noise in the EHI+ system.
Figure 10 shows the quality of the image reconstruction depending on the
orbital separation for observations with the EHI+ system at 560 GHz on the top
panel and 5 THz on the bottom panel. All simulations of observations with the
EHI+ at 560 GHz have a better quality of reconstructed movies compared to
observations with the EHI since EHI+ data have a higher signal-to-noise ratio.
Moreover, all displayed simulations demonstrate a similar quality regardless
of the reconstruction method and the temporal resolution. This shows that,
despite the issue of the $uv$-coverage per frame discussed above, the system
provides coverage that is dense and uniform enough for the very sharp and
accurate source reconstruction. The best results for all displayed simulations
of observations with the EHI+ at 560 GHz are obtained for orbital separations
between 300 and 500 km, based on the quality of individual frames in the
movies. Further increase of the observation frequency to 5 THz exacerbates the
$uv$-coverage issue due to the extra reduction of its density. Additional
imaging difficulties are introduced by the onward total flux decrease (see
Table 1). As demonstrated in Figure 10, data received in 10 $t_{\mathrm{G}}$
(3.7 days) becomes insufficient to ensure satisfactory quality of movies for
the source variability imaging. The temporal resolution of 100
$t_{\mathrm{G}}$ (37 days) improves image quality to some extent. In the case
of observations at 5 THz, the best quality range of separations lies between
60 and 400 km. The quality of movies obtained with orbital separations out of
this range is insufficient for the detailed imaging of the changes in the M87*
environment due to the unstable quality of frames in movies. Nevertheless, a 5
THz observation frequency can provide a resolution of $\sim 0.5$ $\mu$as,
which significantly enhances the imaging capabilities of the system.
Figure 11 shows several frames of the movies from the simulated observation at
560 GHz with 400 km orbital separation and temporal resolutions of 1
$t_{\mathrm{G}}$ and 10 $t_{\mathrm{G}}$. Although variability is hard to
notice even in the model on 1 $t_{\mathrm{G}}$ timescale, the reconstructed
movie reproduces features of the source with high accuracy. For 10
$t_{\mathrm{G}}$ temporal resolution, the snapshot imaging reconstruction
method gives slight artefacts, which mostly disappear with the dynamical
imaging method. The temporal resolution of 100 $t_{\mathrm{G}}$ provides a
lower quality of the movies since reconstruction artefacts intensify.
Observations with 200 km orbital separation and 100 $t_{\mathrm{G}}$ temporal
resolution reached the best quality among all at 5 THz. As seen in Figure 12,
both methods are able to reconstruct the main features of the source in
regions of high intensity. In low-intensity regions, movies demonstrate
artefacts, however, dynamical imaging significantly reduces them. Additional
improvement of the reconstruction quality may be reached with the dynamical
imaging regularizer $\mathcal{R}_{flow}$, which assumes that an image evolves
according to a steady flow of the flux density over time (Johnson et al.,
2017). In summary, all the reconstructed movies demonstrate that the EHI
concept allows for observations of the M87* dynamics from space at frequencies
up to 5 THz with the signal-to-noise ratio as the main limiting parameter.
Moreover, the angular resolution at 5 THz can provide exceptionally precise
measurements of black hole parameters and remarkably deep tests of different
theories of gravity, as well as the additional increase in the number of
sources with resolvable shadows.

Figure 13: Quality of black hole jet movies obtained with different orbital
separations of the EHI system setup at 43 GHz. The movie quality is shown with
the averaged NXCORR, for two sources: (1) M87*, shown in the top panel; (2)
NGC 1052, shown in the bottom panel. Green and magenta lines correspond to
reconstructed movies with a 3.7-day temporal resolution; blue and cyan lines
correspond to reconstructed movies with a 37-day temporal resolution. The red
semitransparent area indicates orbital separations with the best quality of
reconstructed movies, based on the quality of individual frames in the movies.

Figure 14: Reconstruction of a simulated M87* jet observation with the EHI at
43 GHz for 60 km orbital separation. From left to right: the middle
theoretical emission map in the period of 37 days (the model with 3.7 days
between frames); frames of the simulated movie, each lasting 37 days,
reconstructed using the snapshot imaging method. For the EHI observing at 43
GHz, this movie demonstrates the best quality, according to the NXCORR, among
all reconstructions that image source dynamics. The source is varying during
the simulated observation. Colours indicate brightness/pixel in mJy (square
root scale). The full reconstruction is available as an online movie.

Figure 15: Same as Fig. 14, but for NGC 1052, observed with 200 km orbital
separation. The full reconstruction is available as an online movie.

###  4.2 Imaging of an AGN jet

The possible secondary EHI science goal is imaging the variability in the
extended jets of various AGNs at 43 GHz. In this subsection, we research the
EHI imaging capability and the optimal orbital separation for this goal on the
example of M87* and NGC 1052. Simulations of imaging the M87* jet at 43 GHz
were produced similarly to the simulations of the M87* shadow imaging at 230
and 560 GHz. Observations of models with 1 and 10 $t_{\mathrm{G}}$ time
intervals between frames were reconstructed with temporal resolutions of 10
and 100 $t_{\mathrm{G}}$, respectively, considering the variability of the
source during the observation. For coarse models of NGC 1052 (see Sec. 3.1),
simulations were produced analogously. Observations of models with 8.9 and 89
hours between frames were reconstructed with temporal resolutions of 3.7 and
37 days per frame, respectively. All simulations were reconstructed using the
snapshot imaging method.
The image quality of reconstructions becomes more sensitive to the
$uv$-coverage issue (see Sec. 4.1) for the imaging of extended jets. Figure 13
shows the quality of reconstructed movies depending on the orbital separation
for M87* on the top panel and NGC 1052 on the bottom panel. As seen from M87*
and NGC 1052 jet imaging simulations, a relatively high temporal resolution
(3.7 days in this project) is insufficient for adequate observations. The
$uv$-coverage obtained during this period lacks the required amount of
baselines, and only a few frames in a movie can be reconstructed with
satisfactory quality. At the same time, movies with approximately a month-long
(37 days) temporal resolution demonstrate significantly better quality. In the
case of M87*, the best quality of reconstructed movies corresponds to orbital
separations in a range from 50 to 200 km, based on the quality of individual
frames in the movies. For NGC 1052, high-quality movies are produced with
orbital separations between 100 and 400 km. Orbital separations out of this
range provide abundant reconstruction artefacts.
Figure 14 shows several frames of the movie reconstructed from simulated
observations of M87* with 60 km orbital separation and 100 $t_{\mathrm{G}}$
temporal resolution. The spatial resolution at 43 GHz frequency is
significantly lower, but changes in the brightness distribution at the
beginning of the extended jet are visible and can be used to reconstruct jet
dynamics. The noticeable deficit of the observed intensity in the central
region of images can be confused with the shadow but it corresponds to the
spine of the jet directed away from the observer. There are limitations in the
modelling of extended jets since actual sources demonstrate much more extended
emission compared to our model. Thus, the signal-to-noise ratio is reduced in
our simulations, which decreases the quality of reconstructed movies. The best
result for simulated observation of the NGC 1052 jet is reached with 200 km
orbital separation. Figure 15 shows several frames of the corresponding movie
with clear changes in the jet shape and the brightness distribution along it
from frame to frame. The changes are most noticeable in the brightest parts of
the jet. For a more detailed study, physically correctly scaled models are
required in order to consider the timescales of the source variability.

Figure 16: Orbital separations and temporal resolutions providing the best
quality imaging of the source. The EHI system consists of three 4.0-metre
antennas with standard system parameters; the EHI+ system consists of three
15.0-metre antennas with more optimistic noise parameters. Simulated
observations at 230, 560 GHz, and 5 THz image the M87* shadow; simulated
observations at 43 GHz image M87* and NGC 1052 jets.

##  5 Conclusions and summary

In this paper, we have presented simulations of the EHI SVLBI system
consisting of three satellites in circular MEOs with slightly different radii
(Martin-Neira et al., 2017; Kudriashov et al., 2019; Roelofs et al., 2019;
Kudriashov et al., 2021b). The proposed accurate relative positioning of the
satellites and interchange of local oscillator signals allow for the use of
complex visibilities. The absence of atmospheric data corruptions allows for
imaging at frequencies significantly higher compared to ground-based
observations. In this work, we investigated the effects of different orbital
separations on the reconstructed image quality. The orbital separation can be
changed during the mission to the optimal separation for the currently
observed source and frequency.
The EHI setup provides a spiral-shaped sampling of the $uv$-plane and can
perform high-fidelity imaging of the M87* shadow without additional changes in
the satellite motion during observations. Nevertheless, dividing the
$uv$-coverage into snapshots, considered in this paper, requires a balance
between density and homogeneity in each snapshot that influences the selection
of the orbital separation. Too small separations provide coverage of the
$uv$-plane per frame in an overly narrow range of baselines and too large ones
produce excessively sparse $uv$-coverage, which results in distortions and
artefacts in the reconstructed movies. An increase in the observation
frequency also enlarges the $uv$-coverage sparsity.
Using GRMHD simulations of M87* and model system parameters, similar to those
discussed in Roelofs et al. (2019), we performed simulated long-lasting
observations to assess the quality of movies that can be expected. For the
imaging of the structural variations of the M87* environment, an important
potential EHI science goal, simulations were performed with temporal
resolutions of 1, 10, and 100 $t_{\mathrm{G}}$ (8.9 hours, 3.7 days, and 37
days) at 230 and 560 GHz for the standard EHI system, 560 GHz and 5 THz for
the EHI+ system with more optimistic noise parameters. Also, simulations for
the EHI were performed with temporal resolutions of 3.7 and 37 days at 43 GHz
for another potential goal, imaging the structural variations of extended
jets. The simulated reconstructed movies with temporal resolutions of 3.7 and
37 days included source variability during each frame of the resulting movies.
Figure 16 summarizes results of our simulations. The M87* shadow high-fidelity
imaging is possible at 230 and 560 GHz with the orbital separation in the
range from 200 to 500 km for the EHI. The best quality of the reconstructions
has been obtained by observations with 400 km orbital separation and 100
$t_{\mathrm{G}}$ (37 days) temporal resolution at both frequencies. Since the
resolvable variability of the source has approximately a monthly timescale, a
37-day temporal resolution is optimal for observations with the standard EHI
system. Simulations of the M87* shadow imaging at 560 GHz with the EHI+ system
demonstrate the highest accuracy and fidelity with the orbital separation in
the range from about 300 to 500 km and the best performance for 400 km
separation. The signal-to-noise ratio of observations with the EHI+ system at
560 GHz allows for the M87* environment variability imaging with 1
$t_{\mathrm{G}}$ (8.9 hours) temporal resolution. Changes in the M87*
environment are hardly resolvable on the gravitational timescale in the model
itself, nevertheless, a decrease in the temporal resolution is not favourable
since it produces artefacts in the reconstructed movies. In the case of
observations with the EHI+ at 5 THz, simulations only with orbital separations
in the range from 60 to 400 km and 100 $t_{\mathrm{G}}$ (37 days) temporal
resolution demonstrate the sufficient quality of movies for the detailed
imaging of the changes in the M87* environment. The best result at 5 THz,
obtained with 200 km orbital separation, shows accurate imaging of the source
main features in high-intensity regions, however, the $uv$-coverage sparsity
at such a high frequency leads to strong reconstruction artefacts in low-
intensity regions. Nevertheless, the angular resolution at 5 THz allows for
extraordinarily precise measurements.
The imaged changes in the brightness distribution around the shadow and at the
beginning of the jet are expected to provide information sufficient to
reconstruct the dynamics of the source. Such deep probes of the surrounding
environment of M87* can allow for deeper tests for general relativity and
alternative theories of gravity since these theories make predictions of the
appearance of the radio emission generated by material falling into the black
hole. Moreover, more accurate models of a black hole environment can be
tested. The emission at a higher frequency originates closer to the event
horizon, however, the M87* total flux is gradually decreasing at frequencies
higher than 230 GHz (Davelaar et al., 2019). Simulated observations
demonstrate that the reduction of the system noise level noticeably improves
the quality of reconstructed movies and makes source changes accessible for
imaging on shorter timescales and at higher frequencies. Therefore, the
signal-to-noise ratio is the main parameter limiting the EHI in observations
of the M87* shadow, which needs to be considered to optimize the efficiency
and the design of the system in general.
The resolution at 43 GHz is not high enough to resolve the shadow of M87* but
sufficient for imaging changes of the brightness distribution at the beginning
of the M87* jet. Orbital separations in the range from 50 to 200 km and 100
$t_{\mathrm{G}}$ (37 days) temporal resolution are expected to be the most
favourable for observations with the best result for 60 km separation. The EHI
capability for imaging jets of other sources was performed using NGC 1052 as a
generic example. The M87* models were scaled to match the field of view and
the total flux. This scaling is not physical because the black hole mass and
hence the variability timescales were not properly scaled when changing the
angular size of the source, therefore, additional research is essential.
Nevertheless, we have demonstrated that jet observations are exacting to
orbital separation and require lower temporal resolution to obtain sufficient
$uv$-coverage. Observations with a separation in the range from 100 to 400 km
and a 37-day temporal resolution can capture changes in the jet shape and the
brightness distribution along it. The best result has been achieved with 200
km orbital separation. Therefore, EHI observations at 43 GHz can provide
sufficient information to reconstruct structural variations of the
relativistic jets of M87* and other AGNs, for example, NGC 1052, 3C 273, 3C
84, M81* or Centaurus A. These observations can improve our understanding of
the jet launching and collimation processes.
Apart from M87*, the EHI is capable of imaging Sgr A*, as discussed by Roelofs
et al. (2019). Moreover, shadows of some other sources can potentially be
resolved by the EHI, such as the shadows of supermassive black holes at the
centres of the Sobrero Galaxy M104 and the elliptical galaxy M84 (Johannsen et
al., 2012). The EHI+ system observing at 5 THz could potentially resolve
shadows of M81* (Johannsen et al., 2012) and Centaurus A (Janssen et al.,
2021). With all these capabilities, the EHI concept will be of great
astrophysical interest.
The implementation of the concept into an actual mission implies several
technical challenges. The maximum possible orbit reconstruction accuracy
depends on laser ranging, accelerometers, orbital modelling, fringe fitting,
and other measures that are currently under investigation. Sending reduced
data to the ground requires the elaboration of the on-board correlation and
processing. Consequently, realistic observations may hinder image quality due
to uncertainties not considered in this paper. Nevertheless, we demonstrate
the general processes in the temporal structure of the EHI coverage. Moreover,
the signal-to-noise ratio is the primary parameter limiting EHI imaging
capabilities as discussed above, therefore, the possibilities of reducing the
system noise should also be investigated.
Technical system requirements for the EHI can be relaxed with the application
of closure phases to accompany visibility amplitudes instead of complex
visibilities. However, closure phase calculation requires data from all three
pairs of satellites for each frame of the reconstructed movies. Depending on
the orbital separation, partial sampling of the $uv$-plane, considered in this
paper, includes information from three pairs of satellites on different
fractions of the frames on the investigated observational period. $35\%$ of
appropriate frames for the 50 km orbital separation increase to $45\%$ for
separations of 100 and 200 km and $75\%$ and $90\%$ for 300 and 400 km ones,
respectively. Therefore, the dependence of the possibility to calculate
closure phases on the orbital separations should be considered in the design
of the system.
Additional consideration should also be given to the effect of the angle
between the orbital plane and the line of sight on the imaging quality.
Moreover, instead of a three-satellite system, two satellites can be
considered, as it is also one of the possible configurations. Another
possibility is investigating a space-space-ground hybrid system that would
provide fast baseline coverage for dynamical imaging of rapidly varying
sources such as Sgr A* (Palumbo et al., 2019).

###### Acknowledgements.

This work is supported by the ERC Synergy Grant “BlackHoleCam: Imaging the
Event Horizon of Black Holes” (Grant 610058). The authors thank Manuel Martin-
Neira, Volodymyr Kudriashov and Daniel Palumbo for their helpful comments and
discussions on this work. We are grateful to the anonymous referee for useful
and constructive comments. AS personally acknowledges Olesya Kuchay for her
unlimited support during the period of writing the manuscript and Jeremy
Tregloan-Reed for his invaluable advice. FR was supported by NSF grants
AST-1935980 and AST-2034306. This work was supported by the Black Hole
Initiative, which is funded by grants from the John Templeton Foundation
(Grant #62286) and the Gordon and Betty Moore Foundation (Grant GBMF-8273) -
although the opinions expressed in this work are those of the author(s) and do
not necessarily reflect the views of these Foundations. JD is supported by
NASA grant NNX17AL82G and a Joint Columbia/Flatiron Postdoctoral Fellowship.
Research at the Flatiron Institute is supported by the Simons Foundation. The
GRMHD simulations were performed on the Dutch National Supercomputing cluster
Cartesius and are funded by the NWO computing grant 16431.

## References

  * Abuter et al. (2019) Abuter, R., Amorim, A., Bauböck, M., et al. 2019, A&A, 625, L10
  * Andrianov et al. (2021) Andrianov, A. S., Baryshev, A. M., Falcke, H., et al. 2021, MNRAS, 500, 4866
  * Arianespace (2018) Arianespace. 2018, Ariane 6 User’s Manual
  * Bacchini et al. (2021) Bacchini, F., Mayerson, D. R., Ripperda, B., et al. 2021, Phys. Rev. Lett., 127
  * Baczko et al. (2019) Baczko, A.-K., Schulz, R., Kadler, M., et al. 2019, A&A, 623, A27
  * Ball et al. (2018) Ball, D., Sironi, L., & Özel, F. 2018, ApJ, 862, 80
  * Bransgrove et al. (2021) Bransgrove, A., Ripperda, B., & Philippov, A. 2021, Phys. Rev. Lett., 127
  * Broderick & Loeb (2009) Broderick, A. E. & Loeb, A. 2009, ApJ, 697, 1164–1179
  * Bronzwaer et al. (2018) Bronzwaer, T., Davelaar, J., Younsi, Z., et al. 2018, A&A, 613, A2
  * Bronzwaer & Falcke (2021) Bronzwaer, T. & Falcke, H. 2021, The Nature of Black Hole Shadows
  * Bronzwaer et al. (2020) Bronzwaer, T., Younsi, Z., Davelaar, J., & Falcke, H. 2020, A&A, 641, A126
  * Chael et al. (2019a) Chael, A., Bouman, K., Johnson, M., et al. 2019a, eht-imaging: v1.1.0: Imaging interferometric data with regularized maximum likelihood
  * Chael et al. (2019b) Chael, A., Narayan, R., & Johnson, M. D. 2019b, MNRAS, 486, 2873–2895
  * Chael et al. (2018) Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23
  * Chael et al. (2016) Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ, 829, 11
  * Chang et al. (2010) Chang, C. S., Ros, E., Kovalev, Y. Y., & Lister, M. L. 2010, A&A, 515, A38
  * Chen et al. (2020) Chen, Y., Shu, J., Xue, X., Yuan, Q., & Zhao, Y. 2020, Phys. Rev. Lett., 124, 061102
  * Crinquand et al. (2021) Crinquand, B., Cerutti, B., Dubus, G., Parfrey, K., & Philippov, A. 2021, A&A, 650, A163
  * Davelaar et al. (2018) Davelaar, J., Mo´scibrodzka, M., Bronzwaer, T., & Falcke, H. 2018, A&A, 612, A34
  * Davelaar et al. (2019) Davelaar, J., Olivares, H., Porth, O., et al. 2019, A&A, 632, A2
  * Dexter et al. (2012) Dexter, J., McKinney, J. C., & Agol, E. 2012, MNRAS, 421, 1517–1528
  * Doeleman et al. (2019) Doeleman, S., Blackburn, L., Dexter, J., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 256
  * Doeleman et al. (2023) Doeleman, S. S., Barrett, J., Blackburn, L., et al. 2023, Galaxies, 11, 107
  * Doeleman et al. (2008) Doeleman, S. S., Weintroub, J., Rogers, A. E. E., et al. 2008, Nature, 455, 78–80
  * Event Horizon Telescope Collaboration et al. (2022a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022a, ApJ, 930, L12
  * Event Horizon Telescope Collaboration et al. (2022b) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022b, ApJ, 930, L13
  * Event Horizon Telescope Collaboration et al. (2022c) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022c, ApJ, 930, L14
  * Event Horizon Telescope Collaboration et al. (2022d) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022d, ApJ, 930, L15
  * Event Horizon Telescope Collaboration et al. (2022e) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022e, ApJ, 930, L16
  * Event Horizon Telescope Collaboration et al. (2022f) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022f, ApJ, 930, L17
  * Event Horizon Telescope Collaboration et al. (2023) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2023, ApJ, 957, L20
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a, ApJ, 875, L1
  * Event Horizon Telescope Collaboration et al. (2019b) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019b, ApJ, 875, L2
  * Event Horizon Telescope Collaboration et al. (2019c) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019c, ApJ, 875, L3
  * Event Horizon Telescope Collaboration et al. (2019d) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019d, ApJ, 875, L4
  * Event Horizon Telescope Collaboration et al. (2019e) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019e, ApJ, 875, L5
  * Event Horizon Telescope Collaboration et al. (2019f) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019f, ApJ, 875, L6
  * Event Horizon Telescope Collaboration et al. (2021a) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021a, ApJ, 910, L12
  * Event Horizon Telescope Collaboration et al. (2021b) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021b, ApJ, 910, L13
  * Falcke et al. (2000) Falcke, H., Melia, F., & Agol, E. 2000, ApJ, 528, L13–L16
  * Fish et al. (2019) Fish, V. L., Shea, M., & Akiyama, K. 2019, Adv. Space Res.
  * Fromm et al. (2021) Fromm, C. M., Mizuno, Y., Younsi, Z., et al. 2021, A&A, 649, A116
  * Gebhardt et al. (2011) Gebhardt, K., Adams, J., Richstone, D., et al. 2011, ApJ, 729, 119
  * Gull & Daniell (1978) Gull, S. F. & Daniell, G. J. 1978, Nature, 272, 686
  * Gurvits et al. (2022) Gurvits, L. I., Paragi, Z., Amils, R. I., et al. 2022, Acta Astronautica, 196, 314
  * Gurvits et al. (2021) Gurvits, L. I., Paragi, Z., Casasola, V., et al. 2021, Experimental Astronomy
  * Hada et al. (2016) Hada, K., Kino, M., Doi, A., et al. 2016, ApJ, 817, 131
  * Hirabayashi et al. (2000) Hirabayashi, H., Hirosawa, H., Kobayashi, H., et al. 2000, PASJ, 52, 955
  * Hirabayashi et al. (1998) Hirabayashi, H., Hirosawa, H., Kobayashi, H., et al. 1998, Science, 281, 1825
  * Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
  * Issaoun et al. (2019) Issaoun, S., Johnson, M. D., Blackburn, L., et al. 2019, A&A, 629, A32
  * Janssen et al. (2021) Janssen, M., Falcke, H., Kadler, M., et al. 2021, Nat. Astron.
  * Johannsen et al. (2012) Johannsen, T., Psaltis, D., Gillessen, S., et al. 2012, ApJ, 758, 30
  * Johnson et al. (2023) Johnson, M. D., Akiyama, K., Blackburn, L., et al. 2023, Galaxies, 11, 61
  * Johnson et al. (2017) Johnson, M. D., Bouman, K. L., Blackburn, L., et al. 2017, ApJ, 850, 172
  * Johnson et al. (2015) Johnson, M. D., Fish, V. L., Doeleman, S. S., et al. 2015, Science, 350, 1242–1245
  * Johnson et al. (2020) Johnson, M. D., Lupsasca, A., Strominger, A., et al. 2020, Science Advances, 6, eaaz1310
  * Kadler et al. (2004) Kadler, M., Kerp, J., Ros, E., et al. 2004, A&A, 420, 467–474
  * Kardashev et al. (2013) Kardashev, N. S., Khartov, V. V., Abramov, V. V., et al. 2013, Astron. Rep., 57, 153–194
  * Kravchenko et al. (2020) Kravchenko, E., Giroletti, M., Hada, K., et al. 2020, A&A, 637, L6
  * Kudriashov et al. (2019) Kudriashov, V., Martin-Neira, M., Barat, I., et al. 2019, CJSS, 39, 250
  * Kudriashov et al. (2021a) Kudriashov, V., Martin-Neira, M., Lia, E., et al. 2021a, Journal of Astronomical Instrumentation, 10, 2150010
  * Kudriashov et al. (2021b) Kudriashov, V., Martin-Neira, M., Roelofs, F., et al. 2021b, CJSS, 41, 211
  * Kuramochi et al. (2018) Kuramochi, K., Akiyama, K., Ikeda, S., et al. 2018, ApJ, 858, 56
  * Kurczynski et al. (2022) Kurczynski, P., Johnson, M. D., Doeleman, S. S., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12180, Space Telescopes and Instrumentation 2022: Optical, Infrared, and Millimeter Wave, ed. L. E. Coyle, S. Matsuura, & M. D. Perrin, 121800M
  * Levy et al. (1986) Levy, G. S., Linfield, R. P., Ulvestad, J. S., et al. 1986, Science, 234, 187
  * Lu et al. (2018) Lu, R.-S., Krichbaum, T. P., Roy, A. L., et al. 2018, ApJ, 859, 60
  * Marshall et al. (2002) Marshall, H. L., Miller, B. P., Davis, D. S., et al. 2002, ApJ, 564, 683–687
  * Martin-Neira et al. (2017) Martin-Neira, M., Kudriashov, V., Barat, I., Duesmann, B., & Daganzo, E. 2017, in Proceedings of the 5th Workshop on Advanced RF Sensors and Remote Sensing Instruments (ARSI’17), ESTEC (the Netherlands), 12-14 September 2017
  * Martin-Neira et al. (2019) Martin-Neira, M., Kudriashov, V., Barat, I., Duesmann, B., & Daganzo, E. 2019, CJSS, 39, 544
  * Mizuno et al. (2018) Mizuno, Y., Younsi, Z., Fromm, C. M., et al. 2018, Nature Astronomy, 2, 585–590
  * Mościbrodzka et al. (2017) Mościbrodzka, M., Dexter, J., Davelaar, J., & Falcke, H. 2017, Monthly Notices of the Royal Astronomical Society, 468, 2214–2221
  * Mościbrodzka et al. (2016) Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A, 586, A38
  * Mościbrodzka et al. (2009) Mościbrodzka, M., Gammie, C. F., Dolence, J. C., Shiokawa, H., & Leung, P. K. 2009, ApJ, 706, 497
  * Nakahara et al. (2019) Nakahara, S., Doi, A., Murata, Y., et al. 2019, AJ, 159, 14
  * Narayan et al. (2003) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69–L72
  * Olivares et al. (2019) Olivares, H., Porth, O., Davelaar, J., et al. 2019, A&A, 629, A61
  * Olivares et al. (2020) Olivares, H., Younsi, Z., Fromm, C. M., et al. 2020, Monthly Notices of the Royal Astronomical Society, 497, 521–535
  * Palumbo et al. (2019) Palumbo, D. C. M., Doeleman, S. S., Johnson, M. D., Bouman, K. L., & Chael, A. A. 2019, ApJ, 881, 62
  * Parfrey et al. (2019) Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122
  * Pesce et al. (2019) Pesce, D. W., Haworth, K., Melnick, G. J., et al. 2019, Extremely long baseline interferometry with Origins Space Telescope
  * Porth et al. (2019) Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26
  * Porth et al. (2017) Porth, O., Olivares, H., Mizuno, Y., et al. 2017, Computational Astrophysics and Cosmology, 4
  * Psaltis et al. (2020) Psaltis, D., Medeiros, L., Christian, P., et al. 2020, Phys. Rev. Lett., 125
  * Psaltis et al. (2021) Psaltis, D., Talbot, C., Payne, E., & Mandel, I. 2021, Phys. Rev. D, 103
  * Raymond et al. (2021) Raymond, A. W., Palumbo, D., Paine, S. N., et al. 2021, ApJS, 253, 5
  * Roelofs et al. (2019) Roelofs, F., Falcke, H., Brinkerink, C., et al. 2019, A&A, 625, A124
  * Roelofs et al. (2021) Roelofs, F., Fromm, C. M., Mizuno, Y., et al. 2021, A&A, 650, A56
  * Ros & Kadler (2008) Ros, E. & Kadler, M. 2008, J. Phys. Conf. Ser., 131, 012056
  * Rudin et al. (1992) Rudin, L. I., Osher, S., & Fatemi, E. 1992, Phys. D: Nonlinear Phenom., 60, 259
  * Ryan et al. (2018) Ryan, B. R., Ressler, S. M., Dolence, J. C., Gammie, C., & Quataert, E. 2018, ApJ, 864, 126
  * Sargent et al. (1978) Sargent, W. L. W., Young, P. J., Boksenberg, A., et al. 1978, ApJ, 221, 731
  * Skilling & Gull (1991) Skilling, J. & Gull, S. F. 1991, Lect. Notes-Monogr. Ser., 20, 341
  * Sparks et al. (1996) Sparks, W. B., Biretta, J. A., & Macchetto, F. 1996, ApJ, 473, 254
  * Thompson et al. (2017) Thompson, A. R., Moran, J. M., & Swenson, George W., J. 2017, Interferometry and Synthesis in Radio Astronomy, 3rd Edition
  * Vagnozzi et al. (2022) Vagnozzi, S., Roy, R., Tsai, Y.-D., et al. 2022, Horizon-scale tests of gravity theories and fundamental physics from the Event Horizon Telescope image of Sagittarius A∗
  * Van der Gucht et al. (2020) Van der Gucht, J., Davelaar, J., Hendriks, L., et al. 2020, A&A, 636, A94
  * Young et al. (1978) Young, P. J., Westphal, J. A., Kristian, J., Wilson, C. P., & Landauer, F. P. 1978, ApJ, 221, 721

##  Appendix A Parameters of the image reconstruction

###  A.1 Snapshot imaging

An RML algorithm minimizes the weighted sum of $\chi^{2}$ and a regularizer
(Gull & Daniell 1978). $\chi^{2}$ is the goodness-of-fit test statistic that
compares the visibilities of the test image to the data. The regularizer
contains prior information of the image, such as image smoothness or sparsity.
The RML algorithm, by its definition, gives an image that has no more
structures than required to fit the data. For snapshot imaging, we chose to
use two regularization terms.
One of the possible regularizers implemented into the `eht-imaging` software
is the Gull-Skilling entropy function (Skilling & Gull 1991). The Gull-
Skilling entropy function is a general form of the entropy function derived
from a Bayesian interpretation of the RML. In combination with the Total
Squared Variation regularizer, it shows the best visible result of the image
reconstruction during the regularizer selection and, therefore, was picked in
this work. The Total Squared Variation regularization (Kuramochi et al. 2018)
is the denoising process. The Total Squared Variation regularization is the
modification of the Total Variation regularization (Rudin et al. 1992). This
modification favours images with smaller variations between adjacent pixels,
which leads to the smoothing of edgelike features. Apart from better
consistency with observational data compared to post-processing smoothing,
edge-smoothed images demonstrate a better representation of the model images
in the case of diffuse astronomical objects.
The selection of the weighting parameters for the data term and both
regularization terms was performed by simulating an observation of the time-
averaged 10 $t_{\mathrm{G}}$ M87* model. The used system setup had a 50 km
orbital separation and an observing frequency of 230 GHz. Firstly, the
weighting parameters were varied from 10 to $10^{4}$. Figure 17 shows the
best-performing combinations of weights according to our image quality
metrics. The NRMSE performs a pixel-to-pixel comparison, while the NXCORR
compares the bulks of intensity patterns. Nevertheless, both metrics show
similar dependencies. The best performance according to image quality metrics
is achieved when the data term weighting $\alpha_{data}$ is bigger than the
weighting for Gull-Skilling entropy $\alpha_{GS}$. However, reconstructed
images are noticeably better when these two terms are equally weighted (Figure
18). The weighting parameter for the Total Squared Variation regularization
$\alpha_{TSV}$ demonstrates little effect on the quality of reconstruction.
When it is significantly less than the other two terms, the metric performance
is better. However, when this condition is reached, further alteration of the
Total Squared Variation regularization weight provides only a vanishingly
small variation in image quality metrics (Figure 19). Thus, the weighting
parameters for the data term $\alpha_{data}$ and the Gull-Skilling entropy
function term $\alpha_{GS}$ were both selected to equal 1000, the weighting
parameter for the Total Squared Variation regularization term $\alpha_{TSV}$
was selected to equal 200. To facilitate comparisons across different
frequencies and orbital separations, the weighting parameters for the data and
regularization terms are identical for all snapshot imaging reconstructions
displayed in this paper.
For both snapshot and dynamical imaging, the image reconstruction was
performed in several (typically 3-7) rounds to help the optimizer converge to
a single image structure. In the first round, a circular Gaussian with a FWHM
of 55 $\mu$as (700 $\mu$as for NGC 1052 jet simulations) and total flux equal
to that of the corresponding model frame were used as the initial and prior
image. In the following imaging rounds, the result of the previous round was
blurred to the nominal array resolution and used as a new initial and prior
image.

Figure 17: Quality of the reconstructions performed using the snapshot imaging
method with different weighting parameters. The image quality is measured in
two ways: (1) the normalized root-mean-square error against the true image, or
NRMSE, shown in the top panel; (2) the normalized cross-correlation against
the true image, or NXCORR, shown in the bottom panel. $\alpha_{data}$ is the
data term weighting parameter, $\alpha_{GS}$ is the Gull-Skilling entropy
function weighting parameter, $\alpha_{TSV}$ is the Total Squared Variation
regularization weighting parameter. Figure 18: Reconstructions, performed
using the snapshot imaging method with different weighting parameters. Time-
averaged M87* observation simulated at 230 GHz, 50 km orbital separation.
$a\\_data$ is the data term weighting parameter $\alpha_{data}$, $a\\_GS$ is
the Gull-Skilling entropy function weighting parameter $\alpha_{GS}$,
$a\\_TSV$ is the Total Squared Variation regularization weighting parameter
$\alpha_{TSV}$.

Figure 19: Same as Fig. 17, but with enlarged scale.

Figure 20: Quality of the M87* movies reconstructed using the dynamical
imaging method with different values of the dynamical regularizer parameter
$\sigma_{\Delta t}$. The image quality is measured in two ways: (1) the NRMSE,
shown in the top panel; (2) the NXCORR, shown in the bottom panel. Green and
blue lines correspond to reconstructed movies with temporal resolutions of 3.7
days and 37 days, respectively. $\alpha_{data}$ is the data term weighting
parameter, $\alpha_{TSV}$ is the Total Squared Variation regularization
weighting parameter, $\alpha_{SM}$ is the Second Moment regularization
weighting parameter, $\alpha_{\Delta t}$ is the dynamical regularizer
$\mathcal{R}_{\Delta t}$ weighting parameter.

Figure 21: Quality of the M87* movies reconstructed using the dynamical
imaging method with different values of the weighting parameter of the
dynamical regularizer $\mathcal{R}_{\Delta t}$. The image quality is measured
in two ways: (1) the NRMSE, shown in the top panel; (2) the NXCORR, shown in
the bottom panel. $\alpha_{data}$ is the data term weighting parameter,
$\alpha_{SM}$ is the Second Moment regularization weighting parameter,
$\alpha_{TSV}$ is the Total Squared Variation regularization weighting
parameter.

Figure 22: Quality of the M87* movies reconstructed using the dynamical
imaging method with different weighting parameters of static regularizers. The
image quality is measured in two ways: (1) the NRMSE, shown in the top panel;
(2) the NXCORR, shown in the bottom panel. $\alpha_{data}$ is the data term
weighting parameter, $\alpha_{\Delta t}$ is the dynamical regularizer
$\mathcal{R}_{\Delta t}$ weighting parameter, $\alpha_{SM}$ is the Second
Moment regularization weighting parameter, $\alpha_{TSV}$ is the Total Squared
Variation regularization weighting parameter.

Figure 23: Same as Fig. 22, but with enlarged scale.

###  A.2 Dynamical imaging

Dynamical imaging assumes that the images in the set are connected and gets
information for image reconstruction from other frames (Johnson et al. 2017).
This allows the missing data to be partially inherited from the rest of the
movie by interpolation. In this work, a combination of the dynamical
regularization with two static regularization terms was used.
A generic regularizer $\mathcal{R}_{\Delta t}$ from Johnson et al. (2017),
used in this paper, assumes that the images in the set are temporally
connected, each being a small perturbation of the previous frame. This
enforces the continuity of features from frame to frame in the reconstructed
movie. For the static regularization terms, we chose the Second Moment
regularizer (Issaoun et al. 2019) and the Total Squared Variation regularizer
(described in Appendix A.1). The Second Moment regularization constrains the
spread of flux density in reconstructed images to a motivated region defined
by the user. In the case of the EHI system, although single frames can have
very limited coverage, the dense $uv$-coverage is produced in the whole range
of baselines considering time-averaged data. Hence source parameters, required
for the Second Moment regularization, can be obtained from the same
observational data by the time-averaged image reconstruction. The
implementation of the Second Moment regularizer into the dynamical imaging
method can grant the highest completeness, accuracy and fidelity of
reconstructed movies.
The dynamical regularizer $\mathcal{R}_{\Delta t}$ computes the similarity
between frames of the movie by calculating the summed difference of the
reconstructed flux density of pixels among all adjacent frames after blurring
using a circular Gaussian with standard deviation $\sigma_{\Delta t}$ (Johnson
et al. 2017). The selection of the parameter $\sigma_{\Delta t}$ was performed
for each temporal resolution independently. This parameter equals 0 for
observations with 1 $t_{\mathrm{G}}$ temporal resolution. The limit
$\sigma_{\Delta t}\to 0$ is appropriate when the expected motion between
consecutive frames is smaller than the finest resolution of reconstructed
features comparable to the nominal array spatial resolution (Johnson et al.
2017). To choose $\sigma_{\Delta t}$ for observations with 10 $t_{\mathrm{G}}$
temporal resolution, half the number of frames of the M87* model with 10
$t_{\mathrm{G}}$ interval between frames was used. For observations with 100
$t_{\mathrm{G}}$ temporal resolution, the full model with 10 $t_{\mathrm{G}}$
interval between frames was used considering the variability of the source
during the observation. The used system setup had 350 km orbital separation
and an observing frequency of 230 GHz. $\sigma_{\Delta t}$ was varied from 1
to 100 and from 1 to 1000 for reconstructions with averaging over 10 and 100
$t_{\mathrm{G}}$ per frame, respectively. According to image quality metrics
averaged over the duration of reconstructed movies, $\sigma_{\Delta t}$ was
selected to equal 10 and 100, respectively, in these two cases. (Figure 20).
The selection of weighting parameters for the dynamical imaging method was
performed by simulating observation of half the number of frames of the M87*
model with 1 $t_{\mathrm{G}}$ interval between frames. The used system setup
had 1000 km orbital separation for the widest representation of baselines on
each of the frames, the observing frequency was 230 GHz. The data term
weighting $\alpha_{data}$ was left default and equals 10. The weighting
parameter for dynamical regularization $\alpha_{\Delta t}$ was ranged from 10
to $10^{4}$. As $\alpha_{\Delta t}\to 0$ the dynamical imaging becomes
equivalent to snapshot imaging with corresponding static regularizers. Taking
$\alpha_{\Delta t}\to\infty$ leads to the reconstruction of the time-averaged
image (Johnson et al. 2017). According to averaged over the duration of
reconstructed movies image quality metrics, $\alpha_{\Delta t}$ was selected
to equal 500 (Figure 21). The weighting parameters for the Second Moment
regularizer $\alpha_{SM}$ and for the Total Squared Variation regularizer
$\alpha_{TSV}$ were ranged between 10 and $10^{8}$, and between 1 and
$10^{3}$, respectively. Figure 22 demonstrates image quality metrics averaged
over the duration of reconstructed movies for different combinations of these
parameters. Dependencies demonstrated by the NRMSE and the NXCORR are
identical, similar to the snapshot imaging method. The best reconstruction
quality corresponds to the smallest weight of the Total Squared Variation
regularization. The weighting of the Second Moment regularization produces
little effect on the quality of reconstruction with a slight improvement for
large values, as shown in Figure 23. Therefore, weighting parameters for the
Second Moment regularization $\alpha_{SM}$ and the Total Squared Variation
regularization $\alpha_{TSV}$ were selected to equal $10^{7}$ and 1,
respectively. The regularization parameters are the same for all dynamical
imaging reconstructions displayed in this paper.


* Zhao et al. (2023a) Zilong Zhao, Robert Birke, and Lydia Y. Chen. 2023a.  GDTS: GAN-Based Distributed Tabular Synthesizer. In _2023 IEEE 16th International Conference on Cloud Computing (CLOUD)_. 570–576.  https://doi.org/10.1109/CLOUD60044.2023.00078
  * Zhao et al. (2021a) Zilong Zhao, Robert Birke, Aditya Kunar, and Lydia Y Chen. 2021a.  Fed-tgan: Federated learning framework for synthesizing tabular data.  _arXiv:2108.07927_ (2021).
  * Zhao et al. (2023b) Zilong Zhao, Han Wu, Aad Van Moorsel, and Lydia Y Chen. 2023b.  Gtv: Generating tabular data via vertical federated learning.  _arXiv:2302.01706_ (2023).
  * Zhao et al. (2023c) Zhuang Zhao, Feng Yang, and Guirong Liang. 2023c.  Federated Learning Based on Diffusion Model to Cope with Non-IID Data. In _Chinese Conference on Pattern Recognition and Computer Vision (PRCV)_. Springer, 220–231.
  * Zhou et al. (2021) Xu Zhou, Xiaofeng Liu, Gongjin Lan, and Jian Wu. 2021.  Federated conditional generative adversarial nets imputation method for air quality missing data.  _Knowledge-Based Systems_ 228 (2021), 107261.  https://doi.org/10.1016/j.knosys.2021.107261
  * Zhou et al. (2023) Yuhao Zhou, Mingjia Shi, Yuanxi Li, Qing Ye, Yanan Sun, and Jiancheng Lv. 2023.  Communication-efficient Federated Learning with Single-Step Synthetic Features Compressor for Faster Convergence.  arXiv:2302.13562 [cs.LG]
  * Zhu et al. (2021b) Hangyu Zhu, Jinjin Xu, Shiqing Liu, and Yaochu Jin. 2021b.  Federated learning on non-IID data: A survey.  _Neurocomputing_ 465 (2021), 371–390.
  * Zhu et al. (2021a) Zhuangdi Zhu, Junyuan Hong, and Jiayu Zhou. 2021a.  Data-free knowledge distillation for heterogeneous federated learning. In _International conference on machine learning_. PMLR, 12878–12889.