Time Domain Filtering of Resolved Images of Sgr A^∗

We evolve the accretion flow using the conservative GRMHD code harm3d (Noble et al. 2006, 2009). The simulation starts from a constant angular momentum equilibrium torus (Fishbone & Moncrief 1976) that is perturbed by a weak magnetic field to seed the growth of the magnetorotational instability (MRI). The disk is integrated in Kerr–Schild spacetime in modified spherical-polar coordinates, with a logarithmically scaled radial grid. The computational domain runs from within the event horizon to $240{{GM}}_{\mathrm{BH}}/{c}^{2}$ in radius and extends over a full $2\pi$ radians in azimuth. The radial and poloidal boundary condition are set to outflow and reflective boundary, respectively. The simulation is terminated after ${\rm{14,000}}{{GM}}_{\mathrm{BH}}/{c}^{3}$, which is long compared to the timescale for the MRI to reach saturation but short compared to the accretion timescale for the torus. Our initial conditions have only a small net vertical magnetic flux through the disk, and so produce a Standard and Normal Evolution (SANE, Narayan et al. 2012) disk model with short timescale variability that is qualitatively consistent with observational results (Chan et al. 2015a).

We run two models with different black-hole spin: a rapidly spinning black hole with a = 0.9375 (hereafter 0.94) with resolution $260\times 192\times 128$ (in radius, colatitude, and longitude respectively), and a nonrotating black hole (a = 0) with resolution $144\times 144\times 144$. The difference in resolution arises because of computational limitations. One might worry that the difference would produce differences in time variability properties, but this does not appear to be the case; as one indication, the azimuthal correlation length of the emissivity, which is tightly correlated to the observational variability as will be discussed in Section 4, differs by only 10% on average at $r\lt 12{{GM}}_{\mathrm{BH}}/{c}^{2}$ (azimuthal correlation length of fluid variables is discussed in Shiokawa et al. 2012), while the Innermost Stable Circular Orbit (ISCO) orbital frequency differs by a factor of ∼4; 34 minutes and 9 minutes for a = 0 (${R}_{\mathrm{ISCO}}=6{{GM}}_{\mathrm{BH}}/{c}^{2}$) and 0.94 (${R}_{\mathrm{ISCO}}\sim 2.04{{GM}}_{\mathrm{BH}}/{c}^{2}$), respectively. This suggests that the artificial error introduced by the resolution difference has only a minor effect on the observational variability in our radiative models.

Hereafter, we describe both the length unit ${{GM}}_{\mathrm{BH}}/{c}^{2}\,\sim 6.64\times {10}^{11}$ cm and time unit ${{GM}}_{\mathrm{BH}}/{c}^{3}\sim 22\,{\rm{s}}$ as M by setting ${{GM}}_{\mathrm{BH}}=c=1$ unless otherwise indicated.

The next step is to construct time-variable images of the model. The total duration of our model is $2500\,M$ (∼15 hr for Sgr A*) taken from the last part of our accretion disk simulation. We use the ray-tracing method of Noble et al. (2007, 2009) that integrates the transfer equation along null geodesics, incorporating synchrotron emission and absorption. This leads to a grid of intensities (pixels) on a "camera" at Earth's distance from Sgr A*. We produce images at $\lambda =1.3\,\mathrm{mm}$, where EHT will observe.

In modeling Sgr A*, it is common to use a "fast light" approximation, in which the emergent radiation is calculated on individual time slices and changes in the fluid on the light crossing time are ignored. But because we are interested in intensity fluctuations on timescales comparable to the light crossing time, and plasma close to the event horizon is moving at relativistic speed, it is necessary to to use a full "slow light" treatment in which the geodesics are integrated through an evolving flow. To do this, we update the background snapshot of the disk every 0.5 M ∼ 11 s, as measured by a distant observer. A photon is advanced by dxⁱ along its geodesics at each time step ${{dx}}^{0}={dt}=0.5\,M$, following the relation

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}^{i}}{{dt}}=\displaystyle \frac{{{dx}}^{i}}{d\lambda }/\displaystyle \frac{{dt}}{d\lambda }=\displaystyle \frac{{k}^{i}}{{k}^{0}}.\end{eqnarray} \tag{ 1 }$$

Here, λ is the affine parameter (not the wavelength; the difference should be clear from the context) and ${k}^{\mu }$ is the wave four vector.

The model parameters are the mass of the black hole, the distance to Sgr A*, the ratio of proton to electron temperature $R\equiv {T}_{{\rm{p}}}/{T}_{{\rm{e}}}$, the accretion rate $\dot{M}$, and the inclination angle. The mass and distance to Sgr A* are set to $4.5\times {10}^{6}\,{M}_{\odot }$ and 8.4 kpc, respectively (Ghez et al. 2008; see also Boehle et al. 2016 for the recent estimate). There are multiple combinations of R and $\dot{M}$ that are broadly consistent with the data (see, e.g., Mościbrodzka et al. 2009; Broderick et al. 2011; Mościbrodzka et al. 2014, for parameter-fitting exercises). Here we simply fix R = 3 and adjust $\dot{M}$ so that the time-averaged 1.3 mm flux matches the observed 3.6 Jy (Bower et al. 2015). Finally, for the inclination, we explore two viewing angles: a face-on ($i=2^\circ$) and edge-on ($i=90^\circ$) view of the disk. Since the main purpose of this paper is to introduce the time-domain filtering technique, we only examine four models, i.e., edge-on and face-on views for spin = 0 and 0.94, to convey a rough sense of how results depend on model parameters.

Figure 1 shows ray-traced images of our disk models (a = 0 and 0.94) for face-on and edge-on views of one time slice. The field of view (FOV) and image resolution are $30\times 30M=160\,\times 160\,\mu \mathrm{as}$ and 128 × 128, respectively, which gives a pixel size of $1.25\times 1.25\,\mu \mathrm{as}$. In the edge-on image, the left side of the images is brighter due to Doppler beaming as the plasma moves toward the camera on the approaching side of the disk. The innermost part of the bright arcs in the edge-on images and the bright rings of the apparent radius ∼25–30 μas in the face-on images are the so-called "photon ring." They are emission from slightly outside the photon orbit in the equatorial plane, where ${R}_{\mathrm{ph}}=3M(1.43M)$ for $a=0(0.94)$ respectively. The apparent size of the photon orbit is increased by gravitational lensing (Bardeen 1973).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is convenient to define a characteristic location of emission along each ray, since our hypothesis is that variability along a ray is tightly correlated with its point of emission. Suppose ${\rm{\lambda }}$ is the affine parameter along a geodesic. Then,

$$\begin{eqnarray}&&\bar{\lambda }\equiv \displaystyle \frac{\int \lambda {\left(\tfrac{{dI}}{d\lambda }\right)}^{2}\,d\lambda }{\int {\left(\tfrac{{dI}}{d\lambda }\right)}^{2}\,d\lambda },\end{eqnarray} \tag{ 2 }$$

where $I(\lambda )$ is specific intensity. We define the point at $\bar{\lambda }$ along the ray as the "emission point," ${\boldsymbol{x}}(\bar{\lambda })={{\boldsymbol{x}}}_{\mathrm{emiss}};$ $r(\bar{\lambda })\equiv {r}_{\mathrm{emiss}}$ is the characteristic "emission radius" and $\theta (\bar{\lambda })\equiv {\theta }_{\mathrm{emiss}}$ is the characteristic emission latitude. The bright regions in our model images have ${\theta }_{\mathrm{emiss}}\sim \pi /2$, i.e., the disk's equatorial plane.

Figure 2 shows r_emiss for the face-on and edge-on views of both spins. For the edge-on view, the inner edge of the photon ring has the smallest r_emiss (the blue rings in Figure 2). Most of the emission of the rings is from the side and behind the black hole, while the region inside the rings is dominated by emission from gas in front of the black hole; these geodesics extend into the horizon rather than going around the hole. Although the shadow angular radius is insensitive to spin, the range of r_emiss that lights up the photon ring is strongly dependent on spin: the minimum r_emiss in the equatorial plane for the a = 0 model is $\sim 4.5M$ and becomes $\sim 8M$ as the image position moves $20\,\mu \mathrm{as}$ outward, while the same measurement gives $2M$ and $6M$ for the a = 0.94 model. This introduces a difference in image variability for the two spin models, as discussed in the next subsection.

Zoom In Zoom Out Reset image size

For the face-on view, emission radius monotonically increases outward from the center of the images except for the region inside the white dotted ring. There, a very low emission from extended regions outside the disk (jet) dominates. There is a noticeable difference between the different spin models inside the photon rings, indicated by the solid white rings, where r_emiss is smaller by ∼1.5–2M for the a = 0.94 model for a given radius from the image center. The use of our r_emiss definition in the photon ring is inadequate for the face-on view since light rays penetrate the disk multiple times at different radii as they go around the black hole.

In practice, the total intensity in a single pixel is built up from an extended region around its emission point. If this emission region is small, the time variability of the intensity reflects dynamics at the emission point well. If the emission region is large, the variability is built up from variations at multiple locations and it is difficult to extract information from the light curve. Nevertheless, a comparison of Figure 2 with time-domain filtered images that are presented in Section 3.2 (Figures 4 and 5) shows a strong correlation between the emission radius and time variability in our model images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We begin by examing the power spectrum density (PSD) of the total (source integrated) mm flux. The PSD is constructed as follows. For a light curve L_n, where $0\leqslant n\leqslant N-1$ is the index of the measured point in time and N is the total number of points, its discrete Fourier transform is

$$\begin{eqnarray}&&\tilde{{L}_{s}}=\displaystyle \sum _{n=0}^{N-1}{w}_{n}{L}_{n}{e}^{2\pi {isn}/N},\end{eqnarray} \tag{ 3 }$$

where $-N/2\leqslant s\leqslant N/2$ is the index for the temporal frequency and w_n is a Hamming window used to reduce the edge artifacts. Then the PSD is defined as

$$\begin{eqnarray}&&{P}_{s}\equiv \tilde{{L}_{s}}\tilde{{L}_{s}^{* }}/W,\end{eqnarray} \tag{ 4 }$$

where $W=N{\sum }_{n=0}^{N-1}{w}_{n}^{2}$ (Press et al. 2002). We divide each light curve into four segments and average the resultant PSDs, where $N\sim 1200$ per segment in our operation. The duration of each segment is ∼3.5 hr (time resolution of our data set is 0.5M ∼ 11 s), which is close to the integration time of the realistic EHT observation. Averaging over segments also improves the signal-to-noise ratio.

Figure 3 shows the power spectral density (PSD) of the total flux variability for our four models at $a=0,a=0.94$ and $i=0^\circ ,i=90^\circ$. The vertical lines indicate the ISCO frequency for a = 0 (broken) and a = 0.94 (solid). The total flux in our models comes from an extended range in radius, so its variability is not sensitive to disk dynamics at any particular radius. Overall, edge-on views have higher power in a wide band around the ISCO frequency than face-on views (compare black and red lines of the same line type), while no systematic difference is apparent between different spins with the same viewing angle (compare the same colors).

All the models exhibit a monotonic decrease in power with frequency. There is weak evidence for a break close to the ISCO frequency, with $d\mathrm{ln}P/d\mathrm{ln}f\simeq -2$ below and $d\mathrm{ln}P/d\mathrm{ln}f\gtrsim -3$ above. Observations (Dexter et al. 2014) show an ∼f⁻² spectrum between ${f}_{\min }\,\simeq \,1.4\times {10}^{-3}\,\mathrm{Hz}$ and ${f}_{\max }\,\simeq \,3\times {10}^{-3}\,\mathrm{Hz}$, with a turnover to a flat spectrum at frequencies below $\sim 3.5\times {10}^{-5}\,\mathrm{Hz}$. In a future investigation, we will consider the long-timescale behavior and the observed break at $\sim 3.5\times {10}^{-5}\,\mathrm{Hz}$ by extending the duration of our simulation. The peak radius of emission in the edge-on views is approximately 5–6M, but there is still a large contribution to the total flux from smaller radii. The peak emission for face-on views comes from the photon ring (see the left column in Figure 1), but also from a large range in radius. These highly extended origin of emission for the face-on models smear out the variability amplitude and lead to an overall lower power than in edge-on models. Evidently, spatially resolved structural variability is essential if we want to extract more detailed information about the disk dynamics. However, Figure 3 implies the total flux PSD may allow discrimination between different spin models.

Next, we consider the effects of filtering the mock images in the time domain. We first obtain the light curve of each pixel in Figure 1 and produce a corresponding PSD following the procedure described in Section 3.1. The PSD is integrated over a frequency band, and the resulting power is assigned to the same pixel to generate a time-domain filtered image. Images reconstructed from observational data will, of course, have coarser angular resolution. After obtaining the sets of PSD of the pixels in the ray-traced images, we analyze their spatial distribution by mapping the power in a band or frequency range of ${f}_{1}\lt f\lt {f}_{2};$ ${ \mathcal P }(i,j,{f}_{1},{f}_{2})={\int }_{{f}_{1}}^{{f}_{2}}P(i,j,f){df}$, where i and j are the spatial indices of a pixel.

Figures 4 and 5 show unnormalized time-domain filtered images ${ \mathcal P }$ for $({f}_{1},{f}_{2})=({10}^{-3.8},{10}^{-3.3})$, $({10}^{-2.8},{10}^{-2.3})$, and $({10}^{-1.8},{10}^{-1.3})$ Hz (notice that ${f}_{\mathrm{ISCO}}={10}^{-2.7}$ and ${10}^{-3.3}$ Hz for a = 0.94 and 0, respectively) together with the time-averaged intensity map (f = 0) for the face-on and edge-on views in a logarithmic color scale. The figures are over-plotted with ellipses to represent the image size and orientation. To obtain these ellipses, we first construct the covariance matrix (second moments) of the image pixels and find its eigenvectors and eigenvalues, ${\lambda }_{1}$ and ${\lambda }_{2}$. The eigenvectors correspond to the principal axes of the images and their length is the square root of the eigenvalues. The length of the vectors is related to the FWHM of the axisymmetric Gaussian model used for the analysis of VLBI observations by FWHM/2.3. The length of the principal axes of the ellipses in Figures 4 and 5 is the FWHM, i.e., ${\mathrm{FWHM}}_{\mathrm{1,2}}=2.3\sqrt{{\lambda }_{\mathrm{1,2}}}$. The eccentricity of the ellipses $\equiv \,\sqrt{1-{\lambda }_{2}/{\lambda }_{1}}$.

The image size decreases as the filtering frequency increases for all the models in Figures 4 and 5. In general, the smaller r_emiss, the higher the power at high filtering frequencies and the lower the power at low filtering frequencies if the intensity is fixed. This naturally arises because the high orbital frequency at the inner radii of the disk produces rapid variability. Emission from large radii has relatively higher power at low filtering frequencies, and hence the image size is larger due to the extended area of emission from outer radii (see Figure 2). On the other hand, emission from small radii occupies a smaller region in the images and the image size gets smaller for high filtering frequencies. This result confirms that our filtering technique is able to visualize a spatial distribution of the local characteristic timescale at r_emiss. Figure 6 shows the relation between the image size and the filtering frequency range for all the cases, where the "size" is the average of FWHM₁ and FWHM₂. The size decreases with increased filter frequency, and ranges from 30% to 50%, depending on the model considered.

Zoom In Zoom Out Reset image size

The high-frequency filtered edge-on images (rightmost column in Figure 4) nicely traces the distribution of small emission radii (see Figure 2) for both of the high-spin models. This suggests that image eccentricity has a dependence on black-hole spin, especially at high filtering frequencies. The eccentricity of the edge-on filtered images, in general, increases with the filtering frequency, and the high-spin model has higher eccentricity at all the filtering frequencies than the low-spin model; see Figure 7 for the dependence of eccentricity on the filtering frequency. Unlike the edge-on views, eccentricity in the face-on images does not depend much on the filtering frequency. In principle, the degeneracy of spin and viewing angle can be avoided by the unique dependence of eccentricity on the filtering frequency. Surprisingly, face-on images are not always very "round," as we see the face-on a = 0 model in Figure 7 has eccentricity $\gt 0.4$ due to the presence of slowly evolving non-axisymmetric structures.

Zoom In Zoom Out Reset image size

The time-domain filtering technique highlights the direction of the disk's orbital axis, as the semimajor axes of the ellipses are aligned with it. Even if the black-hole shadow is not well resolved or apparent in time-averaged data, filtering at high frequency would show anti-symmetric properties in the variable structure as long as the disk is not viewed face-on. We note that our technique might be used to test if an observed time-averaged feature is a true black-hole shadow: if it is, then a similar feature should grow more prominent as the filtering frequency increases.

Real EHT data is sparsely sampled in the spatial Fourier domain (U–V plane) and has superposed atmospheric and thermal noise as well as interstellar scattering. Thus reconstruction of even static images require sophisticated approaches (e.g., Bouman et al. 2016; Chael et al. 2016). Nevertheless, reconstruction of high time resolution images from incomplete UV-coverage is possible using several newly developed techniques (Johnson et al. 2017, ApJ, submitted).

We performed a preliminary test to apply the time-domain filtering technique to time-variable images produced by a synthetic observation of our high-spin edge-on view model, using the array configuration for the EHT 2017 campaign. However, the effects of interstellar scattering and realistic, irregular time sampling (i.e., due to telescopes off-source f or phase calibration) were not included. We found that the technique could reproduce the trend by which the image size increases and eccentricity decreases as the filtering frequency increases.

Here, we investigate the effect of using a "fast light" approximation by comparing fast light and slow light time-domain filtered images. Implementation of "slow light" in the General Relativistic Ray-tracing code has been done by several researchers (e.g., Dolence et al. 2012, J. Dexter 2016, private communication) and it is known to have non-negligible effects on high-frequency variability (C.-k. Chan & T Bronzwaer 2017, private communication).

We found that the approximation introduces non-negligible effects on filtered image morphology for the edge-on view, at the high filter frequencies. As an example, Figure 8 shows a comparison of the time filtered images with and without the fast light approximation, where the properties of the images are the a = 0.94 model, edge-on view, and filter frequency range of ${10}^{-2.3}\mbox{--}{10}^{-1.8}\,\mathrm{Hz}$. The power around the equatorial plane in the fast light image is approximately 50% less than in the corresponding slow light image. Evidently, the use of slow light is essential for our time-domain filtering technique.

Zoom In Zoom Out Reset image size

For the same model, the effect of the approximation on the time-domain filtered image is much less apparent at the low filter frequencies, approximately below ISCO orbital frequency $\sim 1.9\times {10}^{-3}={10}^{-2.7}\,\mathrm{Hz}$. The total flux light curve of the slow light calculation is slightly smoother than that of the fast light but no qualitative difference is apparent.