Imaging the Schwarzschild-radius-scale Structure of M87 with the Event Horizon Telescope using Sparse Modeling

TL;DR

This paper introduces a new sparse modeling imaging technique for interferometry data, enabling super-resolution imaging of black hole shadows with the EHT, surpassing traditional resolution limits.

Contribution

The paper presents a novel regularization-based imaging method using cross validation to optimize resolution, demonstrating super-resolution capabilities for black hole imaging.

Findings

Achieves 20-30% of diffraction limit resolution

Successfully reconstructs black hole shadow in simulations

Enables super-resolution imaging with EHT data

Abstract

We propose a new imaging technique for radio and optical/infrared interferometry. The proposed technique reconstructs the image from the visibility amplitude and closure phase, which are standard data products of short-millimeter very long baseline interferometers such as the Event Horizon Telescope (EHT) and optical/infrared interferometers, by utilizing two regularization functions: the $\ell_1$-norm and total variation (TV) of the brightness distribution. In the proposed method, optimal regularization parameters, which represent the sparseness and effective spatial resolution of the image, are derived from data themselves using cross validation (CV). As an application of this technique, we present simulated observations of M87 with the EHT based on four physically motivated models. We confirm that $\ell_1$+TV regularization can achieve an optimal resolution of $\sim 20-30$% of the diffraction limit $\lambda/D_{\rm max}$, which is the nominal spatial resolution of a radio interferometer. With the proposed technique, the EHT can robustly and reasonably achieve super-resolution sufficient to clearly resolve the black hole shadow. These results make it promising for the EHT to provide an unprecedented view of the event-horizon-scale structure in the vicinity of the super-massive black hole in M87 and also the Galactic center Sgr A*.

Figures

Tables

Table 1. Table 1: Stations used in simulated observations

Telescope	SEFD (Jy)
Phased ALMA	100
Phased SMA and JCMT	4000
LMT	1400
IRAM 30m	1400
NOEMA single dish	5200
ARO/SMT	11000

Equations

$$V(u,v)=\int dxdyI(x,y)e^{-i2\pi (ux+vy)}.$$

$$\mathbf{V}=\mathbf{F}\mathbf{I}.$$

$$V_j={\bar{V}}_j{e}^{i{\varphi }_j},$$

$${\psi }_m^{123}={\tilde{\varphi }}_j^{12}+{\tilde{\varphi }}_k^{23}+{\tilde{\varphi }}_l^{31}={\varphi }_j^{12}+{\varphi }_k^{23}+{\varphi }_l^{31},$$

$$\underset{\mathbf{I}}{\min }C(\mathbf{I})\mathrm{subject}\mathrm{to}{I}_{i,j}\ge 0.$$

$$C(\mathbf{I})={\chi }^2(\mathbf{I})+{\eta }_l||\mathbf{I}|{|}_1+{\eta }_t||\mathbf{I}|{|}_{\mathrm{tv}},$$

$$||\mathbf{I}|{|}_p={\left(\sum _i\sum _j|I_{i,j}{|}^p\right)}^{\frac{1}{p}}(\mathrm{for}p>0),$$

$${\chi }^2(\mathbf{I})=||\bar{\mathbf{V}}-\mathbf{A}(\mathbf{F}\mathbf{I})|{|}_2^2+||\mathbf{\Psi }-\mathbf{B}(\mathbf{F}\mathbf{I})|{|}_2^2,$$

$${\eta }_{\mathrm{l}}\equiv {\tilde{\eta }}_{\mathrm{l}}(N_{\mathrm{amp}}+N_{\mathrm{cphase}})/\max (\bar{\mathbf{V}}),$$

$$||\mathbf{I}|{|}_{\mathrm{tv}}=\sum _i\sum _j\sqrt{|I_{i+1,j}-I_{i,j}{|}^2+|I_{i,j+1}-I_{i,j}{|}^2}.$$

$${\eta }_{\mathrm{t}}\equiv {\tilde{\eta }}_{\mathrm{t}}(N_{\mathrm{amp}}+N_{\mathrm{cphase}})/4\max (\bar{\mathbf{V}}),$$

$$I_{i}^{\rm filtered}=\left\{\begin{array}[]{ll}0&(|I_{i}|<t)\\ I_{i}&({\rm otherwise})\end{array}\right.,$$

$$\mathrm{NRMSE}(\mathbf{I},\mathbf{K})=\sqrt{\frac{{\sum }_i|I_i-K_i{|}^2}{{\sum }_i|K_i{|}^2}},$$

$$||\mathbf{I}|{|}_{\mathrm{t}{\mathrm{v}}^2}=\sum _i\sum _j\left(|I_{i+1,j}-I_{i,j}{|}^2+|I_{i,j+1}-I_{i,j}{|}^2\right).$$

Full text

accretion, accretion disks — black hole physics — Galaxy: center — submillimeter: general — techniques: interferometric

††journal: ApJ††software: MAPS (http://www.haystack.mit.edu/ast/arrays/maps/), CASA

1 Introduction

Supermassive black holes (SMBHs) reside in the majority of the galactic nuclei in the universe. In a subset of such galaxies, accretion drives a highly energetic active galactic nucleus (AGN) often associated with powerful jets. Understanding the nature of these systems has been a central quest in astronomy and astrophysics. The SMBHs at the centers of our galaxy (Sgr A*) and the giant elliptical galaxy M87 provide unprecedented opportunities to directly image the innermost regions close to the central black hole, since the angular size of the event horizon is the largest among known black holes. The angular size of the Schwarzschild radius ($R_{s}$) is $\sim 10$ $\mu$as for Sgr A* for a distance of 8.3 kpc and a mass of $4.3\times 10^{6}$ $M_{\odot}$ (e.g. Chatzopoulos et al., 2015), and $\sim 3-7$ $\mu$as for M87 with a distance for 16.7 Mpc (Blakeslee et al., 2009) and a mass of $3-6\times 10^{9}$ $M_{\odot}$ (e.g. Gebhardt et al., 2011; Walsh et al., 2013). The apparent diameter of the dark silhouette of the black hole is $\sqrt{27}$ $R_{s}$ for the non-rotating black hole. It corresponds to $\sim 52$ $\mu$as for Sgr A* and $\sim 16-36$ $\mu$as for M87, which changes by only 4% with the black-hole spin and viewing orientation (Bardeen, 1973).

Very long baseline interferometric (VLBI) observations at short/sub-millimeter wavelengths ($\lambda\lesssim 1.3$ mm, $\nu\gtrsim 230$ GHz) can achieve a spatial resolution of a few tens of microarcseconds and therefore are expected to resolve event-horizon-scale structures, including the shadow of SMBHs. Indeed, recent significant progress on 1.3 mm VLBI observations with the Event Horizon Telescope (EHT; Doeleman et al., 2009) has succeeded in resolving compact structures of a few $R_{s}$ in the vicinity of the SMBHs in both Sgr A* and M87 (e.g. Doeleman et al., 2008, 2012; Fish et al., 2011, 2016; Akiyama et al., 2015; Johnson et al., 2015). Direct imaging of these scales will be accessible in the next few years with technical developments and the addition of new (sub)millimeter telescopes such as the Atacama Large Submillimeter/millimeter Array (ALMA) to the EHT (e.g. Fish et al., 2013).

Regardless of the observing wavelength, angular resolution (often referred to as “beam size” in radio astronomy and “diffraction limit” in optical astronomy) is simply given by $\theta\approx\lambda/D_{\rm max}$, where $\lambda$ and $D_{\rm max}$ are the observed wavelength and the diameter of the telescope (or the longest baseline length for the radio interferometer), respectively. A practical limit for a ground-based, 1.3 mm VLBI array like the EHT is $\sim 25$ $\mu$as ($=1.3$ mm/$10000$ km), which is comparable to the radius of the black hole shadow in M87 and Sgr A*. Hence, imaging techniques with good imaging fidelity at a spatial resolution smaller than $\lambda/D_{\rm max}$ would be desirable, particularly for future EHT observations of M87 and Sgr A*.

The imaging problem of interferometry is formulated as an under-determined linear problem when reconstructing the image from full-complex visibilities that are Fourier components of the source image. In the context of compressed sensing (also known as “compressive sensing”), it has been revealed that an ill-posed linear problem may be solved accurately if the underling solution vector is sparse (Donoho, 2006; Candes & Tao, 2006). Since then, many imaging methods have been applied to radio interferometry (see Garsden et al., 2015, and references therein). We call these approaches “sparse modeling” since they utilize the sparsity of the ground truth.

In Honma et al. (2014), we applied LASSO (Least Absolute Shrinkage and Selection Operator; Tibshirani, 1996), a technique of sparse modeling, to interferometric imaging. LASSO solves under-determined ill-posed problems by utilizing the $\ell_{1}$-norm (see §2.2 for details). Minimizing the $\ell_{1}$-norm of the solution reduces the number of non-zero parameters in the solution, equivalent to choosing a sparse solution. The philosophy of LASSO is therefore similar to that of the traditional CLEAN technique (Högbom, 1974), which favors sparsity in the reconstructed image and has been independently developed as Matching Pursuit (Mallat & Zhang, 1993) in statistical mathematics for sparse reconstruction. In Honma et al. (2014), we found that LASSO can potentially reconstruct structure $\sim 4$ times finer than $\lambda/D_{\rm max}$. Indeed, it works well for imaging the black hole shadow for M87 with the EHT in simulations.

Our previous work (Honma et al., 2014) has three relevant issues. The first issue is reconstructing the image only from the visibility amplitudes and closure phases (see §2.1), which have been the standard data products of EHT observations (Lu et al., 2012, 2013; Akiyama et al., 2015; Wagner et al., 2015; Fish et al., 2016) and optical/infrared interferometry. The algorithm of Honma et al. (2014) is applicable only for full-complex visibilities, which are the usual data products from longer-wavelength radio interferometers. We have recently developed a fast and computationally cheap method to retrieve the visibility phases from closure phases (PRECL; Ikeda et al., 2016), which can reconstruct the black hole shadow of M87 combined with LASSO in simulated EHT observations. However, since the phase reconstruction in PRECL adopts a different prior on visibilities than LASSO, the resultant image may not be optimized well in terms of $\ell_{1}$-norm minimization and sparse image reconstruction. Another potential approach is to solve for the image and visibility phases with $\ell_{1}$-norm regularization simultaneously, enabling us to reconstruct the image with full advantage of the regularization function.

The second issue is that the $\ell_{1}$-norm regularization might not provide a unique solution and/or could reconstruct an image that is too sparse image if the number of pixels with non-zero brightness is not small enough compared to the number of pixels. This violates a critical assumption in techniques with $\ell_{1}$-norm regularization that the solution (i.e. the true image) should be sparse. Such a situation may occur for an extended source or also even for a compact source if the imaging pixel size is set to be much smaller than the size of the emission structure. Pioneering work has made use of transforms to wavelet or curvelet bases, in which the image can be represented sparsely (e.g. Li et al., 2011; Carrillo et al., 2012, 2014; Garsden et al., 2015; Dabbech et al., 2015). As another strategy to resolve this potential issue, in Honma et al. (2014), we proposed to add another regularization, Total Variation (TV; e.g. Rudin et al., 1992), which is another popular regularization in sparse modeling. TV is a good indicator for sparsity of the image in its gradient domain instead of the image domain (see §2.2 for details) and it has been applied to astronomical imaging (e.g. Wiaux et al., 2010; McEwen & Wiaux, 2011; Carrillo et al., 2012, 2014; Uemura et al., 2015; Chael et al., 2016) that includes optical interferometric imaging without the visibility phases (e.g. MiRA; Thiébaut 2008; and also see Thiébaut 2013 for a review). TV regularization generally favors a smooth image (i.e. with larger effective resolution) but with a sharp edge, in contrast with maximum entropy methods (MEM; e.g. Narayan & Nityananda, 1986), which favor a smooth edge (see e.g. Uemura et al., 2015, for comparison between TV and MEM). Inclusion of TV regularization enables reconstruction of an extended image while preserving sharp emission features preferred by $\ell_{1}$-norm regularization, thereby extending the class of objects where sparse modeling is applicable. Indeed, regularization with both the $\ell_{1}$-norm and TV has been shown to be effective for imaging polarization with full complex visibilities in our recent work (Akiyama et al., 2017).

An important detail is the determination of regularization parameters (e.g. weights on regularization functions), which is common in the vast majority of existing techniques. Since one can not know the true image of the source a priori, one should evaluate goodness-of-fitting and select appropriate regularization parameters from the data themselves. In well-posed problems, one can use statistical quantities considering residuals between data and models as well as model complexity to avoid over-fitting, such as reduced $\chi^{2}$, the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) using the degrees of freedom to constrain model complexity. However, for ill-posed problems like interferometric imaging, degrees of freedom can not be rigorously defined, preventing the use of such statistical quantities.

In this paper, we propose a new technique to reconstruct images from interferometric data using sparse modeling. The proposed technique directly solves the image from visibility amplitudes and closure phases. In addition to the $\ell_{1}$-norm (LASSO), we also utilize another new regularization term, TV, so that a high-fidelity image will be obtained even with a small pixel size and/or for extended sources. Furthermore, we propose a method to determine optimal regularization parameters with cross validation (CV; see §2.3), which can be applied to many existing imaging techniques. As an example, in this paper, we applied our new technique to data obtained from simulated observations of M87 with the array of the EHT expected in Spring 2017.

2 The Proposed Method

2.1 A brief introduction of the closure phase

A goal of radio and optical/infrared interferometry is to obtain the brightness distribution $I(x,y)$ of a target source at a wavelength $\lambda$ or a frequency $\nu$, where $(x,y)$ is a sky coordinate relative to a reference position so called the phase-tracking center. The observed quantity is a complex function called visibility $V(u,v)$, which is related to $I(x,y)$ by two-dimensional Fourier transform given by

[TABLE]

Here, the spatial frequency $(u,v)$ corresponds to the baseline vector (in units of the observing wavelength $\lambda$) between two antennas projected to the tangent plane of the celestial sphere at the phase-tracking center.

Observed visibilities are discrete quantities, and the sky image can be approximated by a pixellated version where the pixel size is much smaller than the nominal resolution of the interferometer. The image can therefore be represented as a discrete vector ${\bf I}$, related to the Stokes visibilities ${\bf V}$ by a discrete Fourier transform ${\bf F}$:

[TABLE]

The sampling of visibilities is almost always incomplete. Since the number of visibility samples ${\bf V}$ is smaller than the number of pixels in the image, solving the above equation for the image ${\bf I}$ is an ill-posed problem.

Here, we consider that the complex visibility $V_{j}$ is obtained from observation(s) with multiple antennas. Let us define its phase and amplitude as $\phi_{j}$ and $\bar{V}_{j}$, respectively, denoted as follows

[TABLE]

where $j$ is the index of the measurement. Each measurement corresponds to a point $(u_{j},v_{j})$ in $(u,v)$-plane and recorded at time $t_{j}$. In actual observations, some instrumental effects and the atmospheric turbulence primary from the troposphere induce the antenna-based errors in the visibility phase, leading that the observed phase $\tilde{\phi_{j}}$ is offset from the true phase $\phi_{j}$ of the true image. In particular, this is a serious problem in VLBI observations performed at different sites (see Thompson et al., 2001).

However, the robust interferometric phase information can be obtained through the measurements of the closure phase, defined as a combination of triple phases on a closed triangle of baselines recorded at the same time. It is known that the closure phase is free from antenna-based phase errors (Jennison, 1958), which can be seen from the following definition of the closure phase,

[TABLE]

where $m$ is the index of the closure phase, and upper numbers ($1,2,3$) mean the index of stations involved in the closure phase or the visibility phase. The closure phase is also known as a phase term of the triple product of visibilities on closed baselines recorded at the same time, $V_{j}^{12}V_{k}^{23}V_{l}^{31}$, known as the bi-spectrum111Data products of visibility amplitudes and closure phases are also sometimes named as “bi-spectrum” in literature (e.g. Buscher, 1994). In this paper, we strictly distinguish them.. Closure phases have been used to calibrate visibility phases in VLBI observations (e.g. Rogers et al., 1974).

In short/sub-millimeter VLBI or optical/infrared interferometry, the stochastic atmospheric turbulence in the troposphere over each station induces a rapid phase rotation in the visibility, making it difficult to calibrate or even measure the visibility phase reliably (e.g., see Rogers et al., 1995; Thiébaut, 2013). Thus, image reconstruction using more robust closure phases, free from station-based phase errors, is useful for interferometric imaging with such interferometers.

2.2 Image Reconstruction from Visibility Amplitudes and Closure Phases

In this paper, we propose a method to solve the two-dimensional image ${\bf I}=\left\{I_{i,j}\right\}$ by the following equations:

[TABLE]

The cost function $C({\bf I})$ is defined as

[TABLE]

where $||{\bf I}||_{p}$ is $l_{p}$-norm of the vector ${\bf I}$ given by

[TABLE]

and $||{\bf I}||_{\rm tv}$ indicates an operator of TV.

The first term of Eq.(6) is the traditional $\chi^{2}$ term representing the deviation between the reconstructed image and observational data (i.e. the visibility amplitude ${\bf\bar{V}}=\{|V_{j}|\}$ and closure phase ${\bf\Psi}=\{\psi_{m}\}$), defined by

[TABLE]

where ${\bf A}$ and ${\bf B}$ indicate operators to calculate the visibility amplitude and closure phase, respectively. Deviations between the model and observational data are normalized with the errors. This form of the residual sum of squares (RSS) is originally proposed in the Bi-spectrum Maximum Entropy Method (BSMEM; Buscher, 1994) and also for modeling EHT data (e.g. Lu et al., 2012, 2013). Note that it could be replaced to a RSS term for bi-spectra (e.g. Bouman et al., 2015; Chael et al., 2016).

The second term represents LASSO-like regularization using the $\ell_{1}$-norm. Under the non-negative condition, $\ell_{1}$-norm is equivalent to the total flux. $\eta_{\rm l}$ is the regularization parameter for LASSO, adjusting the degree of sparsity by changing the weight of the $\ell_{1}$-norm penalty. In general, a large $\eta_{\rm l}$ prefers a solution with very few non-zero components, while $\eta_{\rm l}=0$ introduces no sparsity. In this paper, we use the normalized regularization parameter $\tilde{\eta}_{\rm l}$ defined by

[TABLE]

which is less affected by the number of visibility amplitude and closure phase data points, $N_{\rm amp}$ and $N_{\rm cphase}$, respectively, and also by the total flux density of the target source.

The third term is the TV regularization, defined by the sum of all differences of the brightness between adjacent image pixels. In this paper, we adopt a typical form for two-dimensional images (Rudin et al., 1992) that has been used in astronomical imaging, (e.g. Wiaux et al., 2010; McEwen & Wiaux, 2011; Uemura et al., 2015; Chael et al., 2016), defined as

[TABLE]

TV is a good indicator of image sparsity in its gradient domain instead of the image domain. TV is highly affected by the effective spatial resolution of the image. For instance, an image with a small TV value has blocks of image pixels whose brightness are similar, since the image is sparse in the gradient domain. The size of such blocks is equivalent to the effective spatial resolution, getting smaller for images with higher TV values. Thus, the regularization parameter $\eta_{\rm t}$ adjusts the effective spatial resolution of the reconstructed image. In general, a larger (smaller) $\eta_{\rm t}$ prefers smoother (finer) distribution of power with less (higher) discreteness, leading to larger (smaller) angular resolution. In the present work, we use the normalized regularization parameter $\tilde{\eta}_{{\rm t}}$ defined by

[TABLE]

similar to the LASSO term. A factor of 4 is based on a property of TV that takes a difference in the brightness to all four directions at each pixel. Note that a major difference to maximum entropy methods, which also favor smooth images, is that TV regularization has a strong advantage in edge-preserving; strong TV regularization favors a piecewise smooth structure, but with clear and often sharp boundaries between non-emitting and emitting regions.

The problem described in Eq.(5, 6) is non-linear minimum optimization. In this work, we adopt a non-linear programming algorithm L-BFGS-B (Byrd et al., 1995; Zhu et al., 1997) that is an iterative method for solving bound-constrained nonlinear optimization problems. L-BFGS-B is one of the quasi-Newton methods that approximates the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm using a limited amount of computer memory. In L-BFGS-B, the cost function and its gradient are used to determine the next model parameters at each iterative process. We approximately set partial derivatives to 0 at non-differentiable points for both the $\ell_{1}$-norm and TV. The partial derivatives of $\chi^{2}$ are calculated numerically with central differences. We use the latest Fortran implementation of L-BFGS-B (L-BFGS-B v3.0; Morales & Nocedal, 2011). We note that the problem is non-convex as other imaging techniques using closure quantities (e.g. Buscher, 1994; Thiébaut, 2008; Bouman et al., 2015; Chael et al., 2016), and a global solution is generally not guaranteed.

2.3 Determination of Imaging Parameters

In the proposed method, the most important tuning parameters are the regularization parameters for the $\ell_{1}$-norm ($\tilde{\eta}_{\rm l}$) and TV ($\tilde{\eta}_{\rm t}$), which determine the sparseness and effective spatial resolution of the image, respectively. Smaller regularization parameters generally favor images with larger numbers of non-zero image pixels and more complex image structure, which could give better $\chi^{2}$ values by over-fitting. On the other hand, large regularization parameters provide images that are too simple and that do not fit the data well. To determine optimal parameters, we need to evaluate the goodness-of-fit using Occam’s razor to prevent over-fitting.

In this work, we adopt cross validation (CV) to evaluate goodness-of-fit. CV is a measure of the relative quality of the models for a given set of data. CV checks how the model will generalize to an independent data set by using separate datasets for fitting the model and for testing the fitted model. CV consists of three steps: (1) randomly partitioning a sample of data into complementary subsets, (2) performing the model fitting on one subset (called the training set), and (3) validating the analysis on the other subset (called the validation set). To reduce variability, multiple rounds of cross-validation are performed using different partitions, and the validation results are averaged over the rounds. If the regularization parameters are too small, the established model from the training set would be over-fitted and too complicated, resulting in a large deviation in the validation set. On the other hand, if the regularization parameters are too large, the established model would be too simple and not well-fitted to the training set, also resulting in a large deviation in the validation set. Thus, reasonable parameters can be estimated by finding a parameter set that minimizes deviations (e.g. $\chi^{2}$) of the validation set.

In this work, we adopted 10-fold CV for evaluating the goodness-of-fit. The original data were randomly222In this work, subsamples are obtained with a uniform probability regardless of baselines following the most basic style of the CV. However, there could be more effective way of choosing subsamples for interferometric imaging. The optimum partitioning for CV is in the scope of our next studies in near future. partitioned into 10 equal-sized subsamples. 9 subsamples were used in the image reconstruction as the training set, and the remaining single subsample was used as the validation set for testing the model using $\chi^{2}$. We repeated the procedure by changing the subsample for validation data 10 times, until all subsamples were used for both training and validation. The $\chi^{2}$ values of the validation data were averaged and then used to determine optimal tuning parameters.

An important advantage of this method compared with previously proposed methods is that it is applicable to any type of regularization functions and also imaging with multiple regularization functions. For instance, Carrillo et al. (2012, 2014) and subsequent work solve images by utilizing the $\ell_{1}$-norm on wavelet-transformed image or TV regularization alone. In this case, the parameter can be uniquely determined from the $\ell_{2}$-norm of the estimated uncertainties on observational data (see Carrillo et al., 2012, for details). However, it is not straightforward to extend the idea for the problems with multiple regularization functions. For another example, Garsden et al. (2015) proposes another heuristic method to determine the regularization parameters on $\ell_{1}$-regularization on the wavelet/curvelet-transformed image by estimating its noise level on each scale, which is successful. However, the method would not work for all types of regularization functions. On the other hand, CV is a general technique that can be applied to imaging with any other regularization functions or any combination of them in principle, which include MEM (e.g. Buscher, 1994; Chael et al., 2016) and patch priors (e.g. Bouman et al., 2015). This advantage is particularly important for Sgr A*, which needs to involve a regularization function to mitigate the interstellar scattering effects (scattring optics; Johnson, 2016) in addition to the general regularization function(s) for imaging.

A relevant disadvantage of this method is its computational cost, since $n$-fold CV requires to reconstruct $(n+1)$ images for each set of regularization parameters. Recently, an accurate approximation of CV, which can be derived from a single imaging on full data set for each parameter set, has been proposed for imaging from full complex visibilities with $\ell_{1}$+TV regularizations (Obuchi & Kabashima, 2016; Obuchi et al., 2016). Future development of such heuristic approximations for other types of data and regularization functions could overcome this issue.

3 Imaging Simulations

3.1 Physically Motivated Models

In this paper, we adopt four physical models previously proposed for 1.3 mm emission on event-horizon scales.

The first model is a simple, but qualitatively correct, force-free jet model (hereafter BL09) in the magnetically dominated regime presented in Broderick & Loeb (2009) and Lu et al. (2014). We adopted a model image presented in Akiyama et al. (2015), which is based on the model parameters fitted to the results of 1.3 mm observations with the EHT in Doeleman et al. (2012) and the SED of M87 (Broderick et al. in preparation). The approaching jet is predominant for this model (see Akiyama et al., 2015, for more details).

The second and third models are based on results of GRMHD simulations presented in Dexter et al. (2012). We used the representative models DJ1 and J2, which are based on the same GRMHD simulation but with different energy and spatial distributions for radio-emitting leptons. The dominant emission region is the accretion flow in DJ1 and the counter jet in J2 illuminating the last photon orbit in J2. We adopt model images in Akiyama et al. (2015), where the position angle of the large-scale jet for models is adjusted to $-70^{\circ}$ inferred for M87 (e.g. Hada et al., 2011).

The last model is based on results of GRMHD simulations presented in Mościbrodzka et al. (2016), which models M87 core emission as radiation produced by the jet sheath. We use the image averaged for $\sim$ 3 months for our simulation (hereafter M16). The image has its dominant contribution from the counter jet illuminating the last photon orbit similar to J2 of Dexter et al. (2012), but the M16 model assumes energy distributions of leptons quite different from J2. We rotate the original model image of Mościbrodzka et al. (2016) to adjust the position angle of the large-scale jet to $-70^{\circ}$.

3.2 Simulated Observations

We simulate observations with the EHT at 1.3 mm (230 GHz) using the MAPS (MIT Array Performance Simulator) package333http://www.haystack.mit.edu/ast/arrays/maps/ based on the above models. The simulated observations are performed for the array expected to comprise in Spring 2017.

We assume an array consisting of stations at 6 different sites: a phased array of the Submillimeter Array (SMA) antennas and the James Clerk Maxwell Telescope (JCMT) on Mauna Kea in Hawaii; the Arizona Radio Observatory’s Submillimeter Telescope (ARO/SMT) on Mt. Graham in Arizona; the Large Millimeter Telescope (LMT) on Sierra Negra, Mexico; a phased array of the Atacama Large Millimeter/submillimeter Array (ALMA) in the Atacama desert, Chile; the Institut de Radioastronomie Millimétrique (IRAM) 30m telescope on Pico Veleta, Spain; and a single dish telescope of the Northern Extended Millimeter Array (NOEMA) in France. We adopt the system equivalent flux density (SEFD) of each station shown in Table 1 based on the proposer’s guide of 1-mm VLBI observations in ALMA Cycle 4.

The simulations assume a bandwidth of 3.5 GHz for Stokes $I$, which is half of the standard setting in ALMA Cycle 4444https://science.nrao.edu/observing/call-for-proposals/1mm-vlbi-cycle4/. We assume a correlation efficiency of 0.7, including a quantization efficiency of 0.88 for 2-bit sampling and other potential losses such as bandpass effects and pointing errors.

We simulate observations as a series of 5-minutes scans with a cadence of 20 minutes over a GST range of 13-0 hour, corresponding to the timerange when M87 can be observed by ALMA or LMT at an elevation greater than 20∘. ALMA and LMT are sensitive stations near the middle of the east-west extent of the array, and they may be important anchor stations for fringe detection. This provides an observational efficiency of 25% in time, expected for VLBI observations with ALMA in 2017. Data are integrated for the duration of each scan (i.e. 5 min) following previous EHT observations (e.g. Doeleman et al., 2012; Akiyama et al., 2015). Figure 1 shows the resultant $uv$-coverage of simulated observations. The maximum baseline length of observations is 7.2 G$\lambda$, corresponding to $\lambda/D=28.5$ $\rm\mu$as.

We note that the conditions of our simulation are much worse than previous simulations in Lu et al. (2014) in terms of the baseline sensitivity, $uv$-coverage, angular resolution (i.e. the maximum baseline length) and the exposure time of observations. Nevertheless, our simulated conditions are much closer to the observational conditions in Spring 2017.

3.3 Imaging

We reconstruct images from simulated data-sets based on the method described in §2.1. We adopt a field of view (FOV) of 200 $\mu$as gridded by 100 pixels in both the RA and Dec directions for all models, giving a pixel size of $\sim$1.6 ${\rm\mu}$as corresponding to a physical scale of $\sim 0.21$ $R_{s}$. Images are reconstructed at 4 regularization parameters for both $\tilde{\eta}_{\rm l}$ and $\tilde{\eta}_{t}$, ranging as $10^{-1},10^{0},...,10^{+2}$. As a result, we obtain $4\times 4=16$ images for each model.

Since the problem described in §2.1 is non-convex, our algorithm may be trapped in a local minimum and therefore may end up at an initial-condition dependent solution, similar to other algorithms using techniques minimizing non-convex functions with gradient descent methods. To avoid this, we start reconstructing images at $\tilde{\eta}_{\rm l}=\tilde{\eta}_{\rm t}=10^{2}$, which is expected to derive the simplest image among parameters we adopt, assuming a point source as initial images. This is then used as the initial image at other values of the regularization parameters.

At each parameter, we do iterations until achieving convergence or 1000 iterations, and then filter the output image with a hard thresholding defined by

[TABLE]

where $t$ is a threshold. We repeat this process until the normalized root mean square error (NRMSE) between the previous and latest filtered images becomes smaller than $1\%$. NRMSE is defined by (e.g. Chael et al., 2016)

[TABLE]

where ${\bf I}$ and ${\bf K}$ are the image to be evaluated and the reference image, respectively. We adopt the latest non-filtered image as the final product. Although this procedure makes the computational time longer, we found this works very well for avoiding convergence at some local minima. In this paper, we set 10% of the peak flux as a threshold $t$. For the simulated observational data in this paper, it takes typically about several to ten minutes on a standard desktop computer with six intel core-i7 CPU cores to reconstruct an image at each set of two regularization parameters.

We evaluate the goodness-of-fit for each image and then selected the best-fit images with 10-fold CV as described in §2.3. The quality of the reconstructed images is evaluated with the NRMSE. Since the all model images have finer resolutions with narrower FOVs than the reconstructed images, we calculated these metrics as follows. First, we adjusted the pixel size of the reconstructed image to that of the model image with bi-cubic spline interpolation. Second, we adjusted the position offsets between these two images so that the positions of their centers of mass coincide555Previous works (e.g. Lu et al., 2014; Fish et al., 2014) derived position offsets between the model and reconstructed images by taking cross correlations of the two images. However, we found that the position offsets derived from cross correlations induce an additional error in the metrics due to errors in position offsets. We found that the center of mass for the image is a better indicator of the position offsets than the cross-correlation minimum., because absolute positions cannot be defined from visibility amplitudes and closure phases alone. Finally, the metrics were evaluated. In addition to NRMSE, we also measure structural dissimilarity (Wang et al., 2004) between the model and reconstructed images using the DSSIM metric adopted work by Lu et al. (2014) and Fish et al. (2014). Since both metrics show similar trends, we show only the behavior of the NRMSE in the figures that follow.

3.3.1 Imaging with the Cotton-Schwab CLEAN

For evaluating the performance of our techniques, we also reconstructed images with the most widely-used Cotton-Schwab CLEAN (henceforth CS-CLEAN; Schwab, 1984) implemented in the Common Astronomy Software Applications (CASA) package666https://casa.nrao.edu/ with uniform weighting. Since CLEAN requires complex visibilities, we adopted the simulated complex visibilities with thermal noises. We set a gain of 0.1 and a threshold of 0.08 mJy beam*-1*, comparable to the image sensitivity of simulated observations. Since the fast Fourier transform is often used in CLEAN, a very small FOV can require a grid size in $uv$-plane that is too large, which could cause additional deconvolution errors. Hence, we set 1024 pixels with the same pixel size in each axis for the entire map, and put a CLEAN box in the central 100$\times$100 pixels to put CLEAN components in the same region as other techniques. We use the model image instead of the CLEAN map for calculating metrics, since the residual map, which is generally added to the CLEAN map, cannot be calculated for the proposed method.

4 Results

4.1 The best-fit images

We show the best-fit images selected with CV in Figure 2. A clear shadow feature is well reproduced for the counter-jet- and accretion-flow-dominated models (J2, M16 and DJ1). This demonstrates that the EHT will achieve effective sufficient spatial resolution to image the black hole shadow of M87 if the mass-loading radius of the jet is not too large (Broderick & Loeb, 2009; Lu et al., 2014; Akiyama et al., 2015).

Figure 3 shows the NRMSE metric for reconstructed images. The black curve labeled “Model” shows the NRMSE calculated when the model image is convolved with a circular Gaussian beam with a full width at half maximum (FWHM) as shown on the abscissa and compared against the original (unconvolved) model image. The Model curve effectively quantifies the best-case scenario in which the differences from the original input are due solely to a loss of resolution, not to errors in reconstructing the image. Figure 3 also show the NRMSE of each of the reconstructed images convolved with circular Gaussian beams.

Figure 3 clearly shows that closure-phase imaging with $\ell_{1}$+TV regularization works well compared to CS-CLEAN in particular at finer resolutions, despite the fact that CS-CLEAN uses full complex visibilities with more information and higher SNRs than closure phases. Both CS-CLEAN and $\ell_{1}$+TV images achieve similar NRMSEs on scales comparable to or greater than the diffraction limit. On the other hand, the NRMSEs of the reconstructed images start to deviate from the model images in the super-resolution regime—namely on scales smaller than the diffraction limit. In this regime, the NRMSEs differ by technique. CS-CLEAN has a common trend for all four models, which is broadly consistent with results of Chael et al. (2016). They achieve minimum errors at a resolution of $\sim 30-60$% of the diffraction limit and then show a rapid increase in errors at smaller scales. In contrast, closure-phase imaging with $\ell_{1}$+TV regularizations show much more modest variations in the super-resolution regime. They achieve minimum errors at a resolution of $\sim 20-30$% of the diffraction limit, smaller than CS-CLEAN, and show only a slight increase at smaller scales. $\ell_{1}$+TV reconstructions produce images that have a smooth distribution similar to the model images, resulting in smaller errors than CS-CLEAN, even if the $\ell_{1}$+TV reconstructions are not convolved with a restoring beam.

In super-resolution regimes, the errors in $\ell_{1}$+TV images mostly arise from the presence of tiny substructures in the image. For instance, we show model, reconstructed and residual images filtered with different baseline lengths for DJ1 in Figure 4. As shown in Figure 4, residuals are small when filtering baseline lengths are shorter than the maximum baseline length. On the other hand, for longer filtering baseline lengths, systematic residuals due to tiny substructures much smaller than the diffraction limit start to appear, which can be traced only with baselines longer than the simulated observations.

Although NRMSEs are better than CS-CLEAN at finer resolutions, $\ell_{1}$+TV images have broader emission region sizes regardless of models. This is due to a typical feature of images reconstructed with the isotropic TV, which prefers flat images with sharp edges. Since the simulated data do not have visibilities at baseline lengths long enough to resolve the width of ring- or crescent-like features, TV enlarges their widths until images start to deviate from observed visibilities. This property of the isotropic TV regularization would be useful to constrain the upper-limit size of the emission regions and black-hole shadow (see §5.1 for further discussions). On the other hand, in terms of the image fidelity, these results suggest that regularizations preferring much smoother edges are preferable; smoother gradients in the image lead to images with higher contrast (i.e. brighter/fainter pixels become even brighter/fainter, respectively) to conserve the total flux, which often makes the effective emission region size smaller. In §5.2, we discuss alternative regularizations of sparse imaging reconstruction for improving the image fidelity in super-resolution regime.

4.2 Regularization parameters and Cross Validation

Reconstructed images for all 16 sets of the regularization parameters are shown in Figure 5 for the accretion-flow-dominated model. As shown in Figure 5, the reconstructed images with $\tilde{\eta}_{\rm l},\tilde{\eta}_{\rm t}\lesssim 1$ have noisy artifacts. This is true for all four models. Such artifacts may appear because the images are poorly constrained by both regularization functions at lower regularization parameters.

CV works as an Occam’s razor and prevents over-fitting. In Figure 6 (a), we show the residual $\chi^{2}$ between the validation set and the model image reconstructed from the training set, which is averaged for all 10 trials (henceforth the CV value). As expected in §2.2, the CV value tends to be large for large regularization parameters, since the regularization functions too strongly constrain the image so that the model image is inconsistent with observational data. Once the CV value achieves its minimum value (at a parameter set marked with stars in Figure 6 (a)), it starts to increase again for lower regularization parameters, since the model image is over-fitted to the training set and then shows larger deviations from the validation set.

Since the groundtruth images are known in this work, one can compare CV-selected images with images on other parameter sets, which can not be done for real observational data with unknown ground truth images. Although CV selects the optimal parameter based on the noise level of the data, it does not guarantee that the selected parameter achieves the best imaging fidelity among all parameter sets examined. In Fig 6 (b), we show the NRMSE of non-Gaussian-convolved images for all parameter sets. Although metrics for the image fidelity of best-fit images are slightly larger than the minimum values for the three of four models (J2, M16 and DJ1), they are consistent to within a few percents. These slight differences between the best-fit and best-fidelity parameters would not produce substantial differences in the image, and the resulting images are good enough.

The above results clearly demonstrate that CV is a useful technique to determine the regularization parameters so that the reconstructed image does not overfit noises in the data. We emphasize that CV is a general technique and can be applied to imaging regardless of the specific data products used (e.g. full complex visibilities, visibility amplitudes and/or closure quantities) or chosen regularization (for instance, sparse modeling; MEM: e.g., Buscher 1994, Chael et al. 2016; patch priors: e.g., Bouman et al. 2015; scattering optics: Johnson 2016).

5 Discussions

5.1 Implications for future EHT observations of M87

The proposed method successfully reproduces a clear signature for the black hole shadow for the models J2, M16 and DJ1. This is as expected from the visibility distribution of these models, which have null amplitude regions, created by the shadow feature, at intermediate baseline lengths of $\sim 3-4$ G$\lambda$ (Akiyama et al., 2015). MEM also succeeds in reproducing the black hole shadow in these models (Lu et al., 2014). Presence of a clear shadow feature is tightly connected to where the dominant emission originates, since the silhouette of the black hole is created by the photons produced a few $R_{s}$ from the black hole, illuminating the last photon orbit (see discussion in Akiyama et al., 2015). As demonstrated in previous imaging simulations (Lu et al., 2014), future EHT observations can constrain the loading radius of the high-energy leptons producing synchrotron emission at 1.3 mm via the appearance of the black hole.

The most important implication of this work is that our regularization function TV and parameter selection with CV will enlarge the ring- or crescent-like emission illuminating the last-photon orbit of the black hole as much as possible, within the range that the model image neither over-fits nor deviates too much from the observed visibilities. Hence, the obtained width of the surrounding emission in the reconstructed image is close to an upper-limit on the width of the emission region, and simultaneously, the obtained diameter of the black hole shadow should be interpreted as a reasonable lower limit for it. The clear shadow features in the reconstructed images for models J2, M16 and DJ1, therefore, strongly indicate that the EHT can sample a large enough range of visibilities with appropriately low noise levels to image the black hole shadow. In addition, the raw reconstructed image with the proposed method can be used to constrain the lower-limit size of the shadow. This would be useful to constrain the mass of M87, which has an uncertainty of an factor of about two between the stellar-dynamical (e.g. Gebhardt et al., 2011) and gas-dynamical (e.g. Walsh et al., 2013) modeling, and therefore be informative to clarify that which of modeling methods is desirable to measure the mass of super-massive black holes.

Another important implication is that, at least for the black hole images, post-processing Gaussian convolution would not be required with the $\ell_{1}+TV$ regularization, although the CLEAN techniques do require it to reduce many compact artifacts in the image. As shown in Figure 3, the NRMSE curves for the $\ell_{1}+TV$ regularization are shallow for smaller convolving sizes, and applying a circular Gaussian beam therefore makes only small improvements of a few percent in the NRMSE regardless of the input model images. Similar results are also shown in recent work with the MEM (Chael et al., 2016). Our results support that application of the beam is not required for the recent state-of-art imaging methods utilizing multi-resolution regularization functions for imaging the $R_{s}$-scale structure of M87 and Sgr A* with the EHT.

5.2 Relevant future issues

5.2.1 Other sparse regularization for smoothed images

In this paper, we adopt TV regularization, which favors a smooth image, so that images can be reconstructed with smaller pixel sizes and/or for more extended sources. As shown in $\S$4, the reconstructed images have good image fidelity. In particular, it is noteworthy that CV selected high values of $\tilde{\eta}_{t}=10^{0}-10^{1}$ for all models (see Figure 6 (a)), suggesting that the solutions would be over-fitted without TV regularization. These results demonstrate that inclusion of the TV regularization can extend the range of objects where sparse modeling is applicable.

However, as described in §4.1, the reconstructed images have larger emission regions than the original model, since the isotropic TV (Rudin et al., 1992) adopted in this work prefers a flat brightness distribution in super-resolution regimes where the image cannot be constrained well. Hence, a regularization preferring more smoothed edges is required to improve the fidelity for the black-hole imaging with the EHT.

In the context of sparse reconstruction, there are several candidates for improving the image fidelity as a natural extension of this work. First of all, there are other forms of regularization functions, which prefer sparse images in the gradient domain with smoother edges. For instance, an alternative form, given by,

[TABLE]

is also convex like the isotropic TV term adopted in this work, and prefers images with smoother edges. Furthermore, previous studies of sparse image reconstruction techniques have shown that regularization with $\ell_{1}$+wavelet/curvelet transformation is also a promising approach (e.g. Li et al., 2011; Carrillo et al., 2012, 2014; Garsden et al., 2015; Dabbech et al., 2015). We will test these sparse regularizations in a forthcoming paper.

5.2.2 Enhancing dynamic range with self-calibration

In VLBI, the visibility phase is initially calibrated with fringe fitting (also called as fringe search), which is a self-calibration technique using phase closure (see Thompson et al., 2001). The fringe fitting can mitigate most station-based errors due to atmospheric and instrumental effects, although errors may remain if an incorrect source model is assumed. Traditionally, self-calibration with hybrid/differential mapping (e.g. Walker, 1995) has been employed to solve for residual structural phase errors and images simultaneously, which has been successful for VLBI imaging.

This work and previous works on closure-phase imaging techniques using other regularizations such as MEMs (Lu et al., 2014, 2016; Fish et al., 2014; Chael et al., 2016) and patch priors (CHIRP; Bouman et al., 2015) demonstrate that an image can be reconstructed with high fidelity even from closure quantities. However, since closure phases and other closure quantities have less information about the source structure and also larger thermal noises than full complex visibilities, imaging with closure quantities can limit the dynamic range, image sensitivity and optimal spatial resolution.

A promising approach to improve the dynamic range is to use a reconstructed image from closure imaging techniques as an initial image for hybrid/differential mapping. It can also be used as a model for fringe fitting, as self-calibrated images are often applied to detect more fringes on faint sources (e.g. Hada et al., 2016). In a forthcoming paper, we will evaluate the performance of such a hybrid-mapping technique including closure-phase/full-closure imaging priors. We emphasize that, for this purpose, one does not need to reconstruct the image with pixels much smaller than scales where the brightness distribution cannot be constrained by data and therefore will not affect results of self-calibration.

6 Conclusions

We have presented a new imaging technique reconstructing images from visibility amplitudes and closure phases by utilizing two regularizations of sparse modeling: the $\ell_{1}$-norm and Total Variation (TV). Furthermore, we also propose a method to select optimal regularization parameters with cross validation (CV), which can be applied to most existing imaging algorithms. As an example, we applied our technique to simulated observations of M87 with the EHT at 1.3 mm. Here, we summarize our conclusions.

We find that $\ell_{1}$+TV regularization can achieve an optimal resolution of $\sim 20-30$% of the diffraction limit $\lambda/D_{\rm max}$, which is the nominal spatial resolution of a radio interferometer. This optimal resolution is better than that of the most-widely used Cotton-Schwab CLEAN, which uses full complex visibilities. 2. 2.

We confirm that cross validation (CV) works as an Occam’s razor and prevents over-fitting when selecting the optimal regularization parameters. CV is a general method that can be applied to interferometric imaging more generally, such as imaging with full-complex visibilities and/or using other regularizations. 3. 3.

Using $\ell_{1}+$TV regularization, the reconstructed image maximizes the width of the emission region within the range that it neither over-fits nor deviates too strongly from the data. Hence, the clear reproduction of the black hole shadow in the reconstructed image suggests that future EHT observations will have the $uv$-coverage and sensitivity sufficient for imaging it. In addition, the reconstructed image will be able to constrain the sizes of the black hole shadow and surrounding emission region.

Finally, we remark that all of above results demonstrate the clear promise of the EHT for providing an unprecedented view of the event-horizon-scale structure of the super-massive black hole in M87 and also the Galactic center Sgr A*.

We thank the anonymous referee for his/her useful and constructive suggestions to improve the paper. K.A thanks members Dr. Michael D. Johnson and Dr. Lindy Blackburn for many useful suggestions on this work. K.A. and this work are financially supported by the program of Postdoctoral Fellowships for Research Abroad at the Japan Society for the Promotion of Science. M.M. acknowledges support from the ERC Synergy Grant (Grant 610058). Event Horizon Telescope work at MIT Haystack Observatory and the Harvard Smithsonian Center for Astrophysics is supported by grants from the National Science Foundation (NSF; AST-1440254, AST-1614868) and through an award from the Gordon and Betty Moore Foundation (GMBF-3561). Work on sparse modeling and Event Horizon Telescope at the Mizusawa VLBI Observatory is financially supported by the MEXT/JSPS KAKENHI Grant Numbers 24540242, 25120007 and 25120008.