© The Authors 2025

1 Introduction

Supermassive black holes (SMBH) are found or are expected to be found in the centers of most galaxies (Richstone et al. 1998). The accretion of matter onto the central object causes the regularly observed periods of activity that significantly impact the intergalactic and circumgalactic environments and the evolution of their host galaxy. Additionally, SMBHs are used as a laboratory to study the physics of gravity itself. In this context, a primary laboratory is the SMBH in our own Galaxy, Sagittarius A*, or Sgr A* (Goddi et al. 2017).

Sgr A* is one of the key targets of the Event Horizon Telescope (EHT), a worldwide network of radiotelescopes that employs very long baseline interferometry (VLBI) at millimeter and submillimeter wavelengths (Event Horizon Telescope Collaboration 2019a; Raymond et al. 2024). The unprecedented angular resolution achieved by the EHT enabled us to capture the first images of SMBHs, that is, the SMBH at the center of the giant elliptical galaxy M87 (Event Horizon Telescope Collaboration 2019b,a,c,d,e,f, 2023) and the SMBH at the center of our own galaxy, in total intensity (Event Horizon Telescope Collaboration 2022a,b,c,d,e,f) and in polarized light (Event Horizon Telescope Collaboration 2021a,b, 2024a,b). In both targets, EHT observations revealed a ring-like feature, which was interpreted as the theoretically predicted signature of the black hole shadow. In the case of M87*, a similar ring-like feature has subsequently also been detected with the Global Millimeter VLBI Array (GMVA) at 86 GHz (Lu et al. 2023; Kim et al. 2025). Multifrequency observations of the innermost regions surrounding the central SMBH are crucial for characterizing the physical conditions of the accreting plasma and for constraining models for black hole accretion (Mościbrodzka et al. 2014, 2016).

However, a multifrequency analysis like this remains limited for Sgr A* as a result of an additional effect: interstellar scattering. Interstellar scattering is caused by variations in the electron density along the line of sight to the target. It results in multipath propagation, which broadens radio images and pulse profiles of objects that are located in or behind the Galactic plane. In certain directions, the scattering is significantly stronger than what is predicted by large-scale Galactic electron distribution models (Taylor & Cordes 1993; Cordes & Lazio 2002). This pronounced scattering has long been linked to Galactic HII regions and supernova remnants (Litvak 1971; Little 1973). Being situated at the center of our Galaxy, Sgr A* is chiefly affected by intense scattering. The size of its radio image increases with λ² (e.g. Bower et al. 2004; Issaoun et al. 2019), consistent with predictions for a thin scattering screen (Blandford & Narayan 1985).

The interferometric phase information is strongly affected by the presence of scattering. As a consequence, scattering-induced distortions on the observed data can significantly degrade the quality of the images, particularly at frequencies below 86 GHz. In this line, several strategies have been developed to mitigate the scattering effects and improve the quality and robustness of the image (Fish et al. 2014; Johnson 2016; Kouroshnia et al. 2025). These algorithms characterize the scattering kernel and incorporate a functional within the imaging deconvolution process to correct for the associated distortions. First, Fish et al. (2014) proposed to apply a Wiener filter (Wiener 1949) to the scattering kernel to increase the signal-to-noise ratio (SNR). In contrast, Johnson (2016) proposed to find the source and screen discrete reconstructions that fit the data best. As a byproduct, this algorithm is able to recover an image of the scattering-corrected observed source and an approximation of the scattering screen. The latter holds significant astrophysical information on its own and supports efforts to detect transient phenomena toward the Galactic Center (Caleb et al. 2022; Beniamini et al. 2023). One of the main challenges in this context is the potential degeneracy between the scattering screen and the image. This led to the development of imaging algorithms based on multimodal image posteriors. In particular, we base our analysis on the multiobjective evolution proposed by Müller et al. (2023) and then was extended by Mus et al. (2024b,a).

We present in this work the design, implementation, and testing of novel algorithms to mitigate the impact of (potentially dynamic) scattering screens that affects VLBI observations by recovering the intrinsic black hole structure images at lower frequencies and the screen itself.

Compared to traditional noise-inflation and deblurring techniques, such as those used by Event Horizon Telescope Collaboration (2022c, 2024a) or even standard unimodal exploration Johnson (2016), our novel strategy effectively explores the family of local solutions. This is an essential feature for addressing the highly nonlinear and ill-posed problem of VLBI imaging and screen modeling. Unlike methods that are constrained to a specific screen model (and thus, to a fixed intrinsic source structure), the flexibility of our approach provides a mathematically rigorous solution to the modeling and mitigation problem.

This paper is structured as follows. We present the background theory on VLBI imaging and scattering theory in Sect. 2. In Sect. 3 we discuss the problem modeling and imaging algorithm, and we verify the reconstruction with a static and dynamic screen and a static source in Sect. 4. In Sect. 5, we discuss the potential of recovering the scattering screen and the black hole shadow after proper descattering at 86 GHz. We show the marginal distributions of the scattering screen regularizers in Sect. 6, and we finally discuss future prospects in Sect. 7.

2 VLBI and imaging background

2.1 VLBI measurements

In a VLBI array, multiple antennas observe the same source at the same time, and the recorded signals are then correlated. These correlated signals from the antennas (visibilities) are approximately the Fourier transform of the true sky brightness distribution at a spatial frequency determined by the baseline separating the antennas in an antenna-pair projected on the sky-plane. In particular, the visibility function 𝒱(u, v), representing the correlated signal between an antenna pair, is governed by the van Cittert-Zernike theorem (Thompson et al. 2017), $$\mathcal{V}(u,v)=\int \int I(l,m)e^{-2\pi i(lu+mv)}dldm.$$(1)

This equation expresses that the true sky brightness distribution, I(l, m), and the observed visibilities are related as a Fourier pair. Here, (u, v) coordinates are defined by the relative baselines between antennas in the plane of the sky and (l, m) direction cosines measured with respect to the (u, v) axes. However, it is crucial to recognize that in practical scenarios, the Fourier domain is not fully sampled, leading to sparse coverage in the (u, v)-plane, particularly in VLBI contexts. This incomplete sampling makes the process of reconstructing the sky brightness distribution from the observed visibilities an ill-posed inverse problem. Additionally, these visibilities are subject to calibration errors and thermal noise. To mitigate these challenges, it is essential to account for calibration effects. Direction-independent calibration errors are modeled through station-based gain factors, g_i. The observed visibilities at a given time t for an antenna pair i, j are thus expressed as $$V(i,j,t)=g_{i,t}{g}_{j,t}^{\ast }\mathcal{V}(i,j,t)+N(i,j,t),$$(2)

where g_i and g_j are the station-based gains, and N (i, j, t) represents the random noise associated with the baseline.

These difficulties are exasperated in sparse VLBI networks such as the EHT, which additionally operates at high radio frequencies, rendering the data calibration and imaging processes particularly challenging. All these difficulties have required the development of a wide range of novel algorithms for image reconstruction that we describe in detail in the next subsection.

2.2 Outline of the imaging

The imaging algorithms employed in the EHT can be roughly categorized into five groups: techniques based on Bayesian exploration (Broderick et al. 2020; Tiede 2022), the regularized maximum likelihood technique (RML) (Akiyama et al. 2017b,a; Chael et al. 2018), global search techniques by multiobjective evolution (Müller et al. 2023; Mus et al. 2024a), compressive sensing (Müller & Lobanov 2022, 2023a), and inverse modeling with novel developments, for instance, by Müller & Lobanov (2023b). The majority of these algorithms achieves a moderate degree of super-resolution, that is, they robustly recover structures smaller than the original resolution limit of CLEAN (Högbom 1974; Clark 1980; Cornwell 2008). This has proven to be crucial for a reliable reconstruction of structures at the spatial scales of the event horizon.

The RML method is one of the most frequently used imaging techniques in the EHT (Event Horizon Telescope Collaboration 2019d, 2021a, 2022c, 2023, 2024c,a). Multiple data-fidelity terms and regularization terms are balanced against each other by a relative weighting. This weighted sum is then minimized. The solution of this optimization problem is the model that describes the data better. The data term takes into account how well one solution is described by the data. Frequently used data terms include χ² metrics for observed data products, for example, the visibilities or closure products. The regularization terms in turn describe the plausibility of the recovered solution. Standard choices include the total variation (promoting piecewise smooth functions), the total squared variation (promoting smoothness), the 11-norm (promoting sparsity), the 12-norm (reducing the overfitting), or an entropy functional.

It is crucial to use the correct relative weighting of the data and regularization terms for RML methods. The EHT has applied a strategy of thorough surveying a sufficient subset of possible parameter combinations. Multiobjective optimization aims at presenting a different strategy. The concept of optimality of a weighted sum of objectives is replaced by the concept of Pareto optimality in a multiobjective problem formulation: A solution is called Pareto optimal when the further optimization of one objective (e.g., the data term) automatically has to worsen the scoring in another objective (e.g., in a regularization term). The set of all Pareto-optimal solutions is referred to as the Pareto front.

The algorithm MOEA/D¹ is based on the genetic evolution that approximates the Pareto front (Zhang & Li 2007; Li & Zhang 2009). It has been adapted for VLBI objectives by Müller et al. (2023), who also showed that the Pareto front in VLBI experiments has a special structure: It is separated into multiple clusters, which are interpreted as locally optimal modes of the potentially multimodal imaging problem. Furthermore, the cluster with the best solution has in practice been found to be the cluster that lies closest to the ideal point (see Müller et al. 2023, for more details on this construction). This algorithm is in active use for the imaging in extremely weak constrained settings or in problems in which internal degeneracy may cause a multimodal posterior. This includes the reconstruction of the polarization structure from the closure traces (Müller 2024) or the imaging of solar flares with the Spectrometer/Telescope for Imaging X-rays (STIX) instrument (Müller et al. 2024).

The variant MO-PSO² is much faster and more accurate than MOEA/D that was proposed by Mus et al. (2024a). Instead of reconstructing the full image structure by an evolutionary algorithm, the code navigates the Pareto front by an evolutionary algorithm (particle swarm optimization) and solves the scaled problems by fast L-BFGS-B (limited memory Broyden–Fletcher–Goldfarb–Shanno) minimization. The Pareto front is searched in order to minimize the distance of a Paretooptimal solution to the ideal point. Thus, instead of recovering the full Pareto front, MO-PSO directly provides its final best image reconstruction, for which it spends less time and fewer computational resources.

2.3 Scattering

In this section, we give an overview of the theory of and motivation for the thin-screen scattering. We focus on the relevant aspects related to the goals of our strategy. For further details, we refer to the reviews by Rickett (1990); Narayan (1992); Johnson & Gwinn (2015); Thompson et al. (2017).

When a source is observed through a medium with spatial variations in its refractive index, it becomes distorted. Refractive inhomogeneities cause different regions of an image associated with the source to be steered and focused differently while surface brightness is maintained (Born & Wolf 1980). In particular, radio-wave scattering in the ionized interstellar medium (ISM) is due to density inhomogeneities because the refractivity in a plasma is approximately proportional to the local electron density (Jonson 1999).

In general, scattering is often described as resulting from turbulent media localized in a single thin screen between the source and observer (see Bower et al. 2014; Dexter et al. 2017; Psaltis et al. 2018, for instance). We work throughout under the assumptions of one unique thin screen that can be approximated by n layers, although the number of screens in the line of sight is still an opened question (e.g. Dexter et al. 2017).

When r is a two-dimensional transverse coordinate on the screen, this screen introduces a stochastic position-dependent phase shift ϕ(r) of the incoming radiation, without altering the amplitude of the waves. If this screen follows a homogeneous Gaussian distribution, the scattering can be quantified in two complementary ways: by the structure function of the phase fluctuations, D_ϕ(r) = ⟨[ϕ(r + r′) − ϕ(r′)]²⟩, or by the power spectrum of the phase fluctuations, P_ϕ(q) ∝ |q|^–(α+2). We assumed α to be 1.38³. These approaches are related by a Fourier transform, P_ϕ(q) ∝ q² D_ϕ(q), with a prefactor that makes P_ϕ(q) dimensionless and independent of the observing wavelength (λ).

Diffractive and refractive scattering effects on images can be treated as distinct phenomena (Blandford & Narayan 1985). Diffractive effects dominate at small scales (approximately r₀) and are well approximated by their ensemble average, that is, the expectation value of the measured complex visibilities when averaging over many realizations of the scattered phases (Narayan & Goodman 1989; Goodman & Narayan 1989). This averaging effectively blurs the intrinsic image with a kernel G(r) (the seeing disk). This kernel is most naturally expressed in the visibility domain as $\tilde{G}(\mathbf{b})=e^{-\frac{1}{2}{D}_{\varphi }(\mathbf{b})}$, where b corresponds to the physical length of an interferometric baseline. Refractive effects can be modeled using a geometrical optics framework, where gradients of the large-scale refractive modes of the phase screen steer and focus the ensemble-average image.

The single-epoch scattered image I_a(r) is related to the unscattered image I_src(r) via Johnson & Narayan (2016); Johnson (2016) $$I_a(r)=I_{\mathrm{src}}(r)\ast G(r)+(\mathrm{\nabla }{\varphi }_r(r))\cdot [I_{\mathrm{src}}(r)\ast \mathrm{\nabla }G(r)].$$

In these expressions, ∇ represents the two-dimensional transverse gradient on the phase screen. The term I_src(r) * G(r), where the asterisk denotes spatial convolution, represents the blurring effect due to scattering and is referred to as the ensemble-average image, I_ea(r) = I_src(r) * G(r). The second term, (∇φ_r(r)) · [I_src(r) * ∇G(r)], accounts for additional distortions introduced by variations in the phase screen. This term captures the influence of the turbulent phase gradient across the scattering screen, further modulating the observed image through spatially dependent distortions.

For each image, r is a transverse coordinate at the distance of the scattering screen D (not the distance of the source D + R), so that the corresponding angular scales are $\theta =\frac{r}{D}$. For simplicity, we use ϕ(r) throughout to denote the refractive phase screen ϕ_r(r). We denote the forward operator applying the scattering with Φ for the remainder of this paper, that is, $$\mathrm{\Phi }:[I_{\mathrm{s}\mathrm{r}\mathrm{c}},\varphi ]\mapsto I_a.$$(3)

The Strehl ratio, 0 < S ≤ 1, is a commonly used metric in the optical literature to quantify image degradation due to scattering. It is defined as the ratio of the peak intensity of the measured point-spread function (PSF), which includes the effects of scattering, to the peak intensity of the ideal diffraction-limited PSF, which assumes no scattering.

This concept was extended to radio interferometry by Johnson (2016), and it is defined as $$S\simeq \frac{{\theta }_{uv}}{\sqrt{{\theta }_{uv}^2+{\theta }_{scatt}^2}},$$

where θ_uv is the full width at half maximum (FWHM) of the beam of the array, and θ_scatt is the FWHM of I_a. Johnson (2016) Fig. 1 shows the expected S for the different arrays.

3 Problem statement

In this section, we introduce the mathematical optimization formulation we used to model the imaging problem, along with the strategy we employed. We adopted a multiobjective optimization approach as presented by Mus et al. (2024a). For the sake of simplicity, we refer to the aforementioned paper for the mathematical formulation of the regularizers.

3.1 Modeling

Let F (x) := (f₁ (x), . . ., f_n (x)) be an n-vector of scalar objective functions f_i : ℝ^p → ℝ, and let $\mathcal{X}\subseteq {\mathbb{R}}_{+\mathrm{\infty }}^p$ denote the feasible set of decision vectors x. We then consider the multiobjective optimization problem

Following the strategy presented by Müller et al. (2023); Mus et al. (2024b,a), we solved Problem (MOP) by scalarization, that is, by rewriting the problem as the unconstrained minimization problem,

We recovered the images by using a wavelet dictionary Ψ, that is, we modeled the total intensity by I = Ψω_I⁴ (see Sect. 3.2 for more details). For the dynamic reconstruction of a movie in full polarization and a movie of the scattering screen, we therefore recovered the following quantities at every frame:

1. The wavelet coefficients describing the image in total intensity ω_I. They are related to the image at every frame by the relation I = Ψω_I.

2. The phase screen ϕ.

A movie is represented by a sequence of single individual images. In Prob. (MOP) x is therefore the vector, $$x=[{\omega }_I,{\varphi }^1,{\varphi }^2,{\varphi }^3,\dots ],$$(4)

where the superscript denotes the index of the frame in the movie.

The individual functionals are $$f_1:=R_{l1}\left(\mathrm{\Psi }{\omega }_I\right),$$(5) $$f_2:=R_{tv}\left(\mathrm{\Psi }{\omega }_I\right),$$(6) $$f_3:=R_{tsv}\left(\mathrm{\Psi }{\omega }_I\right),$$(7) $$f_4:=R_{l2}\left(\mathrm{\Psi }{\omega }_I\right),$$(8) $$f_5:=R_{flux}\left(\mathrm{\Psi }{\omega }_I\right),$$(9) $$f_6:=R_{entr}\left(\mathrm{\Psi }{\omega }_I\right),$$(10) $$f_7:=\sum _{j\in frames}{R}_{SO}\left({\varphi }^j\right),$$(11) $$\begin{array}{rl}& \sum _{f_8:=j\in frames}{S}_{vis}\left(\mathrm{\Phi }(\mathrm{\Psi }{\omega }_I,{\varphi }^j)\right)+S_{amp}\left(\mathrm{\Phi }(\mathrm{\Psi }{\omega }_I,{\varphi }^j)\right)\\ & +S_{clp}\left(\mathrm{\Phi }(\mathrm{\Psi }{\omega }_I,{\varphi }^j)\right)+S_{cla}\left(\mathrm{\Phi }(\mathrm{\Psi }{\omega }_I,{\varphi }^j)\right).\end{array}$$(12)

Here, f₁ to f₆ are the standard static imaging regularizers. For their exact form, we refer to Müller et al. (2023), f₇ is the stochastic optics functional (Johnson 2016), and f₈ is the data functional. These are χ² of the visibilities, amplitudes, closure phases, and closure amplitudes in total intensity. The data-fidelity terms are evaluated on the scattered movies, while the single regularizers are evaluated on the unscattered guess solutions.

By solving Prob. (MOP) via Prob. (MOP Scalar), we obtain the (static/dynamic) reconstruction that along with its associated power spectrum provides the best fit to the data. In Sect. 4.1 we model the power spectrum in detail.

3.2 Pipeline

The pipeline we used for the reconstruction employs naturebased optimization strategies: MOEA/D (Müller et al. 2023; Mus et al. 2024b) and MO-PSO (Mus et al. 2024a). The core concept of these strategies involves scalarizing the objective vector F and iteratively solving a series of subproblems. The two approaches offer significant advantages, which we discussed throughout. For a more in-depth understanding, we refer to the cited references.

3.2.1 MOEA/D

To explore the Pareto-front reconstructions, the interferometric reconstruction problem is modeled by $${\tilde{f}}_i=f_i+f_8,$$

where f_i are those defined from Eq. (5) to Eq. (11), and f₈ is the corresponding data term functional associated with the problem.

With this formulation, every functional ${\tilde{f}}_i$ gives a set of solutions that fit the data and are only affected by the regularizer f_i.

To solve the full MO, we applied a genetic algorithm. That is to say, we first gridded the multidimensional space of functions and drew a random sample of individuals (in our case, images or movies). Then, we applied a mutation and cross-over operation over the individuals producing new populations (Müller et al. 2023).

3.2.2 MO-PSO

For every subproblem, we applied MO-PSO (for instance Du & Swamy 2016), considering a swarm composed by the hyperparameters α_i that appear in Problem (MOP Scalar) and that act on every regularizer. In this way, we let the swarm converge to their optimal position (equivalently, to the optimal hyperparameter combination), and then we solved the imaging RML problem with the optimal weights.

In contrast to MOEA/D, the MO-PSO formulation allowed us to obtain the marginal contribution of every regularizer. In other words, instead of exploring the whole effects of each regularizer over the data terms, it seeks the best compromise among all the functionals.

By letting the population evolve long enough, individuals converge to a set of nondominate solutions (see Pardalos et al. 2017; Müller et al. 2023, for definition) that approximates the Pareto front. Therefore, an explicit dependence (manifold) between data terms and regularizers can be obtained.

4 Reconstruction of the scattering screen

In this section, we validate the algorithm performance at high frequencies, in particular, at 230 GHz. In this regime, scattering corruptions are minimal, thereby providing a clear framework for assessing the efficacy of the mitigation strategy. In the subsequent Sections 5 and Appendix 5.2, we investigate the algorithm behavior at lower frequencies, where the corruptions are stronger and the problem is more constrained.

4.1 Discretizsation of the stochastic optics scattering screen

First, we present the Johnson (2016) theoretical framework on stochastic optics, followed by the strategy we used to mitigate the scattering.

The Johnson (2016) framework represents ϕ(r) in the Fourier domain, with uncorrelated complex Gaussian variable components $$\tilde{\varphi }(\mathbf{q})=\int \varphi (\mathbf{r})e^{-i2\pi \mathbf{q}\cdot \mathbf{r}}d\mathbf{r}.$$

The time-averaged power spectrum is given by $${⟨|\tilde{\varphi }(\mathbf{q})|}^2⟩=A_{\varphi }Q(\mathbf{q}),$$

where A_ϕ is the screen area over which the Fourier transform is computed (see, e.g. Blandford & Narayan 1985).

The power spectrum associated with the observed visibilities is given by $$\begin{array}{rl}Q(\mathbf{q})=& 2^{\alpha }\pi \alpha \frac{\mathrm{\Gamma }(1+\alpha /2)}{\mathrm{\Gamma }(1-\alpha /2)}{\lambda }^{-2}{\left(r_{0,x}{r}_{0,y}\right)}^{-\alpha /2}\\ & \times {[\left(\frac{r_{0,x}}{r_{0,y}}\right)q_x^2+\left(\frac{r_{0,y}}{r_{0,x}}\right)q_y^2]}^{-(1+\alpha /2)},\end{array}$$(13)

where x, y are the coordinates aligned with the major and minor axes of the diffractive kernel. We considered the screen to be located at M = 0.43, as stated by Bower et al. (2014).

The scattering phase screen (in radians) was parameterized using an N × N grid of Fourier coefficients ${\tilde{\varphi }}_{o,s}$, $${\varphi }_{l,m}=\frac{1}{\sqrt{A_{\varphi }}}\sum _{o,s=0}^{N-1}\sqrt{Q(o,s)}{ϵ}_{o,s}{e}^{i2\pi (lo+ms)/N},$$

where ${ϵ}_{o,s}=\frac{{\tilde{\varphi }}_{o,s}}{\sqrt{Q(o,s)}}$. In this way, Johnson (2016) parameterized the unknown phase screen as a set of $\frac{N^2-1}{2}$ variables.

After ϵ_o,s was specified for a particular field of view, the corresponding scattering was computed for any desired observing wavelength. We adopted this framework to model Eq. (11).

Therefore, the mitigation can be summarized as follows: For the observation data, we tried to find the best image that fits it and whose power spectrum is close to Eq. (13). For every reconstruction, we computed its power spectrum and minimized along the χ² axis. We denote this functional as R_{S O}. Mathematically, $$R_{SO}:=\sum {\left|{ϵ}_{o,S}\right|}^2.$$

Johnson et al. (2017) solved the optimization problem

where R_entr, f₉ are the functionals described in Eqs. (10) and (12), respectively, and only one frame (j = 1) was fixed and was applied without wavelet coefficients⁵. There are some limitations, as listed below.

• The lack of multimodal exploration in the context of highly nonlinear nonconvex functionals and the use of a gradientbased approach limits the exploration and screen modeling. To overcome this limitation, the standard approach performs a survey of the different hyperparameters of the optimization problem.

• The inclusion of an additional hyperparameter increases the dimensionality of the problem, and this might reduce the efficiency of the survey strategy. This is translated into a much slower algorithm. In addition, the SO algorithm is slower than regular imaging because the screen is discretized and the power spectrum is computed, which complicates and increases the cost of the imaging problem. As a result, the associated optimization process converges more slowly or requires more time to complete the same number of iterations as in a standard imaging problem. This challenge is particularly pronounced for large datasets.

• The screen effects are complex. To understand these effects, a full marginalization would be better.

In contrast, the flexibility of multiobjective optimization enables the immediate integration of R_{S O} within complex frameworks. This allows the method to leverage several advantages. These include avoiding costly hyperparameter explorations and achieving a balanced optimization with other regularizers beyond R_entr. For example, MOEA/D facilitates exploration of the Pareto front and identification of different local minima, which is particularly crucial for very sparse observations that are weakly affected by refractive noise, such as those from the EHT or ngEHT (Raymond et al. 2021; Chavez et al. 2024) at 230 GHz. Similarly, MO-PSO offers a straightforward way to determine the marginal contribution of the regularizer, such as the screen effect on the data. In this way, the three limitations of the method proposed by Johnson (2016) were addressed.

4.2 Synthetic data generated at 230 GHz

To assess the performance of our modeling, we tested it on various synthetic structures created with eht-imaging library (Chael et al. 2018), with the assumption of static and dynamic screens at 230 GHz. The simulated datasets include different levels of realistic gain corruptions. In the following sections, we show the results.

There are no updates to the modeling of the screen itself in this paper, but we update the way in which the optimization problem associated with its reconstruction is solved. Our main improvement stems from the correct modeling of the multimodality of the problem, that is, of the degeneracies associated with the fact that different realizations of the scattering screen and the on-sky emission might produce similar visibilities. We identified three scenarios for which this strength can be exceptionally beneficial. First, the scattering effect at 230 GHz observations is weak even with a static source and a static screen, and it is expected to be degenerate with the noise realization. Second, when a speed is associated with the scattering, the dynamics of the screen might be misidentified with small-scale corruptions in the actual image. Finally, at 86 GHz, the screen obscures the image and significantly worsens the instrumental resolution. In this situation, small-scale structures at 20 μas scales that are typical for the EHT might be attributed to the effect of the scattering screen or the internal features of the emission. The test cases we discuss here address these three scenarios.

As a first and simplest example, we simulated a static ring with an asymmetric brightness distribution observed during the best time window (Fig. 1, top left panel). This time window spans approximately 100 minutes and presents the best uv coverage, isotropy, and (u, v) coverage fraction for Sgr A* in April 7, 2017, during the EHT campaign (more details can be found in Farah et al. 2022; Event Horizon Telescope Collaboration 2022c)⁶.

Gain corruptions of 4% and 10% were included in the data. We additionally corrupted the observations with several synthetic screens that were simulated using the StochasticOptics (Johnson 2016) module of the eht-imaging library. The bottom panel of Fig. 1 bottom shows the effects of the refractive noise on the longer baselines for the two levels of gain corruption. We demonstrate that we are able to recover the intrinsic descattered source and the scattering screen, and also that the screen prior is important in this poorly constrained problem.

In a second step, we increased the difficulty by simulating a population of different geometric sources that are affected by scattering screens with different velocities. We demonstrate its capability to recover the (dynamic) scattering screen even assuming a gain corruption of 10%. We would like to note that this scenario is only a first step toward constructing an algorithm that can adapt to real data. In EHT data, Sgr A* shows significant variability that stems from the evolution of the source itself and not from the dynamics of the scattering screen (e.g. Event Horizon Telescope Collaboration 2022a; Wielgus et al. 2022). Any algorithm would need to separate the effect of the source dynamics from the effect of the scattering and the internal screen variability. Correct modeling of the dynamics of Sgr A* is a challenging but independent problem (even without scattering), and a wide range of solutions has been proposed (Bouman et al. 2017; Event Horizon Telescope Collaboration 2022c; Müller & Lobanov 2023a; Roelofs et al. 2023; Mus & Martí-Vidal 2024; Mus et al. 2024b,a). We instead studied the simpler question whether we can correctly recover the speed of the screen for a static source that is scattered by a dynamically evolving screen. Despite the simplicity of the question, this yielded two important outcomes. First, the successful demonstration of this capability helped us to constrain the currently broad range of expected screen velocities (~50 km/s to ~200 km/s Gwinn et al. 1991; Reid et al. 2019). Second, it needs to be demonstrated before a framework is established that solves the more complex scenario in which the source and screen are both dynamic.

4.3 Recovering the screen: Algorithm

The strategy we used is summarized in the Algorithm 1. First, we defined a prior image and a prior screen. Throughout this paper, we understand the concept of a prior as the initial point of a minimization algorithm. For the case of 230 GHz data, we used a 60 μas⁷ symmetric Gaussian as a prior. Then, we computed the wavelet coefficients, Ψw_I, as described by Müller & Lobanov (2022), to help us overcome the sparsity of the data. These inputs were used to define a composite vector x that encapsulated the image and scattering information, with an overall dimension of 2npix² − 1. In the first optimization step, the algorithm applied MO-PSO, using the prior image as the starting point. This optimization was carried out over P₀ particles and K₀ iterations, and it yielded an intermediate solution x_img that was selected based on optimal regularization criteria.

The next stage involved redefining the prior as a combination of the optimized image x_img and the initial prior scattering screen. This updated prior served as the initialization for a second MO-PSO run, which aimed to refine the solution further by simultaneously optimizing the image and the scattering screen. This stage used P₁ particles and K₁ iterations, and it ultimately produced a final image with minimized scattering effects.

4.4 Results

Throughout this paper, we evaluate the results by using the normal cross-correlation⁸, or nxcorr metric (see Event Horizon Telescope Collaboration 2019d; Farah et al. 2022, for formal definition). In Table B.1 of the Appendix B, all the nxcorr values for all images can be found.

Fig. 1 Top left panel: uv coverage for Sgr A* in the best time window of April 7, 2017. The different colors represent different baselines (indicated in the legend). Top right panel: simulated ring in this coverage, with a brightness asymmetry. Bottom panel: log amplitude vs. uv distance for the static scattered asymmetric ring with different gain corruptions. Blue shows 4%, and red shows 10%.

4.4.1 Different gain corruptions

In these paragraphs, we assumed a screen of 50 km/s in the east direction during the whole observation.

Figure 2 presents the results of the Algorithm 1 for the cases of 4% and 10% gain corruptions. Since the screen at this high frequency is poorly constrained and the uv-coverage very sparse, a uniform screen (in contrast with the proposal of Johnson 2016) is not a valid starting point.

In both cases, the recovered scattering screen are compatible, and they show some similarities to the real screen in the sense that the positions of the brightest spots are approximately recovered and depressions are kept as well, as is visible in the projection along the y-axis of the recovered screen with the true screen (second and fourth row). However, we would like to note that none of the finer structure in the simulated screen are recovered correctly, and the localization of the peaks and troughs is not very accurate either. This is probably to be expected because of the challenges of an array such as the EHT and the small effect of the screen at 230 GHz. Moreover, since the effect of the screen is very small at 230 GHz, the resulting, descattered image is similar to the on-sky image, that is, the on-sky image already shows a ring structure.

4.4.2 Independence of the screen prior

As discussed above, a physically feasible screen prior is needed to recover the screen in this case with high degrees of freedom and very sparse uv-coverage. To ensure that the algorithm is not prior biased, we present two tests in which the asymmetric ring was scattered with the same true screen, but with two different initial screens.

Figure 3 shows the two cases assuming systematic errors of 10%. The reconstructed screen is indeed clearly different. However, the main features, such as the brightest areas and the depressions, are recovered. In particular, the prior screen (even if it is not similar to the real screen) can help the algorithm to converge to a better screen. Nevertheless, the refractive noise is not the dominating noise and is only weakly constrained, and the differences on the intrinsic structure reconstructions are therefore negligible. This difference is also represented in the slightly different values of the nxcorr shown in Table B.1.

Fig. 2 Scattering reconstruction using MO-PSO of the static ring with an asymmetric brightness distribution with thermal noise and a gain corruption of 4% (top row), and with thermal noise and a gain corruption of 10% (bottom row). From left to right: simulated blurred ring, ground-truth screen, recovered intrinsic structure, recovered screen, and prior screen. The screens were convolved with the instrumental beam. The subpanels below the second and fourth panels depict the projection on the y-axis of the brightness distribution along the right ascension.

4.4.3 Different screen velocities

Next, we attempted a more challenging test case. We tried to recover the speed of the scattering screen (on a static source). As discussed above, this was an artificial experiment because in practice, Sgr A* also evolves dynamically. However, we defer the full analysis of the source dynamics together with a moving screen to a consecutive work. In the dynamic regime, it is essential to investigate whether the temporal evolution of the scattering screen impacts each frame of the observation. To do this, we simulated a screen at different velocities that affect static geometric sources based on the best time-window coverage and corrupted by thermal noise. In particular, we simulated two screens at two different velocities, one velocity of 50 km/s based on observations of Gwinn et al. (1991), and the other of 200 km/s, Reid et al. (2019). For a scattering screen velocity of 50 km/s, the effects of the kinematics of the screen over a period of 100 min are minor, but a velocity of 200 km/s produces noticeable effects. In any dynamic study, the speed of the screen therefore needs to be properly constrained. To study these effects of the speed, we solved MO-PSO + SO for every keyframe.

The weak refractive noise effect on the data implies that the speed velocity does not affect the intrinsic structure recovery. Fig. 4 depicts a comparison of the same scattered source (first panel), a crescent, that is affected by the same screen at different velocities (upper row 50 km/s, lower row 200 km/s).

To illustrate the screen evolution, we projected the pixel brightness distribution along the y-axis and plotted it as a function of time, following a similar approach as was described by Mus & Martí-Vidal (2024). This is called the cylindrical plot. The third and fourth panels display the recovered intrinsic source structure and the corresponding recovered evolution of the scattering screen, respectively. The simulated screen moved along the x-axis, so that the projection was done orthogonal to the direction of motion. The cylindrical plots (second and fourth column) can be understood as the logical extension of the projections shown in Fig. 2, that is, every column of the plotted matrix corresponds to the projection along the y-axis of the true and recovered scattering screen at the respective frame in time.

We recall from our discussion above that while the reconstruction of any structural feature finer than roughly 20 μas is not recovered, even in the static case. However, the broad overall structure of the approximate positioning of troughs and peaks across the field of view can be identified. Fig. 4 now demonstrates that we can even track this evolution in time, and we can hence effectively measure the speed of the screen from the reconstruction by the gradient in the cylindrical plots.

Fig. 3 Comparison of the scattering reconstructions of the ring with an asymmetric brightness with different screen priors. The screens were convolved with the instrumental beam.

Fig. 4 Recovery of the scattering screen velocity for a screen that affects a static crescent. Top row: true scattering screen evolves at 50 km/s (center left panel), while the recovered velocity is 51.4 km/s (right panel). Bottom row: true scattering screen evolves at 200 km/s (center left panel), while the recovered velocity is 193.8 km/s (right panel). The slight deviation from the exact velocity arises from numerical errors in the fitting procedure. The left and center right panels show the corresponding reconstructed crescent for each case.

Fig. 5 Dynamic scattering reconstruction of the intrinsic source and the evolving screen equivalent to that in Fig. 4. Top row refers to the disk model and bottom row to the ring model.

4.4.4 Different intrinsic structures

Figure 5 presents the results of the same experiment with different observed sources. It illustrates that the geometry of the sources significantly affects the recovery of the (dynamic) screen. In the case of the disk (bottom row), the cylindrical plot for the recovered screen is noticeably less noisy than that of the ring (top row). The null in the uv-plane is affected by the screen; it is harder to solve the imaging problem. Nevertheless, we succeeded in correctly measuring the speed of the screen, with an error smaller than 10 km/s.

5 Prospects of imaging the Sgr A* ring at 86 GHz

While at 230 GHz our method successfully recovered the overall dynamics and morphology of the screen, it was unable to capture small-scale structures (as predicted and shown in Johnson 2016) because the refractive noise was poorly constrained. In contrast, this noise at 86 GHz is expected to be more constrained on longer baselines, which might allow us to recover finer structural details.

In this section, we explore the prospects of recovering ring structures in the GMVA frequency regime. To achieve this, we exploited the multimodality capabilities of MOEA/D. We simulated observations using the 2025 GMVA configuration, incorporating Atacama Large Millimeter/submillimeter Array (ALMA) stations, including the Northern Extended Millimeter Array (NOEMA) and a double bandwidth.

5.1 Static 86 GHz. Including gain errors

Using ngehtsim (Pesce 2022), we simulated static ring observations with a shadow diameter of ~50 μas, emulating that reported by the EHT (Event Horizon Telescope Collaboration 2022a). We assumed bad weather conditions to decrease the SNR, and we included thermal errors, as well as gain errors of 4% in phase and amplitude. The aim was to verify whether the intrinsic ring structure could be recovered. Figure 6 shows the uv coverage corresponding to the simulated array (left panel) and the SNR versus uv-distance (right panel).

In the range of 86 GHz, scattering effects dominate the noise, and so the image is strongly affected. Its Strehl ratio is S ~ 0.6 (see Fig. 1 from Johnson 2016). That is to say, refractive noise at longer baselines is much more constrained. This leaves fewer degrees of freedom to reconstruct the screen. In any case, a careful mitigation strategy is crucial.

The complex structure of the screen and the steep brightness gradients that can occur between nearby pixels mean that genetic algorithms may struggle to converge rapidly to an optimal solution. To address this issue, we ran MO-PSO using as a starting point the cluster that contained the most local minima or the most probable solution, which is also the most probable solution based on the RML method (although any other cluster could also be chosen; see Müller et al. 2023; Mus & Martí-Vidal 2024, for further details). By starting in a neighbor of a local minimum, we exploited gradient-based algorithms to quickly find a solution. Moreover, MO-PSO ensures convergence to the solution with the optimal regularizer contributions.

The strategy we used is described below and summarized in Algorithm 2. We initialized the MOEA/D algorithm with a population of 340 individuals and 51 neighbors. For the entropy prior, we assumed a uniform Gaussian of 60, μas. The first generation of individuals was initialized with a Gaussian whose beam parameters were taken from the least-squares (LSQ) Gaussian fit reported by Issaoun et al. (2019), with values of (120, μas, 100, μas) and a position angle (PA) of 90 degrees. We flagged data with an SNR below 10, which is equivalent to approximately 1.25% of the maximum SNR. We considered the data terms of MOEA/D to be closure quantities (phases and logamplitudes). We let the population evolve over 1000 generations. After the population evolution was complete, we clustered the 340 individuals (solutions) by similarity, requiring 10% similar features to belong to the same cluster. We selected a representative image from the cluster that contained more solutions (the most probable local minimum), we self-calibrated N times, running MOEA/D in every iteration, and we used this image as a starting point and let the new population (again composed of 340 individuals and 51 neighbors) evolve over 500 generations⁹. When it was required, we recentered the source to the (0,0) since the relative position information can be lost due to the use of closures. Fig. C.1 shows the set of clusters for three iterations. The clusters surrounded by a green box have the highest nxcorr with respect to the model. The red box marks the cluster that contained more solutions, and the blue box shows the cluster with better hyperparameter combination. The solution we obtained after the last self-calibration iteration was used as starting point to solve MO-PSO. The faster convergence of MO-PSO allowed us to obtain a more accurate solution faster. Fig. 7 shows the final source and screen reconstruction. The intrinsic structure has a nxcorr of 0.97, and the nxcorr of the recovered screen is 0.81 (blurred) and 0.88 (nonblurred). These metrics demonstrate how successful the algorithm was in recovering a screen in a more constrained case.

In the next section, we apply MOEA/D with different starting points and report the probabilities of obtaining a ring solution under optimal conditions using a priori calibrated data. Our findings indicate that for certain initial conditions, no convergence to a ring solution is reached, and the probability for others, such as a Gaussian starting point, remains low, but a ring can indeed be recovered.

Fig. 6 uv-coverage for the GMVA + ALMA array observing Sgr A* at 86 GHz (left panel) and the signal-to-noise ratio in log scale vs. uv distance (right panel). The bad weather conditions and the gain corruption of 4% cause a a low SNR even in the case of ALMA baselines.

5.2 Static 86 GHz. A priori calibrated data with thermal noise

In this section, we explore the chances of obtaining a ring structure using different starting points when the a priori calibrated data are only corrupted with thermal noise. The goal of this section is to demonstrate the importance of the global exploration and that the assumed starting geometry may condition the problem, so that no ring can be recovered even in favorable situations. The strategy we used is summarized in Algorithm 3. The three starting points were the Gaussian observed by Issaoun et al. (2019) using LSQ before applying SO (major FWHM of 228 μas, minor FWHM of 143 μas and PA of 86 grad), a disk of 60 μas, and a ring with the expected shadow size seen by the EHT of 50 mas (Event Horizon Telescope Collaboration 2022a). Fig. A.1 shows the images (top row) and their amplitude versus uv distance in log scale (bottom row).

For every case, we explored ~360 local minima and let the population to evolve during 10 000 generations. We allowed each individual to interact (crossover) with 15% of the population. The data-term functional was the χ² of the visibilities.

Figure C.2 presents the Pareto fronts of the solutions originating from three distinct starting points, the Gaussian, disk, and ring, each clustered using a uniform similarity threshold of 10%. The title of each cluster indicates the percentage of solutions contained within it, along with the χ² value, which quantifies the discrepancy between the true unscattered ring and the reconstructed image.

The different number of clusters contained in every starting point is the first thing to note. The starting point conditions the diversity of the recovered solutions. Then, the ring starting from a Gaussian can be recovered, but not the ring from a disk. This means that the optimization algorithm tends to find disk structures faster than a ring. Therefore, the optimization needs to start from a distant point. Another interesting fact is χ² is systematically lower in all the clusters with ring-like structures. Therefore, even though the remaining clusters may approximate the Pareto front (i.e., are local minima, valid images), those with a ring present the best fit to the data.

Table 1 summarizes the results of the MOEA/D exploration with the different starting points. For the Gaussian starting point, 25.45% of the solutions are ring-like (~85 local minima). This means that the probability of finding a local solution that represents a ring is just ~0.25. Therefore, unimodal optimization methods (based on gradients or Hessian) are likely to find nonring solutions. In contrast, using a ring with the expected size of Sgr A* as the initial point increases the probability to almost to 40%. Although considering a ring as initial point maybe a strong bias, we can claim that with two completely different starting points, MOEA/D is able to recover a ring.

Thus, we conclude that the choice of prior complicates the accurate recovery of the true source structure. In this context, standard local searches struggle to explore alternative candidates, and they therefore are suboptimal strategies for identifying lower-probability solutions in a highly ill-posed problem.

Fig. 7 Scattered ring observed at 86 GHz (first panel) with a gain corruption of 4%. The second panel shows the simulated screen that affects the source. The third panel depicts the recovered solution using the Algorithm 2. The white internal circle shows the diameter of the shadow of the simulated ring (50 μas), and the dashed green circle is 1.6 times the size of the shadow, as predicted by Lu et al. (2023); Kim et al. (2024) in M87*. The fourth panel shows the recovered screen. The screens are blurred to the GMVA resolution. nxcorr can be found in Table B.1.

5.3 Conclusions

We have investigated the feasibility of recovering the Sgr A* ring structure at 86 GHz in the presence of refractive scattering and various observational challenges. The refractive noise at this lower frequency, although dominant, is more strongly constrained on longer baselines than at 230 GHz. This allowed our global optimization schemes (MOEA/D followed by MO-PSO) to resolve finer structural features. Even under adverse conditions, such as a low signal-to-noise ratio and gain errors, the proposed multistep strategy robustly recovers the intrinsic source geometry and a meaningful approximation of the scattering screen. The final nxcorr values reach 0.97 for the source.

Moreover, we demonstrated that the choice of the initial image geometry significantly influences whether the optimization converges to a ring solution. In idealized simulations with a priori calibrated data, starting from a Gaussian prior yields a probability of approximately 25% of recovering a ring, whereas beginning with a ring prior increases this probability to nearly 40%. By contrast, a disk prior never converged to a ring solution in our tests. These findings underscore the importance of genuinely global exploration methods and of the use of diverse initial conditions when strongly scattered sources such as Sgr A* are imaged. Overall, a comprehensive multimodal approach is essential to mitigate the ill-posed nature of VLBI data at 86 GHz and to reliably recover ring-like structures in the presence of refractive scattering.

Fig. 8 Marginal contribution of the hyperparameter associated with Eq. (11) for observations at 86 and 230 GHz.

6 Marginal contribution of the SO regularizer

To conclude the analysis, we examined the marginal contribution of the Regularizer (11) at frequencies of 230 GHz and 86 GHz. We followed a strategy similar to that outlined by Mus et al. (2024a). We ran MO-PSO with the same initial point, number of iterations, simulated screen, and number of particles, but varied the initialization of these particles by following a uniform distribution. This approach allowed us to assess the impact of the regularizer for different refractive noise conditions. Using MO-PSO, we can determine the relative importance of a given regularizer on the final reconstruction. This information constrains the possible values of the regularizer. The lower the variance in the values, the more constrained and significant the regularizer for the reconstruction.

Figure 8 shows the distribution of the hyperparameter values for both frequencies. As expected, we note a broader distribution for 230 GHz, meaning that the regularizer is less constrained. Table 2 summarizes the main features. The standard deviation (std) for 230 GHz is indeed larger than for 86 GHz.

7 Future work

This work opens several promising directions for future investigation. First, the nxcorr metric was found to be suboptimal for validating our results, as elaborated in Appendix B. While the main features of the reconstructions may be recovered, the metric values often do not accurately reflect the visual quality of the outcomes. In future works, we will aim to find a better and appropriate metric that successfully reflects the fidelity of the reconstructions.

Second, we considered Sgr A* as a static source. However, the reconstruction of the scattering screen in Sgr A* may present an additional challenge due to the intrinsic variability on minute timescales of this source (Event Horizon Telescope Collaboration 2022c). A wide range of movie-making algorithms that account for time-domain correlations have been proposed (Bouman et al. 2017; Müller & Lobanov 2022; Roelofs et al. 2023; Mus & Martí-Vidal 2024; Mus et al. 2024b,a). With the recent studies of the dynamics of Sgr A* at horizon scales (EHTC, in prep.), future research will focus on incorporating Sgr A* as a dynamic source and a dynamic scattering screen to achieve reconstructions of real data by better capturing the time-dependent behavior of the system.

Third, we currently investigate how the intrinsic geometry of the initial conditions influences the multimodality exploration. We observed that starting the exploration with certain geometric models, for example, disk, prevented us from recovering the intrinsic ring, while using a Gaussian as an initial point allowed us to recover ring structures for ~25% of the cases (see Sect 5.2 for further details). Another case we aim to understand is why the reconstruction of the scattering screen appears to be noisier when applied to a ring structure, as indicated in Fig. 4.

Fourth, we intend to apply this strategy to real data at 86 GHz that are observed with the new sites and capabilities of the GMVA+ALMA array, and we will try to obtain the real ring of Sgr A*.

Last, we currently examine the effects of scattering on hotspot tracking. This might provide deeper insights into the complex dynamics of the system.

8 Summary and conclusions

Scattering remains one of the most significant challenges in radio astronomy observations of the Galactic Plane, particularly at low frequencies. It complicates the image reconstruction by introducing substantial stochastic distortions that can severely degrade the quality of the recovered images (Johnson & Narayan 2016) and corrupt signals. This poses considerable difficulties for tasks such as pulsar searches (Narayan 1992).

Several methods have been developed to mitigate its effects (Fish et al. 2014; Johnson et al. 2018) and even to reconstruct the affecting screen. However, these methods are constrained by their tendency to find only local minima and are limited to a restricted set of regularizers due to constraints of the computational power.

We introduced a new strategy to mitigate scattering and reconstruct the screen that integrates the modeling of stochastic optics within the framework of a multiobjective optimization. Our method provides a mathematically rigorous solution to the problem without being restricted to a specific screen model or intrinsic source structure. This flexibility is particularly valuable to model the highly nonlinear and ill-posed nature of the VLBI imaging and screen. In contrast to traditional noise-inflation and deblurring techniques (Event Horizon Telescope Collaboration 2022c, 2024a) and to standard unimodal exploration (Johnson 2016), our strategy thoroughly explores the family of local solutions. This crucial aspect is frequently over-looked in conventional approaches. In this way, this approach solves the challenges faced by previous methods, such as (1) the high computational cost, (2) the inefficiency due to the increased dimensionality of the problem by the inclusion of a new regularizer, and (3) the loss of information about the effects of scattering mitigation on the image. By solving Prob. (MOP Scalar) using nature-inspired optimization techniques (specifically, genetic algorithms and particle swarm optimization), we effectively explored the multimodal nature of the problem.

We demonstrated that this novel strategy successfully mitigates scattering in very sparse EHT observations of approximately 100 minutes at 230 GHz, where scattering noise is weak and poorly constrained, and in high SNR GMVA+ALMA observations at 86 GHz that lasted about 12 hours, where scattering is stronger and better constrained.

We applied our algorithm to synthetic EHT observations to test various screen velocities, intrinsic geometric sources, levels of gain corruption, and different priors. In all cases, the recovered unscattered source and the corrupting scattering screens were well reconstructed. Notably, even under these poorly constrained conditions, our method succeeded in recovering the screen velocity. This is an improvement over earlier approaches, and it is crucial for constraining the expected range (Gwinn et al. 1991; Reid et al. 2019). However, in this high-frequency domain, the prior screen proved crucial for aiding convergence, as expected because the nature of the problem is poorly constrained.

Additionally, we modeled a ring structure and applied the scattering kernel predicted by Bower et al. (2015); Psaltis et al. (2018) for Sgr A* at 86 GHz. This introduced various gain corruptions, and we consistently assumed adverse weather conditions. Our objective was to recover the intrinsic ring structure despite these challenges. Leveraging the significant exploration capabilities of genetic algorithms, we were able to successfully reconstruct the intrinsic ring and the scattering screen.

To this end, we proposed two different strategies based on the level of gain corruption under two representative gain corruption conditions. One strategy was designed to only explore multimodality, and the second strategy selected the most probable solution and self-calibrated to it. The strategy iterated this procedure as many times as imposed by the user.

We analyzed the different local minima we recovered and determined, many of which had a ring morphology based on different initial configurations (Gaussian, ring, and disk). Our findings indicate that, for instance, starting from a disk configuration presents the highest likelihood of failure to accurately recover the ring. Because of the more constrained nature of this problem, we were able to recover a better image of the scattering screen.

We finally explored the marginal contribution of the stochastic optics regularizer. Our findings confirmed that at lower frequencies, where long baselines are dominated by refractive noise, the constraints are tighter. This reduces the degrees of freedom of the regularizer.

In conclusion, multiobjective optimization together with exploration algorithms is a tool that might enhance our ability to see the ring of Sgr A* in real GMVA+ALMA observations at 86 GHz. It has not been detected at this frequency so far.

Data availability

Our imaging pipeline and our software is included in the second release of MrBeam¹⁰. Our software makes use of the publicly available eht-imaging library (Chael et al. 2016, 2018), regpy (Regpy 2019), MrBeam (Müller & Lobanov 2022, 2023b,a; Müller et al. 2023) and pyswarms¹¹ packages.

Acknowledgements

This work was supported by the Italian Ministry of University and Research (MUR) -– Project CUP F53D23001260001, funded by the European Union – NextGenerationEU. T.T. acknowledge the “Center of Excellence Severo Ochoa” grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033 awarded to the Instituto de Astrofísica de Andalucía. H.M., G.Y. and A.L. acknowledge the M2FINDERS project funded by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (Grant Agreement No. 101018682) and by the MICINN Research Project PID2019-108995GB-C22. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The authors gratefully acknowledge Y. Kovalev for his insightful comments and valuable suggestions, and the anonymous referee for their constructive feedback.

Appendix A Additional figures

Fig. A.1 To explore the multimodality of the solutions at 86 GHz we set the starting point for MOEA/D to be a Gaussian (corresponding to LSQ solution in Issaoun et al. 2019), disk and ring. Top row of the figure shows the images, bottom row, the amplitude vs uv-distance including the simulated scattered ring (black dots).

Appendix B Comparison metrics

In this section, we present the nxcorr values obtained from all the reconstructions performed. Table B.1 contains the nxcorr values of the reconstructions of the intrinsic descattered source and the recovered (blurred and non-blurred) scattering screen appearing in Figures 2, 3 and 7. It is important to note that the nxcorr tends to improve at lower frequencies, which is an expected outcome due to the increased constraints at those frequencies. However, it should also be noted that the nxcorr may not be the most suitable metric for evaluation. The nxcorr assumes a Pearson correlation, which may not be optimal for non-normally distributed data and because of the given the high variance in pixel brightness values (ranging between −3000 to 3000 radians in some instances), which significantly reduces the accuracy of the comparisons. $${\rho }_{\mathrm{N}\mathrm{X}}(X,Y)=\frac{1}{N}\sum _i\frac{(X_i-⟨X⟩)(Y_i-⟨Y⟩)}{{\sigma }_X{\sigma }_Y}.$$(B.1)

Eq. (B.1) provides the definition of nxcorr, used to compare two images, X and Y. Due to the high variance in this context, the application of the mean and normalization by the standard deviation (σ) may be less appropriate. As discussed in Sect. 7, we are developing a more sophisticated metric to address these limitations.

Appendix C Pareto fronts

In this appendix we show the Pareto fronts of the different sources discussed in the main text.

Fig. C.1 Pareto front for the intrinsic source at 86 GHz across each self-calibration iteration. Arranged from left to right, the panels correspond to the first, second, and third iterations, respectively. In the initial iteration, the ring structure is absent from the solution set. A hint of the ring appears in the second iteration, meaning the null in the visibilites is presented in the set of local minima. In addition, we note how the diversity of solutions increases. The red box encapsulates the most frequently repeated (most probable) solution, the blue box highlights the solution closest to the ideal, and the green box identifies the one with the lowest χ2.

Fig. C.2 Pareto front for the intrinsic source at 86 GHz. (a) Gaussian (LSQ Issaoun et al. 2019) starting point, (b) disk starting point, (c) ring starting point. The red box encapsulates the most frequently repeated (most probable) solution, the blue box highlights the solution closest to the ideal, and the green box identifies the one with the lowest χ2.

References

All Tables

All Figures

	Fig. 1 Top left panel: uv coverage for Sgr A* in the best time window of April 7, 2017. The different colors represent different baselines (indicated in the legend). Top right panel: simulated ring in this coverage, with a brightness asymmetry. Bottom panel: log amplitude vs. uv distance for the static scattered asymmetric ring with different gain corruptions. Blue shows 4%, and red shows 10%.
In the text

Fig. 2 Scattering reconstruction using MO-PSO of the static ring with an asymmetric brightness distribution with thermal noise and a gain corruption of 4% (top row), and with thermal noise and a gain corruption of 10% (bottom row). From left to right: simulated blurred ring, ground-truth screen, recovered intrinsic structure, recovered screen, and prior screen. The screens were convolved with the instrumental beam. The subpanels below the second and fourth panels depict the projection on the y-axis of the brightness distribution along the right ascension.

In the text

	Fig. 3 Comparison of the scattering reconstructions of the ring with an asymmetric brightness with different screen priors. The screens were convolved with the instrumental beam.
In the text

Fig. 4 Recovery of the scattering screen velocity for a screen that affects a static crescent. Top row: true scattering screen evolves at 50 km/s (center left panel), while the recovered velocity is 51.4 km/s (right panel). Bottom row: true scattering screen evolves at 200 km/s (center left panel), while the recovered velocity is 193.8 km/s (right panel). The slight deviation from the exact velocity arises from numerical errors in the fitting procedure. The left and center right panels show the corresponding reconstructed crescent for each case.

In the text

	Fig. 5 Dynamic scattering reconstruction of the intrinsic source and the evolving screen equivalent to that in Fig. 4. Top row refers to the disk model and bottom row to the ring model.
In the text

	Fig. 6 uv-coverage for the GMVA + ALMA array observing Sgr A* at 86 GHz (left panel) and the signal-to-noise ratio in log scale vs. uv distance (right panel). The bad weather conditions and the gain corruption of 4% cause a a low SNR even in the case of ALMA baselines.
In the text

Fig. 7 Scattered ring observed at 86 GHz (first panel) with a gain corruption of 4%. The second panel shows the simulated screen that affects the source. The third panel depicts the recovered solution using the Algorithm 2. The white internal circle shows the diameter of the shadow of the simulated ring (50 μas), and the dashed green circle is 1.6 times the size of the shadow, as predicted by Lu et al. (2023); Kim et al. (2024) in M87*. The fourth panel shows the recovered screen. The screens are blurred to the GMVA resolution. nxcorr can be found in Table B.1.

In the text

	Fig. 8 Marginal contribution of the hyperparameter associated with Eq. (11) for observations at 86 and 230 GHz.
In the text

	Fig. A.1 To explore the multimodality of the solutions at 86 GHz we set the starting point for MOEA/D to be a Gaussian (corresponding to LSQ solution in Issaoun et al. 2019), disk and ring. Top row of the figure shows the images, bottom row, the amplitude vs uv-distance including the simulated scattered ring (black dots).
In the text

Fig. C.1 Pareto front for the intrinsic source at 86 GHz across each self-calibration iteration. Arranged from left to right, the panels correspond to the first, second, and third iterations, respectively. In the initial iteration, the ring structure is absent from the solution set. A hint of the ring appears in the second iteration, meaning the null in the visibilites is presented in the set of local minima. In addition, we note how the diversity of solutions increases. The red box encapsulates the most frequently repeated (most probable) solution, the blue box highlights the solution closest to the ideal, and the green box identifies the one with the lowest χ2.

In the text

	Fig. C.2 Pareto front for the intrinsic source at 86 GHz. (a) Gaussian (LSQ Issaoun et al. 2019) starting point, (b) disk starting point, (c) ring starting point. The red box encapsulates the most frequently repeated (most probable) solution, the blue box highlights the solution closest to the ideal, and the green box identifies the one with the lowest χ2.
In the text