Millimeter Light Curves of Sagittarius A* Observed during the 2017 Event Horizon Telescope Campaign

The Sgr A* image field of view can be split into two components at angular scales of arcseconds probed by ALMA:

• 1.  
  An extended structure with a low brightness temperature, primarily originating from thermal emission from ionized gas and dust infalling into the central region of the Galaxy, the so-called "minispiral" (e.g., Lo & Claussen 1983; Goddi et al. 2021). From our ALMA observations, the integrated extended flux density of the minispiral is ∼1.1 Jy. Given the physical origin of this emission and the spatial scales involved (several tens of parsecs), we can assume the brightness distribution of the minispiral to remain constant during the few days of the EHT observing campaign.
• 2.  
  An unresolved and highly variable component corresponding to the compact source Sgr A*, with a flux density typically ranging between 2 and 5 Jy at 230 GHz.

With the superb sensitivity of ALMA (M87* Paper III; Paper II) and the sufficiently high integrated extended flux density of the minispiral, it is possible to detect the whole field-of-view structure (i.e., the minispiral plus Sgr A*) in each ALMA 4 s snapshot. Therefore, one can assume a two-component Fourier domain model, ${V}_{t}^{\mathrm{mod}}$, composed of (1) ${{ \mathcal F }}^{e}$, a Fourier transform of the static extended minispiral, F^e , corrupted with a time-dependent amplitude gain, G_t , accounting for atmospheric and instrumental effects (following the QA2 calibration, these effects can be modeled with a single function of time, representing the effective gain of the interferometric array); and (2) F_t , an unresolved Sgr A* compact component with a time-dependent flux density (still corrupted by the G_t gain at this stage), so

$$\begin{eqnarray}&&{V}_{t}^{\mathrm{mod}}={G}_{t}{{ \mathcal F }}^{e}+{F}_{t}.\end{eqnarray} \tag{ 1 }$$

If we denote the visibility observed at a time t on a baseline i as ${V}_{i,t}^{\mathrm{obs}}$ and the model sampled at the same Fourier plane location as ${V}_{i,t}^{\mathrm{mod}}$, the model can be fitted to the data by minimizing

$$\begin{eqnarray}&&{\chi }_{t}^{2}({G}_{t},{F}_{t})=\sum _{i}{\omega }_{i,t}{\left|{V}_{i,t}^{\mathrm{obs}}-{V}_{i,t}^{\mathrm{mod}}\right|}^{2}\ \end{eqnarray} \tag{ 2 }$$

for each time t, with S/N-based baseline weights, ω_i,t , and the summation extending over all baselines available at a given time t (Martí-Vidal et al. 2014).

Since the true integrated flux density of the minispiral is assumed to be constant, we can use the values of G_t to remove the residual corruption effects in the Sgr A* flux density estimates, F_t . Hence, we produce a corrected estimate of the Sgr A* flux density, F^c _t , using the equation

$$\begin{eqnarray}&&{F}_{t}^{c}=\displaystyle \frac{{F}_{t}}{{G}_{t}}.\end{eqnarray} \tag{ 3 }$$

In practice, we also need to solve for the image domain minispiral model, F^e . We use the CLEAN algorithm (e.g., Högbom 1974) implemented in the Common Astronomy Software Application (CASA) framework (McMullin et al. 2007) as the tclean task, iteratively reconstructing the image of the minispiral, recalibrating the data with G_t , and updating G_t and F_t . While the minispiral structure is assumed to be constant across observed frequencies, the absolute flux density scale is allowed to vary between the subbands. The procedure runs until convergence. The times with unphysical or unconverged (G_t , F_t ) are flagged. Special attention is given to the minispiral total flux density, which is fixed per subband to the median of the flux densities estimated from all of the snapshots obtained throughout the EHT campaign.

In Figure 1, we show two model images of the field around Sgr A*; the left panel corresponds to the original QA2 calibrated data (the initial condition for the iterative procedure) and shows the complete source structure (minispiral and Sgr A*); the right panel is the final minispiral model obtained after the convergence of the intrafield calibration. The high noise level seen in the QA2 image (i.e., the artifacts that are distributed across the whole field of view) is due to the effects of the time variability of Sgr A*. By modeling the source variability with Equation (1), the noise level in the final minispiral image (Figure 1, right) is reduced. We notice that the size of the ALMA primary beam at 230 GHz is of the order of the size of the minispiral structure. This implies that different points across the image are affected by a different ALMA sensitivity. Only emission from regions with a primary-beam response >5% can be considered as detections above 5σ of the ALMA sensitivity. As we approach this primary-beam threshold, marked with a dashed line in Figure 1, the noise effects in the brightness distribution increase. Regions far from the phase center have a negligible contribution to the visibilities, since such a contribution also scales with the primary-beam response. The effects of these regions on the calibration of the Sgr A* light curves will thus be small.

Zoom In Zoom Out Reset image size

In Figure 2, we show the visibility amplitudes for a representative ALMA snapshot from 2017 April 7, band B1. The contribution of the extended minispiral model (green crosses) at short baselines is clearly visible. As a final product, we obtain converged light curves of Sgr A* with a snapshot cadence of 4 s. Data corresponding to a source elevation below 25°, exhibiting significant quality loss, are flagged.

Zoom In Zoom Out Reset image size

The time-dependent factor, G_t , can also be used to correct the amplitudes of the VLBI visibilities related to the phased ALMA. Based on Martí-Vidal et al. (2016), the scaling gain factor to correct ALMA amplitudes on VLBI baselines is $\sqrt{{G}_{t}}$. This approach has been employed for a priori amplitude calibration of the EHT VLBI Sgr A* data (see Appendix A and Paper II).

A custom script was prepared to process the nonstandard array data acquired during the phased ALMA observations of Sgr A* using measurements of the system equivalent flux density (SEFD) of each ALMA antenna. While similar to the standard ALMA QA2 pipeline, it includes additional calibration steps necessary to produce the time-dependent light-curve data. The ALMA observations are grouped in scans, which consist of subscans of 18 s cadence, with 16 s of correlated data (Goddi et al. 2019). A nearby phase calibrator, J1744–3116 (J1744–312), was observed for 30 s every 20 minutes. Observations of two bright quasars, NRAO 530 (J1733–1304, B1730–130) and J1924–2914 (B1921–293), were also included for the amplitude calibration; see Table 2. First, the phase delays associated with the atmospheric water vapor were estimated from measurements of the 183 GHz water line, performed with high time cadence using radiometers located at each ALMA antenna. The radiometer measurements allowed us to estimate the column of water vapor above each ALMA antenna, which were then converted into a phase correction related to the atmospheric optical path. Conversion from the relative visibility correlation amplitude to a flux density scale was performed by applying the system temperature measurements performed routinely at each antenna. The corrected data were concatenated and reduced to produce a single CASA measurement set (McMullin et al. 2007) for each observing day, containing relevant data for all four subbands.

The second step was the bandpass calibration of all of the frequency channels in each subband. We used NRAO 530 to generate the bandpass calibration tables, choosing a scan when the source was nearest to the zenith. The chosen reference antenna was located near the array center and not shadowed by neighboring antennas. Lower-sensitivity channels near the edges of each subband, as well as data from shadowed antennas, were flagged. Following bandpass calibration, all channels within each subband were averaged.

In the third step, we determined the amplitude scale of the observations on each day, and we applied the phase-referencing calibration. Since all three calibrators were observed as unresolved point sources, anomalously low amplitudes on some antennas were apparent. The few low-amplitude data points, with less than about 70% of the nominal sensitivity of the majority of antennas, were subsequently flagged. Next, the flux density scale over the entire observation was set using NRAO 530 as a flux density calibrator, assuming a flux density of 1.56 Jy at 213.1 GHz and a spectral index of −0.72, obtained from the ALMA calibrator catalog. ¹⁵⁰ The flux density calibration was applied using the average gain of the two polarizers of each antenna (X and Y), so that there was no gain bias caused by the source linear polarization.

The flux densities of the other two quasar calibrators were verified to be constant over each observation day to within a few percent. The last calibration step, phase referencing between J1744–3116 and Sgr A*, was performed by deriving the antenna-based phase for each J1744–3116 scan with the CASA task gaincal and then interpolated to each Sgr A* scan, which completed the calibration cycle.

The above steps provide calibrated complex visibilities from which an image containing a strong point source, Sgr A*, and the extended minispiral emission can be obtained. Since the phase calibrator J1744–3116 was observed only every 20 minutes, the interpolated phase correction can deviate from the true values. To determine the time variability of Sgr A* in the presence of extended emission and calibration phase errors, we first flag baselines shorter than 70 m (about 50 kλ; Figure 2). Since virtually all of the extended emission is resolved out on longer baselines, the remaining data are reasonably consistent with a point-source model, although its position may vary in time. Subsequently, given a large number of available long baselines, we perform phase self-calibration to remove the residual phase errors remaining after the calibration with J1744–3116. The phase self-calibration algorithm determines phase corrections for each antenna and time segment, producing a set of visibilities consistent with a point source at a fixed location.

We then reconstruct CLEAN (Högbom 1974) images corresponding to the phase self-calibrated long-baseline data on timescales of individual subscans. These images correspond to a near-perfect point source with a flux density equal to that of Sgr A* at each short time period. The relevant flux density and error estimates were obtained by fitting the CLEAN image using the CASA task imfit. The sequence of these short-time flux density measurements defines the time-dependent light curve of Sgr A* for each observing day. Finally, data corresponding to source elevations below 30° were found to be of poor quality and self-consistency and were subsequently flagged in the final data set.

An initial pass through the SMA data was performed with a custom MATLAB ¹⁵¹ based reduction pipeline, primarily responsible for preliminary flagging and bandpass calibration. Bandpass calibration was performed using various bright calibrators, given in Table 2. After these steps, the bandpass corrections were applied to the data, after which they were spectrally averaged down by a factor of 128, to a channel resolution of 17.875 MHz.

After averaging, a second round of bandpass solving was performed, and the solutions were inspected to verify that the gain corrections were consistent with unity (as the data had already been bandpass corrected). The absolute flux density scale was set by using the flux density calibrator observed closest to the time of the Sgr A* observations, also noted in Table 2, using the Butler–JPL–Horizons 2012 models, ¹⁵² on a spectral-window-by-spectral-window basis. Next, amplitude gains for individual antennas were derived using bright quasars, NRAO 530 and J1924–2914. Analysis of the gain solutions showed that the most significant trends are correlated with the elevation of the gain calibrator, consistently with known issues with antenna pointing on the SMA at very low elevations (∼15°; data corresponding to lower elevations were flagged). In light of this, gain amplitudes were interpolated using a third-order polynomial fitted based on the elevation of Sgr A*, with the median amplitude elevation-dependent correction below 25° being approximately 5%.

Due to the rapid fluctuations in the instrumental phase arising from the real-time phasing loop used in VLBI beam forming, the first three integrations (≈30 s in total) were flagged whenever the telescopes moved onto Sgr A* from a calibrator source. Additionally, due to strong line absorption, presumably arising from CN foreground absorption (see Appendix H.1. of Goddi et al. 2021), spectral channels between 226.6 and 227.0 GHz were flagged.

After amplitude-only gain calibration and the aforementioned flagging, a round of phase-only gains were derived and applied using self-calibration of Sgr A* itself. Data for these observations were collected while the array was in a compact configuration including baselines with lengths spanning ∼5–50 kλ, which at 230 GHz are sensitive to structures up to ∼20'' in size, picking up extended minispiral emission surrounding Sgr A*. Examination of the SMA data shows a strong uptick in visibility amplitudes at (u, v)-distances below 15 kλ (about 20 m). Therefore, visibilities from shorter baselines are flagged prior to self-calibration and further analysis. The remaining long-baseline data are mostly sensitive to the flux density from the unresolved Sgr A* point source; see Figure 2.

Once phase self-calibration corrections were applied, an Sgr A* light curve was generated by taking the naturally weighted vector average of all baselines, for each spectral window observed by the SMA. The resultant light curves were evaluated for large fluctuations in amplitude, under the assumption that, over a 30 s interval, changes in the brightness of Sgr A* should be subdominant to the instrumental noise. Where fluctuations greater than 3σ were seen, the measurement in question was flagged, with the total volume of data flagged in this way amounting to ∼ 1%. Finally, to help improve the S/N of the data, they were time-averaged over 62 s intervals (six integration steps).

The data sets were correlated through a Locally Normalized Discrete Correlation Function (LNDCF), as defined by Lehar et al. (1992), which revised the standard algorithm proposed by Edelson & Krolik (1988),

$$\begin{eqnarray}&&\mathrm{LNDCF}({\rm{\Delta }}t)=\displaystyle \frac{1}{{M}_{{\rm{\Delta }}t}}\sum _{i,j}\displaystyle \frac{\left({a}_{i}-{\overline{a}}_{{\rm{\Delta }}t}\right)\left({b}_{j}-{\overline{b}}_{{\rm{\Delta }}t}\right)}{\sqrt{\left({{\sigma }^{2}}_{a{\rm{\Delta }}t}-{e}_{a}^{2}\right)\left({{\sigma }^{2}}_{b{\rm{\Delta }}t}-{e}_{b}^{2}\right)}},\end{eqnarray} \tag{ 4 }$$

where a_i and b_j indicate the flux density measurements of the two compared data sets, e_a and e_b refer to the estimated measurement errors, and M_Δt represents the number of data pairs contributing to the lag bin, Δt. The flux density means and standard deviations, ${\overline{a}}_{{\rm{\Delta }}t},{\overline{b}}_{{\rm{\Delta }}t},{\sigma }_{a{\rm{\Delta }}t},{\sigma }_{b{\rm{\Delta }}t}$, are calculated for each lag, Δt, using exclusively the flux density measurements that contribute to the calculation of the LNDCF(Δt). As a comparison between the data sets, we compute the LNDCF(0) ≡ LNDCF₀, presented in Table 3.

The correlation is generally high between the two ALMA pipelines, A1 and A2, with values higher than 0.8 on each individual day and band. It remains larger than 0.7 if we consider full light curves formed by joining the individual days. Similarly, there is a rather high correlation between the SMA data set and the ALMA pipeline A1, reaching above 0.75 in all cases. The correlation is less satisfactory between the A2 pipeline and the SMA, dropping below 0.5 in some cases on 2017 April 6 and 11, but remaining high for the longest and most informative light curve from 2017 April 7. Some discrepancies between A2 and SM can be directly seen in Figure 3. Note that the ALMA–SMA correlation is calculated only in the short overlapping time of 2–3 hr, when Sgr A* is seen at low elevation by both ALMA (where it is setting) and the SMA (where it is rising), contributing additional difficulty to constraining systematic errors.

Apart from the correlation, which informs us about the consistency of the variable component, we are also interested in the consistency of the absolute flux density scale. We characterize it by comparing the median flux density in the overlapping observing periods. These results are summarized in Table 4. The systematic uncertainties of the absolute flux density scaling can be as large as 10% and vary between the days, although the ratios are quite consistent between the bands. These uncertainties do not affect relative variability metrics such as the light-curve modulation index, σ/μ, defined as the standard deviation divided by the mean, given in Table 1 (see also Section 4.2). Yet another way to quantify the differences between the data pipelines is through a mean flux density absolute difference between A2 and A1. In terms of this metric, the mean A1–A2 light-curve consistency is 3.7% on 2017 April 6, 2.7% on 2017 April 7, and 16.0% on 2017 April 11, the latter being strongly dominated by the constant offset in the flux density measurements.

We see that the overall discrepancy between light curves produced by different pipelines can be substantial. In particular, it can be significantly larger than the formal level of the thermal error in the data. Hence, we conclude that the errors are strongly dominated by the calibration systematics, which we attribute predominantly to the imperfections in the gain calibration of the individual telescopes participating in the connected-element arrays, manifesting themselves as slowly varying differences between the pipelines. If the pipelines were considered fundamentally equal, the consistency metrics provided in this section could serve as a proxy for quantifying systematic errors. However, since the A1 pipeline relies less heavily on a priori sensitivity estimates than the other ones, it is expected to be more robust against the relevant sources of corruption. In Figure 3, we observe the decorrelation between A1 and A2 to increase toward the end of the observing tracks, which suggests that these inconsistencies are related to the low source elevation, which is exactly when more severe gain-related corruptions are to be expected. In the subsequent analysis, we stress the A1 pipeline results in particular, supplementing them with the A2 and SM results, allowing us to assess our confidence in the obtained result.

The LNDCF₀ coefficient computed between the different frequency bands within the same pipeline is remarkably high in all cases, above 0.99. This consistency across the frequency bands can be seen in Figure 4. We also verify the ratio of medians in the overlapping observing periods, using the HI band as a reference; see Table 5. We notice that the ratio is very close to unity, which is consistent with the flat mm spectrum of Sgr A* (Marrone et al. 2006; Bower et al. 2015) and expected given the narrow fractional band, Δν/ν ≲ 0.1. There is a persistent systematic effect of 4% missing flux density in the LO frequency band (227.1 GHz), seen in both of the ALMA pipelines; see Table 5 and the right panel of Figure 4. This could be a systematic processing/scaling error shared by both of the ALMA reduction pipelines, or an effect of absorption in the LO band. A similar effect is not seen in the SMA data, for which spectral channels possibly affected by CN absorption were flagged within the LO band (Section 2.3). However, the absorption alone was estimated to be too small to be responsible for a 4% effect (Appendix H.1. of Goddi et al. 2021). As a result of this discrepancy, we refrain from using the ALMA LO band for applications such as the spectral index estimation.

Zoom In Zoom Out Reset image size

We model the frequency dependence of the flux density with a power law, F_ν ∝ ν^α , thus defining the spectral index as α. Subsequently, we compute α for each pair of simultaneous flux density measurements in the bands (B1, H i) and (B2, H i). We show the results with sample standard deviation error bars in Figure 5. We conclude that the spectral index measured between 213.1 and 229.1 GHz is consistent with zero, α₂₂₀ = 0.0 ± 0.1. Figure 5 implies that the calibration-related systematic uncertainties and short-timescale fluctuations of the spectral index dominate the associated error budget. In Appendix B, we confirm these findings with the full-bandwidth SMA data analysis.

Zoom In Zoom Out Reset image size

Combining our flux density measurements at 220 GHz of 2.4 ± 0.2 Jy with the compact flux of 2.0 ± 0.2 Jy at 86 GHz reported by Issaoun et al. (2019) based on semi-simultaneous observations on 2017 April 3, we find a spectral index of α₁₅₀ = 0.19 ± 0.13 at ν₀ = 0.5 × (86 + 220) ≈ 150 GHz. Hence, we find a small positive spectral index at about 150 GHz that becomes consistent with zero at about 220 GHz. These findings are generally consistent with a flat spectral index at mm wavelengths reported by Bower et al. (2015) and Iwata et al. (2020), as well as with our broad understanding of the Sgr A* spectral energy distribution, with a flattening spectrum in the mm approaching a peak in the submm ("submm bump"; Zylka et al. 1995; Melia & Falcke 2001). The mean light-curve spectral index may be an important discriminant of the theoretical models of Sgr A* (Ricarte et al. 2022), as this quantity is sensitive to physical properties such as temperature, magnetic field strength, optical depth, and the electron distribution function.

One can also resolve the measured spectral index as a function of time, obtaining the results presented in Figure 6. These measurements show large fluctuations of the spectral index and swings in a range between −0.2 and 0.1 on a timescale of ∼1 hr. This can be interpreted as rapid fluctuations of the effective optical depth of the compact system, possibly related to the turbulent character of the accretion flow. Interestingly, both pipelines indicate that α was more negative immediately after the 2017 April 11 X-ray flare (ALMA observations begin at 9.0 UT, about 10–15 minutes after the peak of the X-ray flare reported by Chandra in Paper II; see also Figure 3), reaching −0.23 ± 0.05 and subsequently recovering to values consistent with zero on a timescale of 1–2 hr. This suggests an increased contribution of the optically thin component to the total intensity immediately after the X-ray flare. Indeed, since the synchrotron self-absorption decreases with decreasing magnetic field, B, and increasing plasma temperature, T _e (Rybicki & Lightman 1979), a flaring event injecting energy of the magnetic field into electrons through magnetic reconnection (Yuan et al. 2003) is expected to reduce the effective optical depth of the system. Additionally, we note that the first scan by ALMA, which marginally overlaps with the X-ray flare, indicates a decrease in the 1.3 mm emission, while all subsequent scans for the next 2 hr show a growing flux density, in total by about 60%. Such an evolution of the flux density and spectral index suggests a particle acceleration event where magnetic reconnection heats up electrons to a power-law distribution (Guo et al. 2014; Sironi & Spitkovsky 2014; Werner et al. 2015), thus shifting the emission to near-IR and X-ray wavelengths and causing an inverted spectrum (i.e., α < 0). As the electrons cool down radiatively, and subsequently the reconnection layer powering the flare depletes (Ripperda et al. 2021), the optically thin emission shifts back to mm and radio wavelengths (e.g., Brinkerink et al. 2015), and the source eventually settles back to the state before the flare.

Zoom In Zoom Out Reset image size

An elevated X-ray activity was also reported in Chandraobservations on 2017 April 7 at 11–13 UT (Paper II). Here we see that this X-ray event was accompanied by a decreased spectral index period in our A1 pipeline data (as seen in Figure 6) and a total flux density decrease at 11–13 UT in the A1 pipeline and the SMA observations. This was then followed up by a flux density recovery seen in the SMA data around 14 UT (as seen in Figure 3). All of these observations further strengthen the presented interpretation.

In Table 6, we present the previously published Sgr A* light-curve data sets at frequencies close to 230 GHz (that is, closest in frequency to our H i band). We only consider observations with radiointerferometric arrays, where reliable extraction of the compact source light-curve component is feasible. Compared to data sets published in this paper, summarized in Table 1, the archival data sets typically have lower cadence and a far lower number of collected data points. Thus, more reliable studies of the source variability, particularly on short timescales, are enabled by our new data sets.

All data sets given in Tables 1 and 6, spanning a total period of about 14 yr, can be divided into 18 observing epochs no longer than 16 days, where the EHT observations constitute a single epoch of 2017 April 5–11. Normalized histograms of the 230 GHz flux density observed in these epochs are shown in the top left panel of Figure 7. The flux density remains remarkably consistent across all these epochs, with all measurements in agreement with 4.0 Jy within about 50%.

Zoom In Zoom Out Reset image size

We also show a (differently normalized) generalized λ distribution (GλD; Freimer et al. 1988) fit to all of the 2005–2019 data sets, computed using the gldex package (Su 2007). It approximates the full distribution of the Sgr A* flux density at 230 GHz across all of the observing epochs. The GλD fit corresponds to flux density values within ${3.24}_{-0.60}^{+0.68}$ Jy at 68% confidence, indicating a weak positive tail driven primarily by the record high flux densities observed by Fazio et al. (2018). All measurements given in Tables 1 and 6 are also shown in the bottom panel of Figure 7 as a function of the observing date.

These ranges are also consistent with Sgr A* monitoring with ALMA and SMA in 2013 June–2014 November presented in Bower et al. (2015). ¹⁵⁴ The relative calmness of Sgr A* is in strong contrast to the X-ray and IR behavior, where flux densities may vary by orders of magnitude during the flaring events (Porquet et al. 2003; Do et al. 2019b). We notice that the 2017 April epoch is characterized by the lowest mean flux density among all 2005–2019 observations. Within several hours of a single observing epoch, the mm flux density of Sgr A* may fluctuate by ∼1 Jy.

We quantify the variability with the modulation index, σ/μ, corresponding to the ratio between the signal standard deviation, σ, and its mean value, μ. It is related to the rms–flux relation, abundantly used in IR and X-ray studies of variability. We notice, however, that unlike in the IR (e.g., Gravity Collaboration et al. 2020), we did not find strong indications of a linear rms–flux relation in the mm light curves. For a red-noise stochastic process, which we expect to describe the variability of Sgr A* mm light curves well (Dexter et al. 2014), most variability manifests on the longest timescales. Hence, the modulation index, σ_T /μ_T , calculated from a chunk of data of a finite duration T, is biased with respect to the asymptotic modulation index, σ_∞/μ_∞. On the other hand, for a particular light-curve realization, the influence of sparse or nonuniform sampling on σ/μ is small. We verified the robustness of the σ/μ estimation against the measurement noise and irregular sampling by comparing our findings with the results of the intrinsic modulation index algorithm of Richards et al. (2011), finding an excellent agreement. The influence of the light-curve duration T can be seen in the top right panel of Figure 7, where light curves of longer duration generally exhibit a larger modulation index. The figure presents all of the observations listed in Tables 1 and 6. For data sets spanning several days, we also show the modulation index calculated for the entire campaign, corresponding to the histograms in the top left panel of Figure 7. The variability measurements are compared with the expectations from a GP model (damped random walk; black line indicating the expected values, with 1σ, 2σ, and 3σ error ranges plotted as blue bands) best-fitting the full combined 2005–2019 data set; see Section 5 for details and discussion.

The 230 GHz light curves collected in 2005–2019, in particular the high-quality EHT light curves from 2017 April, indicate rather low modulation index, σ/μ, typically below 0.10. Hence, we conclude that on 2017 April 6 and 7 the source displayed an amount of variability consistent with historical measurements. On 2006 July 17 (Marrone et al. 2008), 2015 May 14 (Fazio et al. 2018), and 2017 April 11 (this work and Paper II) increased variability metrics can be connected to flares detected in the X-ray; however, the variability enhancement is particularly clear only in the case of the 2017 April 11 observations. We expect that modulation index values above ∼0.15 seen in the top right panel of Figure 7 may possibly be outliers suffering from calibration errors—it is generally far easier to increase the apparent variability with the calibration errors than to reduce it (e.g., erroneous amplitude gains, coherence losses, pointing issues).

The modulation index measured in general relativistic magnetohydrodynamic (GRMHD) simulations was found to be generally larger than what the observations indicate (Chatterjee et al. 2021; Paper V). For comparisons between observations and simulations, a T = 3 hr ≈ 540 GM/c³ window for computing the modulation index was used in Paper V. This duration is justified by the synthetic observations decorrelation argument—separate 3 hr segments are expected to behave like statistically independent draws from the modulation index statistic. In Table 7, we give nonoverlapping values of ${\left(\sigma /\mu \right)}_{3\,\mathrm{hr}}$ from all days and sites/pipelines (three nonoverlapping samples for ALMA 2017 April 7, a single 3 hr modulation index measurement for all of the other light curves). The measurements presented in Table 7 show a factor of 2 enhancement of the 3 hr modulation index on the X-ray flare day of 2017 April 11. On the remaining days, the modulation index varies between 0.024 and 0.051, while the damped random walk model fitted to all of the 2005–2019 data sets predicts ${\left(\sigma /\mu \right)}_{3\,\mathrm{hr}}={0.03}_{-0.01}^{+0.02}$, as shown in the top right panel of Figure 7.

To investigate the possible existence of characteristic variability timescales in the Sgr A* light curves, a second-order SF analysis (Simonetti et al. 1985) has been applied to the data. The SF of a time series {x_i } = x₁, x₂,...,x_n , observed at times {t_i } = t₁, t₂,...,t_n , at time lag Δt, is defined as

$$\begin{eqnarray}&&\mathrm{SF}({\rm{\Delta }}t)=\displaystyle \frac{1}{{M}_{{\rm{\Delta }}t}}\sum _{i,j}{\left({x}_{i}-{x}_{j}\right)}^{2},\end{eqnarray} \tag{ 5 }$$

where the sum is extended to all M_Δt pairs (t_i , t_j ) for which Δt − Δt₀/2 < (t_i − t_j ) < Δt + Δt₀/2, and Δt₀ is the shortest time lag for which the SF is calculated. The SF informs us about the signal variance across a range of timescales. A noise contribution has been neglected in Equation (5), given the very high reported data S/N. For this analysis, we use the data cadence reported in Table 1 as Δt₀. Assuming that the observed variability can be described as a sum of the random error (measurement/calibration error) and the true signal, possibly resulting from a complex superposition of processes with different spectral properties and characteristic timescales, the SF is expected to show the following:

• 1.  
  A flat slope at the shortest timescales, when the random error amplitude dominates over the source signal.
• 2.  
  A steepening increase on a range of timescales for which the random error amplitude is nonnegligible compared to the signal.
• 3.  
  A steep increase with a constant slope, on timescales for which the contribution of the random error to the flux density measurements is negligible compared to the variations induced by the source. If the signal can be modeled in the spectral domain as a power law with the PSD exponent (α_PSD) steeper than ≈−1.5 but not steeper than −3 (Emmanoulopoulos et al. 2010), the SF slope (α_SF) should have a value of α_SF ≈ − (1 + α_PSD).
• 4.  
  A change of slope at the characteristic timescale of the source signal, which corresponds to 0.5/f _b , where f _b is the frequency at which the power-law PSD shows a break. In case the signal is a superposition of multiple components, each characterized by its own timescale and power-law slope, these should be reflected as slope changes in the SF.
• 5.  
  A plateau at a time lag corresponding to the maximum characteristic timescale of the source.
• 6.  
  A flat slope at larger timescales, where the SF should oscillate around a value of twice the sum of the variances of the signal and the measurement/calibration noise.

An SF analysis is prone to identifying spurious characteristic timescales resulting from random fluctuations in a finite realization of a red-noise process, possibly interacting with the sampling window function. To determine whether a detected characteristic timescale is real, it is necessary to sample at least several cycles of the variability (that is, to observe for a duration of at least several times longer than the timescale in question). The significance of timescales larger than 0.2 times the total duration of the observations is low; it slowly increases as this ratio becomes smaller. At the shorter timescales, the SF results reflect quite accurately the properties of the signal realized in the light curves. The SF slope can therefore be a faithful estimator for the PSD slope (α_PSD).

4.3.1. Estimating Intrinsic Noise

As a first step of the SF analysis, we isolated the random noise component by applying a denoising algorithm, which works as a low-pass filter with a cutoff timescale of 0.01 hr < 2GM/c³. To verify the correct separation between the source signal and the random noise contribution to the variability, we applied the SF to both the denoised signal and the noise component. In the first case, we checked that the slope at the shortest timescales follows the same trend as at the intermediate ones, where the random noise is negligible. For the noise component, we verified that the SF slope is approximately zero, in agreement with properties of white noise. The noise component becomes subdominant for timescales longer than ∼1 minute; see Figure 8. As a by-product of this step, we obtain a realistic estimate of the flux density uncertainties. These uncertainties turn out to be generally larger than the statistical errors reported in the data sets, with values on the order of ∼0.5% of the measured flux densities, or about 0.01 Jy. A cross-correlation of the ALMA noise component extracted from the two independent calibration procedures shows that for all epochs and frequencies there is no correlation between them. This result allows us to conclude that the random noise is mainly due to the calibration-specific uncertainties.

Zoom In Zoom Out Reset image size

4.3.2. SF Analysis Results

In the case of stationary signals, there is a unique relationship between the autocorrelation and the SF (see Section 5). Nevertheless, apart from the uncertainty in the stationarity assumption, studying autocorrelations separately offers a different perspective into the data. We study the signal autocorrelation using the LNDCF method (Lehar et al. 1992), Equation (4). A summary of these results is shown in Figure 9. In this plot, we indicate all of the contributing autocorrelation measurements (circles), with the running mean (colored line) and the running one standard deviation bands. In the first panel of Figure 9, we show the autocorrelation calculated using all of the available observations in each pipeline, in the H i band (the other bands are very consistent). The data indicate autocorrelations decreasing roughly on a timescale of ∼1 hr; the black dashed line corresponds to $\exp (-{\rm{\Delta }}t/1\,\mathrm{hr})$. This is significantly less than one would expect based on the results of Dexter et al. (2014). In the subsequent panels of Figure 9, we show autocorrelations for 2017 April 6, 7, and 11. The nonmonotonic structure of the autocorrelation functions is not detected confidently, given the associated uncertainties. The persistence of such features may be established with more observations, e.g., pipeline A1 results show a bump at ∼30 minutes for both 2017 April 7 and 11. This resembles the innermost stable circular orbit period for a Schwarzschild black hole with 4 × 10⁶ M_⊙ mass, but there is very little confidence in such an association at this point. Significant biases that may affect the autocorrelation measurement prevent us from drawing strong conclusions based on this analysis; see the discussion in Section 5.5.

Zoom In Zoom Out Reset image size

Several authors have recently considered signatures of multiple-path propagation of photons traveling through a strongly curved spacetime in a black hole vicinity and reaching the observer with a delay (e.g., Moriyama et al. 2019; Chesler et al. 2021; Hadar et al. 2021; Wong 2021; or Wielgus et al. 2020 for the case of exotic spacetimes of black hole mimickers). While observing a feature related to a photon shell around a black hole is a tantalizing possibility, such an observation does not seem feasible yet, given the model simplifications and the limited duration of the observations. There are no significant signatures of autocorrelations detected at the relevant time lags of ∼ 20GM/c³ ≈ 400 s in our presented data sets.

The presence of time lags between frequency bands in observations of Sgr A* has been theoretically predicted. In the optically thick regime of radio observations at frequencies below 100 GHz, theoretical models attribute those lags to the adiabatic spherical expansion of plasma blobs (van der Laan 1966; Yusef-Zadeh et al. 2008; Eckart et al. 2012), or to a bulk outflow (Falcke et al. 2009). Such lags could be detected across the spectrum, with the higher-frequency signal typically leading the lower frequency signal (Yusef-Zadeh et al. 2009; Brinkerink et al. 2015, 2021). Hints of a similar delay structure have been seen in some numerical GRMHD models of Sgr A* compact emission (Chan et al. 2015). In this theoretical framework, we could expect lags of 1–2 minutes between the HI and B1 bands, easily detectable with the cadence of the ALMA data. However, no indication of a correlation lag between the bands in any of our data sets is found, with the uncertainty no larger than 20 s. The cross-correlation function very clearly peaks at zero for all days and for both ALMA reduction pipelines, as shown in Figure 10. We interpret this lack of detectable delays as a signature that the emission region in the 213–229 GHz range is already optically thin all the way to the horizon, possibly with patches of higher optical depth material formed in the turbulent accretion flow, necessary to explain the intermediate spectral index α ≈ 0 (identified in Section 3.3). This is particularly likely given that in 2017 April Sgr A* was in a rather low mm flux density state (Section 4; Mościbrodzka et al. 2012). Historically, no delays were reported by Iwata et al. (2020) across similar frequency bands; there was also no significant delay between 230 and 345 GHz reported by Marrone et al. (2008), and no delay between 134 and 146 GHz or between 230 and 660 GHz reported by Yusef-Zadeh et al. (2009). Finally, the conclusion of a low optical depth is consistent with the interpretation of the first EHT images of Sgr A*, reported in Paper III, as the observable shadow of a supermassive black hole.

Zoom In Zoom Out Reset image size

DRW (or the Ornstein–Uhlenbeck process) is a unique Markovian stationary GP (Rasmussen & Williams 2006). Application of the DRW as a mathematical model to describe the optical variability of quasars was proposed by Kelly et al. (2009). The DRW with variance σ² is characterized by a covariance function,

$$\begin{eqnarray}&&{k}_{\mathrm{DRW}}({\rm{\Delta }}t)={\sigma }^{2}\exp \left(-\left|\displaystyle \frac{{\rm{\Delta }}t}{\tau }\right|\right),\end{eqnarray} \tag{ 6 }$$

where a characteristic timescale, τ, is a model parameter. For stationary processes there is a general relation between the covariance and the SF defined in Equation (5),

$$\begin{eqnarray}&&\mathrm{SF}({\rm{\Delta }}t)=2{\sigma }^{2}-2k({\rm{\Delta }}t).\end{eqnarray} \tag{ 7 }$$

The corresponding PSD is related to the covariance function through the Wiener–Khinchin theorem, and in the case of a DRW process it becomes

$$\begin{eqnarray}&&\mathrm{PSD}(f)=\displaystyle \frac{4{\sigma }^{2}\tau }{1+{\left(2\pi f\tau \right)}^{2}},\end{eqnarray} \tag{ 8 }$$

which corresponds to a red-noise spectrum with an index of α_PSD = − 2 in a high-frequency limit (f τ ≫ 1) and a flat white-noise spectrum at low frequencies (f τ ≪ 1). In the context of the Galactic Center, Dexter et al. (2014) modeled the Sgr A* variability at mm wavelengths with the DRW, following the procedure outlined by Kelly et al. (2009). By fitting several years of Sgr A* observations (Section 4), they identified a poorly constrained DRW timescale of τ ∼ 8 hr.

Our model differs from that of Kelly et al. (2009) and Dexter et al. (2014) only by the inclusion of an additional parameter, σ₀, representing the noise floor. Then, the full model consists of four variables,

$$\begin{eqnarray}&&\theta =(\tau ,\mu ,\sigma ,{\sigma }_{0}),\end{eqnarray} \tag{ 9 }$$

the correlation timescale τ, the mean value of the process μ, the standard deviation σ, and the noise floor σ₀. The timescale τ is not necessarily related to the timescales estimated in Section 4.3—the SF timescales may be indicative of the presence of multiple stochastic components in the real signal. Because the DRW is a Markovian process, the likelihood function for observations, {x_i } = x₁, x₂,...,x_n , observed at times {t_i } = t₁, t₂,...,t_n , can be calculated directly as

$$\begin{eqnarray}&&p\left(\{{x}_{i}\}| \theta \right)=\prod _{i=1}^{n}\displaystyle \frac{\exp \left[-0.5{\left({\hat{x}}_{i}-{x}_{i}^{* }\right)}^{2}/{\widetilde{{\rm{\Omega }}}}_{i}\right]}{{\left(2\pi {\widetilde{{\rm{\Omega }}}}_{i}\right)}^{1/2}},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\widetilde{{\rm{\Omega }}}}_{i}={{\rm{\Omega }}}_{i}+{\sigma }_{0}^{2}+{\sigma }_{i}^{2},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{x}_{i}^{* }={x}_{i}-\mu .\end{eqnarray} \tag{ 12 }$$

The indexed σ_i represents measurement uncertainties and is distinct from the estimated process standard deviation σ and the noise floor σ₀. The quantities Ω_i and ${\hat{x}}_{i}$ are calculated through an iterative procedure,

$$\begin{eqnarray}&&{\hat{x}}_{i}={a}_{i}{\hat{x}}_{i-1}+\displaystyle \frac{{a}_{i}{{\rm{\Omega }}}_{i-1}}{{\widetilde{{\rm{\Omega }}}}_{i-1}}\left({\hat{x}}_{i-1}-{x}_{i-1}^{* }\right),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{i}={{\rm{\Omega }}}_{1}\left(1-{a}_{i}^{2}\right)+{a}_{i}^{2}{{\rm{\Omega }}}_{i-1}\left(1-\displaystyle \frac{{{\rm{\Omega }}}_{i-1}}{{\widetilde{{\rm{\Omega }}}}_{i-1}}\right),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{a}_{i}=\exp \left[-({t}_{i}-{t}_{i-1})/\tau \right],\end{eqnarray} \tag{ 15 }$$

with the initial conditions

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{1}={\sigma }^{2};{\hat{x}}_{1}=0.\end{eqnarray} \tag{ 16 }$$

Note that we use a slightly different parameterization of the DRW model than Kelly et al. (2009), with σ² representing the variance of the DRW, related to their parameterization by ${\sigma }_{\mathrm{Kelly}}=\sigma \sqrt{2/\tau }$. Since the procedure outlined above allows us to explicitly compute the best-fitting DRW realization for a given vector of parameters, θ, we can assess the fit quality by computing the reduced-χ² statistic for the residuals,

$$\begin{eqnarray}&&{\chi }_{n}^{2}=\displaystyle \frac{1}{n}\sum _{i=1}^{n}\displaystyle \frac{{\left({x}_{i}^{* }-{\hat{x}}_{i}\right)}^{2}}{{\widetilde{{\rm{\Omega }}}}_{i}}.\end{eqnarray} \tag{ 17 }$$

The DRW model fixes the high-frequency limit PSD slope to α_PSD = − 2. This is a rather strong assumption, and there are indications of a steeper PSD slope in both the context of optical variability of quasars (Mushotzky et al. 2011; Zu et al. 2013) and the variability of X-ray binaries (e.g., Tetarenko et al. 2021). The high-frequency PSD slope may be a relevant parameter to extract, less affected by the sampling and observation duration limitations than the timescale τ, and having the potential to constrain theoretical models of Sgr A*. Moreover, observations presented in this paper sample the high-frequency regime, relevant for constraining α_PSD uniquely well. Hence, we employ a more general statistical model of a GP with a Matérn covariance function (see, e.g., Rasmussen & Williams 2006),

$$\begin{eqnarray}&&{k}_{\mathrm{Mat}}({\rm{\Delta }}t)={\sigma }^{2}\displaystyle \frac{{2}^{1-\nu }}{{\rm{\Gamma }}(\nu )}{\left(\sqrt{2\nu }\displaystyle \frac{{\rm{\Delta }}t}{\tau }\right)}^{\nu }{K}_{\nu }\left(\sqrt{2\nu }\displaystyle \frac{{\rm{\Delta }}t}{\tau }\right),\end{eqnarray} \tag{ 18 }$$

where K_ν is the modified Bessel function of the second kind. The parameter ν defines the order of the Matérn process and subsequently controls the smoothness of the resulting curve. The PSD of the Matérn process is

$$\begin{eqnarray}&&{\mathrm{PSD}}_{\mathrm{Mat}}(f)\propto {\left[1+\displaystyle \frac{{\left(2\pi f\tau \right)}^{2}}{2\nu }\right]}^{-(\nu +1/2)},\end{eqnarray} \tag{ 19 }$$

so in the high-frequency limit f τ ≫ 1 we find the PSD slope of α_PSD = − 2ν − 1. The DRW is recovered as a special case of the Matérn process with ν = 0.5. As an arbitrary ν, the Matérn covariance represents a non-Markovian process, and the likelihood cannot be evaluated explicitly as in the case of the DRW. Instead, we evaluate it numerically using the Stheno library. ¹⁵⁵

Given the low computational cost of the DRW model fitting, we were able to perform a survey of different modeling parameters, such as the type and range of priors, subsets of data to be used, and treatment of the systematic uncertainties and the noise floor. Our general conclusion is that the timescale τ cannot be well constrained and its posterior distributions are dominated by the assumed priors. As an example, Dexter et al. (2014) used log-uniform priors, reducing the distribution tails for large τ. We find that for our data sets τ remains poorly constrained, and with uniform priors very large timescales are permitted. As noted by Kozłowski (2017), the duration of the light curve needs to be significantly longer than the timescale τ to constrain it reliably. The duration and sampling of the 2017 April data may not be sufficient. On the other hand, when fitting data spanning several years (such as in the case of Dexter et al. 2014), one needs to consider whether the underlying process can be assumed to be stationary on such long timescales (e.g., as a consequence of the mass accretion rate modulation).

As a result of the DRW survey, we selected the following set of priors, π(θ):

$$\begin{eqnarray}&&\begin{array}{l}\tau \ (\mathrm{hr})\sim {{ \mathcal N }}_{{\rm{T}}}(0,8),\\ \mu \ (\mathrm{Jy})\sim {{ \mathcal N }}_{{\rm{T}}}(2,1),\\ \sigma \ (\mathrm{Jy})\sim {{ \mathcal N }}_{{\rm{T}}}(0,1),\\ {\sigma }_{0}\ (\mathrm{Jy})\sim { \mathcal U }(0.0,0.1),\end{array}\end{eqnarray} \tag{ 20 }$$

where ${{ \mathcal N }}_{{\rm{T}}}(a,b)$ is a normal distribution of mean a and standard deviation b truncated to positive values, and ${ \mathcal U }(a,b)$ is a uniform distribution with a range between a and b. For the Matérn model fitting, we adopt identical priors as given by Equation (20), with an additional prior on the PSD slope, α_PSD,

$$\begin{eqnarray}&&{\alpha }_{\mathrm{PSD}}=-2\nu -1\sim { \mathcal U }(-9,-1).\end{eqnarray} \tag{ 21 }$$

An overview of the fitting results for different data subsets is shown in Table 10, where the maximum likelihood (ML) estimator parameters are given, along with the 68% confidence intervals. These results establish a good consistency between bands, less so between the pipelines. The estimated noise floor, σ₀, is comparable to the noise amplitudes estimated in Section 4.3.