# Rapid Variability of Sgr A* across the Electromagnetic Spectrum

G. Witzel Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121, Bonn, Germany  G. Martinez Department of Physics and Astronomy, University of California, Los Angeles, Box 951547, Los Angeles, CA 90095-1547, USA  S. P. Willner Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138 USA  E. E. Becklin  Department of Physics and Astronomy, University of California, Los Angeles, Box 951547, Los Angeles, CA 90095-1547, USA  SOFIA Science Center, Moffett Field, CA USA  H. Boyce  Department of Physics, McGill University, 3600 University St., Montreal, QC H3A 2T8, Canada  McGill Space Institute, McGill University, Montreal, QC H3A 2A7, Canada  T. Do Department of Physics and Astronomy, University of California, Los Angeles, Box 951547, Los Angeles, CA 90095-1547, USA  A. Eckart I. Physikalisches Institut, Universität zu Köln Zülpicher Str. 77 D-50937 Köln, Germany  Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121, Bonn, Germany  G. G. Fazio Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138 USA  A. Ghez Department of Physics and Astronomy, University of California, Los Angeles, Box 951547, Los Angeles, CA 90095-1547, USA  M. A. Gurwell Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138 USA  D. Haggard Department of Physics, McGill University, 3600 University St., Montreal, QC H3A 2T8, Canada  McGill Space Institute, McGill University, Montreal, QC H3A 2A7, Canada  R. Herrero-Illana European Southern Observatory (ESO), Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile  J. L. Hora Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138 USA  Z. Li (李志远) School of Astronomy and Space Science, Nanjing University, 163 Xianlin Avenue, Nanjing 210023 Jiangsu, China  J. Liu (刘俊) Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121, Bonn, Germany  N. Marchili Istituto di Radioastronomia—INAF, Via Piero Gobetti 101, 40129, Bologna, Italy  Mark R. Morris Department of Physics and Astronomy, University of California, Los Angeles, Box 951547, Los Angeles, CA 90095-1547, USA  Howard A. Smith  Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138 USA  M. Subroweit I. Physikalisches Institut, Universität zu Köln Zülpicher Str. 77 D-50937 Köln, Germany  J. A. Zensus Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121, Bonn, Germany  I. Physikalisches Institut, Universität zu Köln Zülpicher Str. 77 D-50937 Köln, Germany

(Received ??, 2020; Accepted ??, 2020)

###### Abstract

Sagittarius A* (Sgr A*) is the variable radio, near-infrared (NIR), and X-ray
source associated with accretion onto the Galactic center black hole. We have
analyzed a comprehensive submillimeter (including new observations
simultaneous with NIR monitoring), NIR, and 2–8 keV dataset. Submillimeter
variations tend to lag those in the NIR by $\sim$30 minutes. An approximate
Bayesian computation (ABC) fit to the X-ray first-order structure function
shows significantly less power at short timescales in the X-rays than in the
NIR. Less X-ray variability at short timescales combined with the observed
NIR–X-ray correlations means the variability can be described as the result of
two strictly correlated stochastic processes, the X-ray process being the low-
pass-filtered version of the NIR process. The NIR–X-ray linkage suggests a
simple radiative model: a compact, self-absorbed synchrotron sphere with high-
frequency cutoff close to NIR frequencies plus a synchrotron self-Compton
scattering component at higher frequencies. This model, with parameters fit to
the submillimeter, NIR, and X-ray structure functions, reproduces the observed
flux densities at all wavelengths, the statistical properties of all light
curves, and the time lags between bands. The fit also gives reasonable values
for physical parameters such as magnetic flux density $B\approx 13$ G, source
size $L\approx 2.2R_{S}$, and high-energy electron density $n_{e}\approx
4\times 10^{7}$ cm-3. An animation illustrates typical light curves, and we
make public the parameter chain of our Bayesian analysis, the model
implementation, and the visualization code.

††journal: ApJS††facilities: Spitzer(IRAC),ALMA,SMA,Chandra(ACIS)

##  1 Introduction

Since the discovery of rapid flaring of Sagittarius A* (Sgr A*) in the X-rays
and near infrared (NIR) in the early 2000s (Baganoff et al., 2001; Genzel et
al., 2003; Ghez et al., 2004), many monitoring programs have been executed to
understand the origin and properties of the variable emission. As a result,
more than 70 papers have been published describing observations and modeling
light curves in the submillimeter (submm), NIR, and X-rays. Most of these
publications have focused on statistical analyses of flux densities and timing
properties, multi-wavelength observations and modeling of the spectral energy
distribution (SED), and NIR and submm polarization. (See reviews by Genzel et
al. 2010 and Morris et al. 2012, the comprehensive Witzel et al. 2018
introduction, and references in their paper.) This fascination with Sgr A*’s
variability has a reason: light-crossing-time arguments link rapid changes in
flux density to spatial scales that—until recently—were not accessible
otherwise. For example, Dodds-Eden et al. (2009) found sudden flux density
changes of a factor ${\sim}2$ in less than 47 seconds, and Do et al. (2019)
saw changes of a factor ${\sim}9$ in less than 2 minutes. These times
correspond to spatial scales of $1.2R_{S}$ and $3R_{S}$, respectively.111Here
$R_{S}$ means Schwarzschild radius, $R_{S}=1.23\times 10^{10}$ m for mass
$M=4.15\times 10^{6}$ M☉ (Gravity Collaboration et al., 2019). The scales
mentioned are not to be taken as distances to the black hole. They are rather
upper limits on the characteristic size of the volume where the radiation
originates, e.g., of a region in the accretion disc or in a jet likely well
away from the event horizon.

Event-horizon spatial scales can now be studied by two types of large
interferometers. One is VLTI/GRAVITY, which operates in the NIR and has
observed the variable emission source moving in the plane of the sky (Gravity
Collaboration et al. 2018a, b, 2019, 2020a, 2020b). The other is very long
baseline interferometry (VLBI) at millimeter wavelengths (reviewed by Boccardi
et al. 2017). The Event Horizon Telescope (EHT) aims to image Sgr A* with a
resolution close to $R_{S}$ at 1.3–3.5 mm wavelengths (Dexter et al. 2014;
Tilanus et al. 2014; Lu et al. 2018). The detection of circular trajectories
of the center of light during flares of Sgr A* by VLTI/GRAVITY and
corresponding loops in the Stokes $(Q,U)$ plane (Gravity Collaboration et al.
2018b) suggests compact source structure of ${<}5R_{S}$ on an orbit around an
average position at ${\sim}9R_{S}$ (Gravity Collaboration et al. 2020a). At
1.3 mm wavelength, VLBI studies also found compact structure with upper limits
on the intrinsic source size of ${\lesssim}4R_{S}$ (Doeleman et al. 2008),
likewise with some evidence that the center of light is not centered on the
black hole itself.

In light of the interferometric results, the variability data obtained over
the last two decades are valuable as a complementary source of information
about the physical processes at event-horizon scales. Many multi-wavelength
campaigns in the submm, NIR, and the X-rays have been organized in the hope of
determining—or at least constraining—the underlying radiative processes
(Baganoff et al. 2001; Eckart et al. 2004, 2006; Gillessen et al. 2006; Yusef-
Zadeh et al. 2006a, b; Eckart et al. 2008a, b; Marrone et al. 2008; Yusef-
Zadeh et al. 2008; Dodds-Eden et al. 2009; Yusef-Zadeh et al. 2009; Trap et
al. 2011; Eckart et al. 2012; Yusef-Zadeh et al. 2012; Haubois et al. 2012;
Mossoux et al. 2016; Rauch et al. 2016; Ponti et al. 2017).

Several models have been proposed to explain the variability at different
wavelengths. All models assume that the NIR is dominated by optically thin
synchrotron radiation, but they differ in the mechanism for the X-ray
emission. Some make the case for optically thin synchrotron radiation with a
cooling break to explain the X-rays (Dodds-Eden et al. 2009; Ponti et al.
2017), while others suggest synchrotron self-Compton scattering
(synchrotron–SSC) (Dodds-Eden et al. 2009; Eckart et al. 2012; Mossoux et al.
2016) or inverse Compton scattering of photons by a second population of
electrons (Dodds-Eden et al. 2009). While some authors have claimed evidence
for a close relation between submm and NIR variability (e.g., linking these
bands through the evolution of the optical depth caused by adiabatic
expansion), the phenomenology and degree of correlation remains inconclusive.
Dexter et al. (2014), for example, found a submm variability timescale of
$\sim$8 hours, significantly longer than in the NIR, and concluded that
different mechanisms might be creating the variability in the two wavelength
regimes. Eckart et al. (2006), Bremer et al. (2011), and Witzel et al. (2018)
found evidence for an (exponential) cooling cutoff of the SED at NIR
frequencies. If the cooling-cutoff energy varies with source luminosity, that
could explain spectral index changes as a function of flux density. Dodds-Eden
et al. (2010) developed the first time-dependent models for Sgr A*
variability, integrating the differential equations for the electron energy
distribution under injection and escape, resulting in sequences of SEDs.

While past studies are very informative about many basic properties of the
covariability of the flux densities at different wavelengths, multi-wavelength
campaigns are difficult to organize and have to rely on favorable weather and
operational conditions at all observing sites and on the luck of Sgr A*
varying above the detection limit during the campaign. As a result, only a
small number of simultaneous multi-wavelength observations are available
today. Many of these have been obtained during 24- and 16-hour observations
with Spitzer/IRAC combined with some of Keck, Chandra, SMA, and ALMA (Fazio et
al. 2018; Witzel et al. 2018; Boyce et al. 2019). Open questions still include
the exact nature of the NIR–submm and NIR–X-ray correlations and the power
spectral density (PSD) of the X-ray variability.

Despite the critical role played by the variable component, it does not
account for all of Sgr A*’s emission. The radio emission at wavelengths
$\gtrsim$3 mm has little or no variability (Genzel et al., 2010), and there is
a non-varying, spatially extended (size?) X-ray component (Baganoff et al.,
2003; Xu et al., 2006). The angular size of the constant radio component is
unknown because of interstellar scattering. This paper addresses only the
variable component of emission. The other components must come from a separate
process or processes, and the spatial extent shows that the constant X-ray
component at least must originate in a different volume of the source than the
location giving rise to the variable emission.

The goal of the present work is to put the statistical properties of Sgr A*’s
variability at all wavelengths into the context of a single radiative model by
using available long-duration observations to best advantage and without
requiring simultaneity of the individual observations. This approach has
already been developed by Witzel et al. (2018) in modeling the covariability
of 2.2 and 4.5 µm. Given the clear evidence of structure with size
${<}5R_{S}$, we consider a one-zone model. The model combines a self-absorbed
synchrotron spectrum with cooling cutoff and synchrotron–SSC scattering. This
model choice is partly motivated by the timing properties of the correlated
NIR and X-ray light curves, which cannot easily be reproduced by the other
models mentioned above. Additionally, Eckart et al. (2012), Mossoux et al.
(2016), and Subroweit et al. (2020) showed convincingly that peak fluxes of
simultaneously observed NIR and X-ray flares as well as the flux density
distributions are well described by a synchrotron–SSC model. Eckart et al.
(2012) pointed out that bright, compact synchrotron sources show the precise
conditions to exhibit self-absorption and self-scattering. In their model, the
observed pairs of flux density peaks suggest source sizes of a few $R_{S}$ and
magnetic flux densities of a several tens of gauss. A source of this size and
magnetic flux density shows a flux density spectrum that peaks in the submm
(naturally contributing to variability at submm wavelengths as well), and
X-ray photons are mainly up-scattered submm photons. In contrast to earlier
studies, our analysis takes into account the auto-correlation and cross-
correlation properties of the light curves and the effects of the cooling
cutoff in the NIR by fitting a time-dependent, analytic model to the body of
submm, NIR, and X-ray data.

This paper has two parts. The first analyzes—for the first time—the PSD of the
X-ray variability. The PSD model is based on a generic statistical (i.e., not
physically motivated) flux density model. The paper’s second part, motivated
by the PSD analysis, proposes a single-zone radiative model and analyses the
variability at submm, NIR, and X-ray wavelengths. Both parts discuss prior
distributions, a model, and the resulting posteriors of the model parameters.
Section 2 describes the datasets used in this analysis. Section 3 shows
evidence that at the shortest variability timescales, power in the X-ray PSD
appears suppressed compared to the NIR PSD. Further, the correlation between
the NIR and the X-rays (characterized by strict correlation in arrival time
and lack of correlation in flux density levels) can be understood as the
effect of the difference between the two PSDs. Section 4 shows that the
separation of timescales is a natural result of the equation for SSC flux
densities (Marscher, 1983) if the source of the fast NIR variability is the
fast-varying cooling cutoff close to NIR frequencies. We additionally discuss
modes of co-development of the synchrotron self-absorption with NIR and X-ray
flux densities and fit a simple, semi-analytic version of the model to the
observed first-order structure functions in all bands. Section 5 presents the
model-fit results including animations. Section 6 discusses the merits and
shortcomings of the model, and Section 7 summarizes our findings.

As a convention, $H$-, $K$-, $L$-, and $M$-band refer to the respective NIR
bands. Radio bands will be denoted by their central wavelengths. We use the
term spectral energy distribution (SED) even when the quantity is expressed as
flux density ($S_{\nu}$) rather than energy ($\nu S_{\nu}$).

Table 1: New SMA and ALMA epochs Date | Start Time | Stop Time | BaselinesaaBaselines $<$40 k$\lambda$ were not used for the final flux densities in order to avoid contamination by extended structure. | # Ant. | Tuning | BWbbEffective spectral bandwidth including all polarizations, not necessarily continuous. | Calibrators
---|---|---|---|---|---|---|---
UT | UT | UT | k$\lambda$ |  | GHz | GHz |
SMA |  |  |  |  |  |  |
2014 June 18 | 07:23:48 | 13:28:47 | 14.9–99.4 | 8 | 345 | 8 | Neptune, NRAO 530
2015 May 14 | 09:52:41 | 15:52:24 |  | 6 | 237 | 8 | Titan, Callisto, NRAO 530
2016 July 13 | 05:34:48 | 12:04:48 | 47–412 | 7 | 236 | 16 | NRAO 530, J1924$-$292ccFlux density calibration was based on secondary calibrations of the gain calibrators listed.
2017 July 16 | 04:46:05 | 12:27:14 |  | 8 | 228 | 32 | NRAO 530, J1924$-$292
2017 July 26 | 05:19:43 | 11:22:07 |  | 7 | 228 | 32 | NRAO 530, J1924$-$292
ALMA |  |  |  |  |  |  |
2016 July 12 | 22:58:12 | 03:38:14 | 9.4–674.4 | 40- | 232 | 7.45 | Titan, PKS 1741$-$312
2016 July 18 | 23:14:51 | 24:10:34 | 9.4–674.4 | 40- | 232 | 7.45 | Titan, PKS 1741$-$312
2016 July 19 | 02:05:37 | 05:43:36 | 9.4–674.4 | 40- | 232 | 7.45 | Titan, PKS 1741$-$312

##  2 The Data

Our model is based on a coanalysis of a comprehensive multi-wavelength
variability dataset of Sgr A*. With the exception of two of the ALMA and four
of the SMA observations, all data were published before. While comprehensive,
this dataset is not complete, but it uses the large NIR and X-ray samples
analyzed in previous publications because they are statistically well
characterized already.

Figure 1: Observed light curves for Sgr A*. The data are presented without
observational gaps between the epochs, but each epoch is in a different color.
Gap durations are hours to years, but they should not affect the present
analysis. Time scales differ among panels. To the accuracy of the
calibrations, all flux densities are for a point source, excluding extended
structures, and all data are presented as observed with no correction for
interstellar extinction. Top panel shows 2.2 µm with the vertical dashed line
separating VLT data on the left from Keck data on the right. Second panel
shows 4.5 µm Spitzer/IRAC data in one-minute bins. Third panel shows
Chandra/ACIS-S/HETG data, and fourth panel shows Chandra/ACIS-I data on the
left and Chandra/ACIS-S data on the right. X-ray data are shown with one-
minute binning and no pileup correction. The different Chandra instruments
cover nearly the same energy range but with different sensitivities, which are
accounted for in the analysis. Fifth panel shows 345 GHz data from APEX/LABOCA
on the left and SMA on the right, and the bottom panel shows 230 GHz data from
SMA on the left and ALMA on the right.

###  2.1 Near Infrared

The NIR dataset used in this analysis is a combination of the extensive 2.12,
2.18, and 4.5 µm data of Witzel et al. (2018) and the 2.12 µm dataset of Do et
al. (2019). The original publications explain the reduction and statistics.
(For the historic $K$-band dataset see also Witzel et al. 2012 and Meyer et
al. 2014, and for the fundamentals of observing Sgr A* with Spitzer/IRAC see
Hora et al. 2014.) Our dataset contains eight 24-hour epochs of Sgr A* at 4.5
µm with IRAC on the Spitzer Space Telescope, 93 epochs of 2.18 µm data from
Naos Conica at the Very Large Telescope, 34 epochs of 2.12 µm data from the
NIRC2 camera at the Keck Observatory, in total 95 307 NIR measurements. Figure
1 shows the light curves of these three observatories. The average cadence of
the Spitzer data is 8.4 s, and typical cadences of the VLT/Keck light curves
are about one image per 1.1–1.2 minutes with integration times of 30–40 s and
28 s, respectively. The typical uncertainties of the individual data points
are 0.66 mJy for Spitzer, 0.033 mJy for the VLT, and about 0.017 mJy for the
Keck data. The light curves of all three instruments are prone to
contamination with an additive background flux density level from the resolved
and unresolved stellar population at the Galactic center. The Spitzer data
are—by nature of the applied data reduction algorithm—a differential
measurement with an arbitrary zero point. The VLT and Keck data show a typical
photometric offset of about 0.06 mJy and 0.03 mJy, respectively. These offsets
play an important role in spectral index measurements but do not affect the
first-order structure function, which quantifies differences in flux densities
rather than absolute levels.

Figure 2: Simultaneous submm and $4.5\micron$ datasets obtained during the
Spitzer campaigns. Blue lines show the 4.5 µm Spitzer light curves with one
minute binning. Flux densities are as observed, not corrected for reddening.
Purple shows 230 GHz ALMA light curves, and orange shows 230 GHz and 345 GHz
SMA light curves as labeled. Times are arbitrary but matched to heliocentric
for all curves. The 345 GHz SMA light curve of 2014 Jun 17 is the same one
published by Fazio et al. (2018). Other submm data are presented here for the
first time.

###  2.2 Chandra Data

Sgr A* has been observed often by Chandra starting in 1999. We include all
data available from the Chandra archive through 2017 having aimpoint within 1
arcmin of Sgr A*. This ensured an optimal point-spread function and best
photometric performance for a compact source like Sgr A*. The Chandra
instruments and observing modes have changed over time:

  * •

In 1999–2011, all data were taken with ACIS-I, and there were two additional
ACIS-I observations taken in 2013. These 49 ACIS-I exposures sum to 1.5 Ms.
Pileup affects these observations.

  * •

In 2012, in the framework of the X-ray Visionary Program (XVP), the data were
taken with the ACIS-S/HETG. There are three additional observations taken with
ACIS-S/HETG in 2013. These 41 ACIS-S/HETG exposures sum to ${\sim}$3 Ms. Only
the zeroth-order image is used here. Pileup is small but not always zero.

  * •

From 2013 to 2017, almost all observations were done with ACIS-S (no grating)
in sub-array mode, making them basically unaffected by pileup. The 39 ACIS-S
exposures sum to $\sim$1.4 Ms. However, during 2013 and 2014 (the first 25
epochs), the magnetar PSR J1745$–$2900 (Coti Zelati et al., 2017) contributed
significantly to the flux. Because this varying source contributed additional
photon noise, we ignored all ACIS-S epochs prior to 2015-10-21. This leaves
about a two-year gap between the last ACIS-I observation (2013) and the first
ACIS-S observation (2015-10), when the magnetar became dim enough not to
affect our photometry. The 14 good epochs give $\sim$0.6 Ms of data.

All told, we have 103 epochs and 5.3 Ms of 2–8 keV Chandra data.

We reduced the Chandra data ourselves to guarantee consistency across all
epochs. All archival data were downloaded and reprocessed using the Chandra
Interactive Analysis of Observations software package (CIAO v4.10) and the
Chandra Calibration Database (CALDB v4.7.8) following the standard procedures
as outlined by Zhu et al. (2019), who used the same Chandra dataset to study a
candidate parsec-scale jet from Sgr A*. Photons were extracted from within a
radius of 1$\farcs$25 from the best-guessed centroid of Sgr A* and within the
2000–8000 eV range to be consistent with most previous work. This gave a total
of 34 000 photons in the useful epochs. We used unbinned data of the arrival
times of individual photons (counts) and corrected the times to the barycenter
of the solar system.

Figure 1 shows the light curves of all three modes. The effective area of each
mode is a weighted mean over the 2–8 keV band, and it assumes an incident
source spectrum. If we changed the assumed spectral model, the absolute values
of the effective area would change, but the relative values among observations
using the same detector should be insensitive to the model. The relative
values between different detectors (I versus S) could be more sensitive but
are assumed to be constant as well.

###  2.3 APEX Data

We used 32 epochs (6641 minutes) of 345 GHz data from the LABOCA bolometer at
the APEX telescope (Subroweit et al., 2017). This dataset was generated using
on-the-fly mapping, resulting in fully sampled maps of
0$\fdg$5$\times$0$\fdg$17 with 280 s integration time. The data were taken
over the course of seven years and have a typical cadence of about 8 minutes.
An average map was created for each epoch and, after subtracting a Gaussian
point source at the position of Sgr A*, subtracted from each individual image.
The flux density of Sgr A* was derived by modeling a Gaussian source at the
position of Sgr A* in each residual image and using two secondary calibrators
(G10.62$-$0.38, IRAS 16293$-$2422). The relative uncertainty of the flux
density calibration is about $4\%$ or 0.1 Jy, while the absolute uncertainty
is expected to be of order $15\%$. Eckart et al. (2008a); García-Marín et al.
(2011) and Subroweit et al. (2017) gave detailed descriptions of the data
reduction, calibration, and atmospheric opacity monitoring. The upper panel of
Figure 1 shows the resulting 345 GHz light curve.

###  2.4 SMA Data

We include one published (365 minutes—Fazio et al. 2018) epoch of 343 GHz data
and four new epochs (1573 minutes) of 230 GHz data from the Submillimeter
Array (SMA). Observation details are in Table 1. The SMA was operated in a
dual-receiver polarization track with double-sideband observations using
sideband separation implemented in the correlator. The continuum visibility
was calculated by averaging the two same-sense polarization signals. Final
flux-density measurements were determined by vector-averaging the measured
visibility data for instantaneous baselines greater than 35–40 k$\lambda$ to
filter out large-scale emission structure around Sgr A*. Water vapor was
$\sim$1.5 mm for most epochs but nearer 0.9 mm in 2016 and higher and unstable
on 2016 Jul 19. The lower panel of Figure 1 shows the resulting 230 GHz SMA
light curve.

###  2.5 ALMA Data

This analysis used 12 epochs (1374 minutes) of data from the ALMA observatory.
The first two epochs (project code 2015.A.00021.S, PI G. Witzel) were taken in
2016 while Spitzer was observing at 4.5 µm. Details are in Table 1. The data
were calibrated with the standard ALMA pipeline using the Common Astronomy
Software Application package (CASA: McMullin et al., 2007). A few spectral
windows showed suspicious absorption features and were excluded from the
analysis. Because ALMA has relatively fewer short baselines than SMA, all
baselines were included, and the source was imaged with uniform weighting
using task clean. Sgr A* is a strong source at mm wavelengths, and its
visibility was mostly flat as a function of baseline length. That means the
source is compact, and with ALMA’s good U–V coverage and uniform weighting,
light curves could be extracted by simply measuring the peak flux densities of
the image. The off-source RMS gave the uncertainty. The robustness of the
results was investigated using 1–30-minute averaging times. The final choice
of 3 minutes time is a good balance between signal-to-noise and sampling the
light curve. We obtained two long light curves, the first 389 minutes with a
gap and the second 280 minutes as shown in Table 1.

The remaining 10 ALMA epochs are 70-minute light curves at 234 GHz (Iwata et
al., 2020). The epochs were observed in 2017 October over ten days, i.e., they
are separated by roughly one day and have average cadence of ${\sim}1.6$
minutes. Iwata et al. (2020) gave details of the data reduction. The lower
panel of Figure 1 shows all 12 epochs of 230 GHz ALMA data.

###  2.6 Simultaneous NIR and Submm Observations

A majority of the submm data were taken during the Spitzer observations,
resulting in the largest set of synchronous submm (230 or 345 GHz) and NIR
(4.5 µm) light curves. (Chandra observed simultaneously as well, but there was
no X-ray event during times of SMA or ALMA coverage.) In total, we obtained
$>$1800 minutes on-source within a total duration of 2247 minutes of
simultaneous data. The combined ALMA and SMA data on 2016 Jun 12 cover 787
minutes with a gap of 116 minutes between the datasets, making this the
longest simultaneous light curve so far. Figure 2 shows all simultaneous data
taken during the Spitzer campaigns.

Figure 3 shows the discrete cross-correlation (DCC) of the entire 230 GHz
dataset with the 4.5 µm light curves. Because the data were not regularly
sampled and show large gaps, we used the algorithm of Edelson & Krolik (1988)
as implemented by Robertson et al. (2015). We calculated the DCC for both the
data as presented in Figure 2 and for a transformed version of the 4.5 µm
data. The transform was to convolve the data with a Hanning smoothing kernel
with size of 61 minutes and then take the logarithm. (Section 6 explains the
reasons for doing this.) Both DCCs have their maxima at time lags of 20–30
minutes, positive lag meaning that the 4.5 µm data lead the submm.

Figure 3: Discrete cross-correlation (DCC) functions of the 230 GHz and 4.5 µm
light curves shown in Figure 2. Orange shows the DCC of the unaltered data,
and blue with shaded error bars shows the DCC of 230 GHz with the logarithm of
the low-pass-filtered 4.5 µm data. The vertical dashed line indicates a time
lag of zero. The DCC reaches its maximum at a time lag of about 33 minutes for
the unaltered data and at $27\pm 8$ minutes (solid vertical line) with the
transformed data. Positive lags indicate the 4.5 µm data lead the submm. The
quoted uncertainty is that of the peak position of the part of the DCC that
lies above the 95% false alarm probability (dashed horizontal curves), which
was derived from our Bayesian posterior (Section 2).

##  3 The X-ray Power Spectrum and the NIR–X-ray Correlation

The first step of our analysis was to determine the PSD of the X-ray
variability of Sgr A* under the preliminary assumption of a log-normal PDF for
the X-ray flux densities. We did this by forward-modeling light curves for a
range of PSDs, comparing the model light curves with observed ones through a
suitable summary statistic (a “distance function”), and exploring the
parameter space by approximate Bayesian computation (ABC). ABC is appropriate
for problems with no analytic likelihood function. All three parts—the forward
modeling, distance function, and implementation of the ABC—were described in
detail by Witzel et al. (2018, their Appendix B). The present analysis closely
follows their procedure but was adapted for X-rays by deriving photon
statistics from the modeled variable flux densities plus a constant source
contribution. Details are given below.

A log-normal PSD was chosen because a) it can—as the results show—describe the
observed data successfully, and b) a log-normal was the preferred model for
NIR light curves (Witzel et al., 2018). That analysis showed that the NIR PSD
is a single broken power law with a slope of 2 (red noise) and a precisely
measured break timescale $\tau_{b}=243^{+82}_{-57}$ minutes. In all such
analyses, the inferred PSD parameters depend on the model for the flux-density
PDF, and using a log-normal PSD lets us directly compare X-ray PSD parameters
with the NIR. The analysis in Section 4 replaces the PSD model assumption with
a physical model.

###  3.1 Forward Modeling of Light Curves

We used the FFT-based method (Timmer & Koenig, 1995) to generate random light
curves from a given PSD and a set of independently drawn random numbers. This
method results in realizations of a Gaussian process $g(t)$ exhibiting
periodograms (i.e., PSD estimators) consistent with the input PSD. We
parameterized the PSD as a broken power law of the form:

$\mathrm{PSD}(f)\propto\left\\{\begin{array}[]{lcl}f^{-\gamma_{0}}&\mathrm{for}&f<f_{b}\\\
f^{-\gamma}&\mathrm{for}&f\geq f_{b}\quad,\end{array}\right.$ (1)

where $f$ is temporal frequency, and we assume $\gamma_{0}=0$ (Meyer et al.
2009). In order to generate realistic light curves, we transformed the values
$g(t)$ to make the resulting distribution of $S(t)=T[g(t)]$ consistent with
the distribution of observed flux densities. With a log-normal distribution as
our target,

$T(g(t))=\exp{(\sigma_{\rm{logn}}\cdot g(t)+\mu_{\rm{logn}})}\;\;,$ (2)

$\mu_{\rm{logn}}$ and $\sigma_{\rm{logn}}$ being the log-normal parameters.

For Chandra data, we modeled the light curves in count rate
$\Lambda(t)=T[g(t)]$ (counts per second, cps) instead of flux density $S(t)$.
The conversion factor between flux density (as observed at Earth, i.e., after
suffering interstellar extinction) and count rate is absorbed by the log-
normal mean $\mu$ and otherwise does not affect log-normality. We assumed an
effective area $e_{G}=1$ for the ACIS-S/HETG data and Gaussian priors for
$e_{I}/e_{G}$ and $e_{S}/e_{G}$, the relative effective areas of the other two
modes.222Our priors on $e_{I}/e_{G}$ and $e_{S}/e_{G}$ are based on the
spectral model used by Zhu et al. (2018), i.e., an absorbed bremsstrahlung
with a plasma temperature of 10 keV and a foreground absorption column density
of $10^{23}$ cm-2. For each mode, we included an additive term
($\chi_{G}$,$\chi_{S}$,$\chi_{\rm{I}}$) with flat priors to account for the
contribution from the extended, non-varying X-ray source at the position of
Sgr A*. The resulting total count rate (before measurement errors) from the
position of Sgr A* is then:

$\Lambda_{\rm{tot}}(t)=e_{i}/e_{G}\cdot\left[\Lambda(t)+\chi_{i}\right]\;\;,$
(3)

where $i$ is $S$, $I$, or $G$ for the three Chandra instruments used, and
$\chi_{i}$ is the count rate of the constant source as seen by the respective
instrument. Pileup was included in the forward modeling for ACIS-I and
ACIS-S/HETG data for all count rates $>$0.02 cps. (It is negligible below
that.) We used Equations 3 and 4 of Yuan & Wang (2016):

$\displaystyle\Lambda_{\rm{out},I}(t)=$
$\displaystyle\Big{[}4.180\Lambda_{\rm{tot}}(t)^{−0.07387}$
$\displaystyle+0.5381\Lambda_{\rm{tot}}(t)^{−1.160}\Big{]}^{-1}\;\;,$ (4)

and

$\displaystyle\Lambda_{\rm{out},G}(t)=$
$\displaystyle\Big{[}3.933\Lambda_{\rm{tot}}(t)^{−0.03541}$
$\displaystyle+0.6564\Lambda_{\rm{tot}}(t)^{−1.107}\Big{]}^{-1}\;\;.$ (5)

For ACIS-I count rates between 0.1 and 0.2 cps, the corrections are 27 to 64%.
We finally modeled the measured count rate $\Lambda_{\rm{meas}}$, sampled once
per minute, by a Poisson process:

$\Lambda_{\rm{meas}}(t)=\rm{Pois}\left[\Lambda_{\rm{out}}(t)\cdot
60\right]/60\;\;.$ (6)

With this 1-minute binning, the effect of the Chandra frame time (3.2 s for
ACIS-I and HETG, 0.4 s for ACIS-S subarray) is negligible.

###  3.2 The Distance Function

Following Witzel et al. (2018, their Appendix B.2), we used the first order
structure function as the distance function. The structure function quantifies
the variance of the flux density at any given time separation and contains
information on the PSD as well as the flux-density distribution of the
variability process. The structure function is defined as:

$\displaystyle V(\tau_{i})$
$\displaystyle=\frac{1}{n_{i}}\sum_{t_{j},t_{k}}\left[F(t_{j})-F(t_{k})\right]^{2}$
(7)
$\displaystyle\quad\mathrm{for}\;\tau_{i}\leq(t_{j}-t_{k})<\tau_{i+1}\quad,$

that is, the sum of $\left[F(t_{j})-F(t_{k})\right]$ over all measurement
pairs whose time lags $(t_{j}-t_{k})$ fall within the bin
$[\tau_{i},\tau_{i+1}]$, there being $n_{i}$ such pairs. The structure
functions for the three modes of Chandra data are shown in Figure 4.

We defined the distance between two light curves as the weighted L2 norm of
the difference between the logarithms of the respective structure function’s
binned values:

$\phi(V_{1},V_{2})=\sum_{i}w_{i}(\log\left[V_{1}(\tau_{i})/V_{2}(\tau_{i})\right])^{2}\quad$
(8)

with $w_{i}$ the weights for the chosen binning. The weights adopted here were
unity for each structure function bin except the single wide bin at large time
lags, which had $w_{i}=3$. These values gave uniform and reasonably quick
convergence of the fit.

Figure 4: Logarithmically binned structure functions (Eq. 7). Upper panel
shows the X-ray SFs, red for ACIS-I, green for ACIS-S, and blue fpr
ACIS-S/HETG. Lower panel shows the 230 GHz and 345 GHz SFs, red the combined
SF from SMA and ALMA data, and blue the APEX data combined with one epoch of
SMA data.

There is a general problem of consistency when comparing variability
timescales quoted in the literature. The definition of “timescale” can depend
on the method used for its estimation, and ignoring that can lead to errors in
the interpretation of time series analysis results, especially in their
conversion into physical quantities. The most substantial difference in the
definition of timescales concerns the domain in which they were estimated. In
the spectral domain, a timescale naturally refers to the period of a
sinusoidal signal. In the time domain, in particular when the estimation is
through a structure function, “timescale” refers to a characteristic interval
of coherent variation. (Roughly speaking, this is the average time of
monotonic increase or decrease of the analyzed quantity.) The
Ornstein–Uhlenbeck-process (OU-process) timescale given by Dexter et al.
(2014) is in the time domain and is consistent with the definition provided by
the structure function. For a sinusoidal function of frequency $f$, the time-
domain timescale equals $0.5/f$. The equation
$\tau_{\mathrm{time\;domain}}=0.5/f_{\mathrm{spectral\;domain}}$ is in fact
valid for any function because, through a Fourier analysis, we can describe
any function as a superposition of sinusoids. As long as the consistency issue
is taken into account, different definitions of timescales are equivalent. The
relevant question is which timescale is more informative about physical
properties. The “memory” or coherence time of a physical process can be better
described by a structure-function timescale rather than by a spectral-domain
one. Structure-function timescales are also the appropriate choice for the
calculation of source sizes through causality arguments. Periodicities, on the
other hand, are naturally quantified in the spectral domain. This paper uses
the structure function to define a distance between mock data and observed
data. However, our timescale parameter $f_{b}$ is defined in Equation 1 and
refers to the spectral domain.

Caution should also be exercised in comparing structure-function slopes. While
originally defined as in Equation 7, the structure function is sometimes
expressed as

$\displaystyle V(\tau_{i})$
$\displaystyle=\sqrt{\frac{1}{n_{i}}\sum_{t_{j},t_{k}}\left[F(t_{j})-F(t_{k})\right]^{2}}$
(9)
$\displaystyle\quad\mathrm{for}\quad\tau_{i}\leq(t_{j}-t_{k})<\tau_{i+1}\quad,$

(e.g., Dexter et al. 2014). Conversion of the structure function slope into a
PSD power-law index depends on the structure function definition used.

###  3.3 Approximate Bayesian Computation

ABC is an iterative method based on prior distributions of model parameters,
two Monte Carlo sampling steps—one for picking a random parameter set and one
for drawing a light curve realization for this parameter set—and an acceptance
step. The acceptance is based on the comparison of the distance between
simulated and observed data with a threshold decreasing with each iteration.
For each iteration, the two sampling steps are executed many times until $n$
(e.g., $n=500$) realizations are accepted. The parameters (called “particles”)
of the accepted realizations represent the current best estimate of the
posterior. This estimate is used to inform the parameter sampling of the next
iteration. If the distance function is informative for the model at hand, the
prior distribution is transformed into a close approximation of the posterior
after a sufficient number of iterations. A detailed description and
explanation of the algorithm and its implementation is given by Witzel et al.
(2018, their Appendix B.4).

Table 2: Priors and Posteriors of Bayesian Analysis |  | Median of |
---|---|---|---
Parameter | Prior | Posterior | Description
Analysis of X-ray PSD, log-normal model
$\gamma$ | flat on [0.2, 5.5] | $\phantom{-}3.70^{+0.75}_{-0.70}$ | PSD slope
$f_{b}$ [$10^{-3}\cdot\rm{minutes}^{-1}$] | flat on [1.0, 600.0] | $\phantom{-}3.4^{+1.5}_{-1.2}$ | correlation frequency
$\mu_{\rm{logn}}$aa$\mu_{\rm{logn}}$, $\sigma_{\rm{logn}}$, and $\chi$ are data descriptive quantities representing the observed count rate in cps at the Chandra detector, i.e., after interstellar extinction (Equations 2 and 3). | flat on [$-$35.0,0.0] | $-15.6^{+3.0}_{-2.9}$ | log-normal mean
$\sigma_{\rm{logn}}$aa$\mu_{\rm{logn}}$, $\sigma_{\rm{logn}}$, and $\chi$ are data descriptive quantities representing the observed count rate in cps at the Chandra detector, i.e., after interstellar extinction (Equations 2 and 3). | flat on [0.0001, 6.0] | $\phantom{-}3.59^{+0.97}_{-0.95}$ | log-normal standard deviation
$\chi_{\mathrm{S}}$ [$10^{-3}\cdot$ cps]aa$\mu_{\rm{logn}}$, $\sigma_{\rm{logn}}$, and $\chi$ are data descriptive quantities representing the observed count rate in cps at the Chandra detector, i.e., after interstellar extinction (Equations 2 and 3). | flat on [0.5,4.5] | $\phantom{-}2.8\pm{0.2}$ | contribution of steady X-ray source, ACIS-S
$\chi_{\mathrm{I}}$ [$10^{-3}\cdot$ cps]aa$\mu_{\rm{logn}}$, $\sigma_{\rm{logn}}$, and $\chi$ are data descriptive quantities representing the observed count rate in cps at the Chandra detector, i.e., after interstellar extinction (Equations 2 and 3). | flat on [0.5,4.5] | $\phantom{-}1.9\pm{0.2}$ | contribution of steady X-ray source, ACIS-I
$\chi_{\mathrm{G}}$ [$10^{-3}\cdot$ cps]aa$\mu_{\rm{logn}}$, $\sigma_{\rm{logn}}$, and $\chi$ are data descriptive quantities representing the observed count rate in cps at the Chandra detector, i.e., after interstellar extinction (Equations 2 and 3). | flat on [0.5,4.5] | $\phantom{-}2.2\pm{0.1}$ | contribution of steady X-ray source, ACIS-S/HETG
$e_{\rm{I}}/e_{G}$ | Gaussian ($\mu=2.70$, $\sigma=0.05$) | $\phantom{-}2.70\pm{0.04}$ | effective area of ACIS-I rel. to ACIS-S/HETG
$e_{S}/e_{G}$ | Gaussian ($\mu=3.0$, $\sigma=0.1$) | $\phantom{-}3.02\pm{0.08}$ | effective area of ACIS-S rel. to ACIS-S/HETG
Analysis of X-ray, NIR, and submm data, synchrotron–SSC model
$\gamma_{e}$ | flat on [1.5, 4.0] | $\phantom{-}2.95^{+0.21}_{-0.19}$ | power-law index of electron energy distribution
$B_{0}$ [G] | log flat on [$10^{-2}$, $10^{3}$] | $\phantom{-}8.5^{+1.8}_{-1.4}$ | minimum magnetic flux density
$L_{0}$ [$R_{S}$] | log flat on [$10^{-2}$, $10^{3}$] | $\phantom{-}2.7\pm{0.3}$ | maximum physical source diameter
$\gamma_{\rm{slow}}$ | flat on [1.2, 6.5] | $\phantom{-}4.9\pm{0.8}$ | PSD slope of slow Gaussian process
$\gamma_{\rm{fast}}$ | flat on [1.2, 9.5] | $\phantom{-}2.1^{+0.4}_{-0.3}$ | PSD slope of fast Gaussian process
$f_{b,\rm{slow}}$ [$10^{-2}\cdot\rm{minutes}^{-1}$] | flat on [0.001, 2.0] | $\phantom{-}0.74^{+0.13}_{-0.12}$ | correlation frequency of slow Gaussian process
$f_{b,\rm{fast}}$ [$10^{-2}\cdot\rm{minutes}^{-1}$] | flat on [0.001, 2.0] | $\phantom{-}1.21^{+0.29}_{-0.27}$ | correlation frequency of fast Gaussian process
$\mu_{\rm slow}$bb$\mu$ and $\sigma$ are the log-normal parameters in Equations 38 and 39. They describe the time series of optically thin synchrotron flux density $S_{\rm{thin}}(t)$ in Jy at $\nu_{\rm{NIR}}=136269$ GHz, and the exponential cutoff frequency $\nu_{2}(t)$ in units of $\nu_{\rm{NIR}}$. | flat on [0,10] | $-5.77^{+0.19}_{-0.26}$ | log-normal mean of slow process
$\sigma_{\rm slow}$bb$\mu$ and $\sigma$ are the log-normal parameters in Equations 38 and 39. They describe the time series of optically thin synchrotron flux density $S_{\rm{thin}}(t)$ in Jy at $\nu_{\rm{NIR}}=136269$ GHz, and the exponential cutoff frequency $\nu_{2}(t)$ in units of $\nu_{\rm{NIR}}$. | flat on [0.001, 1.5] | $\phantom{-}0.616^{+0.075}_{-0.059}$ | log-normal standard deviation of slow process
$\mu_{\rm fast}$bb$\mu$ and $\sigma$ are the log-normal parameters in Equations 38 and 39. They describe the time series of optically thin synchrotron flux density $S_{\rm{thin}}(t)$ in Jy at $\nu_{\rm{NIR}}=136269$ GHz, and the exponential cutoff frequency $\nu_{2}(t)$ in units of $\nu_{\rm{NIR}}$. | flat on [$-$7.5, 4.0] | $-1.65^{+0.66}_{-0.57}$ | log-normal mean of fast process
$\sigma_{\rm fast}$bb$\mu$ and $\sigma$ are the log-normal parameters in Equations 38 and 39. They describe the time series of optically thin synchrotron flux density $S_{\rm{thin}}(t)$ in Jy at $\nu_{\rm{NIR}}=136269$ GHz, and the exponential cutoff frequency $\nu_{2}(t)$ in units of $\nu_{\rm{NIR}}$. | flat on [0.001, 3.0] | $\phantom{-}1.55^{+0.46}_{-0.40}$ | log-normal standard deviation of fast process
$\chi_{I}$ [$10^{-3}\cdot$ cps]aa$\mu_{\rm{logn}}$, $\sigma_{\rm{logn}}$, and $\chi$ are data descriptive quantities representing the observed count rate in cps at the Chandra detector, i.e., after interstellar extinction (Equations 2 and 3). | flat on [0.0,0.005] | $\phantom{-}1.1^{+0.8}_{-0.5}$ | contribution of steady, extended X-ray source
$\sigma_{\mathrm{Keck}}$ [mJy] | Gaussian ($\mu=0.015$, $\sigma=0.008$) | $\phantom{-}0.018\pm{0.006}$ | measurement noise of the Keck observations
$\sigma_{\mathrm{VLT}}$ [mJy] | Gaussian ($\mu=0.031$, $\sigma=0.008$) | $\phantom{-}0.034\pm{0.005}$ | measurement noise of the VLT observations
$\sigma_{\mathrm{IRAC}}$ [mJy] | Gaussian ($\mu=0.65$, $\sigma=0.4$) | $\phantom{-}0.654^{+0.041}_{-0.047}$ | measurement noise of the IRAC observations
$\sigma_{345\,\mathrm{GHz}}$ [Jy] | Gaussian ($\mu=0.1$, $\sigma=0.1$) | $\phantom{-}0.10\pm{0.01}$ | measurement noise of the 345 GHz (APEX and SMA)
$\sigma_{230\,\mathrm{GHz}}$ [Jy] | Gaussian ($\mu=0.15$, $\sigma=0.1$) | $\phantom{-}0.08\pm{0.01}$ | measurement noise of the 230 GHz (SMA and ALMA)
$e_{S}/e_{I}$ | Gaussian ($\mu=1.12$, $\sigma=0.04$) | $\phantom{-}1.12\pm{0.03}$ | effective area of ACIS-S relative to ACIS-I
$e_{G}/e_{I}$ | Gaussian ($\mu=0.38$, $\sigma=0.02$) | $\phantom{-}0.38\pm{0.01}$ | effective area of ACIS-S/HETG relative to ACIS-I
$\Delta A_{{K}}$ | Gaussian ($\mu=0.0$, $\sigma=0.1$) | $-0.01\pm{0.07}$ | modification of $K$-band extinction
$\Delta A_{{M}}$ | Gaussian ($\mu=0.0$, $\sigma=0.14$) | $\phantom{-}0.09^{+0.08}_{-0.07}$ | modification of $M$-band extinction
Posteriors (medians and 16% and 84% quantiles) for the synchrotron-SSC model
derived from 1000 light curves of 700 minutes duration each
$B(t)$ [G] |  | $\phantom{-}12.7^{+4.1}_{-2.9}$ | time-dependent magnetic flux density
$L(t)$ [$R_{S}$] |  | $\phantom{-}2.20^{+0.32}_{-0.30}$ | time-dependent source size
$n_{e}(t)$ [$10^{7}$ cm-3] |  | $\phantom{-}4.3^{+6.9}_{-2.7}$ | time-dependent electron density
$\alpha_{\rm{NIR}}(t)$ |  | $-1.65^{+0.29}_{-0.30}$ | time-dependent NIR spectral index (2.2 to 4.5 $\micron$)
$\nu_{m}(t)$ [GHz] |  | $\phantom{-}308^{+96}_{-75}$ | time-dependent self-absorption turnover
$\nu_{2}(t)$ [THz] |  | $\phantom{-}54^{+111}_{-29}$ | time-dependent cooling cutoff frequency
Figure 5: Results for the Bayesian structure function fit for the X-ray
assuming a log-normal flux-density probability density function. Contours show
the joint (posterior) probability density for each parameter pair, and panels
along the right upper edge show histograms of the marginalized posterior of
each parameter defined in Table 2. For each histogram, the dashed lines mark
the 16%, 50%, and 84% quantiles.

Our custom C++ implementation of the forward modeling and ABC algorithm
includes a fast algorithm for repetitive calculations of structure functions
(Witzel et al. 2018, Appendix C) and is MPI-integrated333MPI means “message-
passing interface,” a standard for parallel computing. to be run on computing
clusters. For this analysis, we modified the model to include photon
statistics according to Equations 3 to 6. We ran the ABC on 200 cores of the
VLBI correlator computer cluster of the Max Planck Institute for Radio
Astronomy (Bonn) with 20 000 initial light curve drawings. From these we
selected the particles having the smallest distance values. The final ABC run
is the result of 30 iterations with $n=10\,000$ particles.

###  3.4 PSD Parameters of the X-ray Variability

Model parameters and their priors and posteriors are listed in Table 2 and
shown in Figure 5. While the posteriors of $e_{I}/e_{G}$, and $e_{S}/e_{G}$
are mere minor alterations of the priors, the posteriors of all other
parameters are constrained, peaked distributions. The log-normal parameters
$\mu$ and $\sigma$ are highly correlated. A more precise determination of
these parameters is not possible because of photon noise and pileup, which
make the log-normal distributions for some $\mu$ and $\sigma$ combinations
indistinguishable. The background levels $\chi_{i}\cdot e_{i}/e_{G}$ are
reasonable when compared to X-ray spectra accumulated over all periods of time
without flares. However, $\chi_{S}\cdot e_{S}/e_{G}$ is somewhat higher, which
likely can be attributed to the contribution of the magnetar.

The parameters of main interest are $\gamma$ and $f_{b}$, shown in Figure 6 in
more detail. Determinations of both parameters are about an order of magnitude
less precise for the Chandra data than for the NIR. This is not surprising
considering that in 18 years only $\sim$34 000 X-ray photons were detected.
However, despite the lack of precision, the comparison with the NIR is
fruitful: the X-ray contours are displaced towards higher slopes and slightly
smaller break frequencies. While the two sets of contours can be interpreted
as marginally consistent, Figure 6 suggests that the X-ray variability does
not show as much power at high frequencies as the NIR. This interpretation of
suppressed power of the fastest variability turns out to be a clue to
understanding the NIR–X-ray correlations and allows us to identify a radiative
model that can explain many aspects of the rapid (i.e., minutes to hours)
variability of Sgr A* from the submm to the X-rays.

Figure 6: Credible contours (68%, 95%, 99%) for the parameters $\gamma$ and
$f_{b}$. Blue shows the results of the present analysis for X-ray data (Table
2), and red shows results of the NIR analysis of (Witzel et al., 2018, their
Case 3). In both analyses, the posteriors have been marginalized over all
other parameters.

###  3.5 Simulations of NIR–X-ray Correlations

Figure 7 demonstrates how PSD parameters affect NIR–X-ray correlations The
figure shows two light curves from the same set of Gaussian random numbers
with two different PSDs. In the previously mentioned method of Timmer & Koenig
(1995), a random number—one for each relevant frequency in the Fourier
representation of the light curve—is multiplied by the square root of the PSD
at the corresponding frequency. By, let’s say, using a broken power-law PSD
model for one of the light curves and a second PSD with the same break
frequency but a steeper slope for the other, we can generate pairs of light
curves. Figure 8a) shows a flow chart for the procedure. Because the curves
use the same random numbers, they are strongly correlated. However, one of
them is the low-pass filtered version of the other. When we apply two, e.g.,
log-normal transformations (with the one applied to the slower process showing
a heavier tail), the slower, more non-linear process (representative of the
X-ray variability) will always have a counterpart in the faster process
(representative of the NIR), while the faster process shows a lot of peaks
having a wide range of levels with (almost) no counterpart. When peaks are
seen in both curves, their arrival times are strictly correlated, but there is
little correlation in peak flux. This phenomenology mimics the NIR–X-ray
observations (Genzel et al., 2010; Fazio et al., 2018) and becomes even more
realistic when we include photon statistics, pileup effect, and a constant
Poissonian background for the slower process.

Figure 7: Multiwavelength mock light curves generated from the X-ray data
model. The lower panel shows a 50 000-minute realization of log-normal NIR
flux densities, and the upper panel presents X-ray count rates during the same
time interval. Count rates are as observed and include the steady X-ray
background, pileup, and Poisson noise. (Flux-density normalization and the
extinction correction are implicit in the model parameter $\mu$.) Figure 8a
illustrates the procedure for generating the mock light curves.

Figure 8: Flow chart of simulation algorithms for modeling light curves. Upper
panel a) shows an example of low-pass correlated Gaussian processes under
transformation with log-normal distributions. The generation of each Gaussian
process from random numbers follows the procedure described by Timmer & Koenig
(1995). The two PSDs applied are both broken power laws, one of which has—by
means of a steeper slope and/or a smaller break frequency—significantly lower
power at short variability timescales. Lower panel b) shows a similar low-pass
correlated set of log-normal trigger processes that are then used to inform
our synchrotron–SSC model. The radiative model combines the fast and slow
processes according to the equations developed in Section 4.4 and lets us
derive time series of flux densities at any given wavelength.

##  4 Synchrotron Self-Compton Scattering and the Evolution of the Optical
Depth

The theory of astronomical synchrotron sources was developed in the 1950s and
1960s (Ginzburg 1951; Shklovskii 1952, 1953; Ginzburg & Syrovatskii 1965).
This paper uses the later, systematic development of synchrotron and
synchrotron-SSC theory by Kardashev (1962), van der Laan (1966), Gould (1979),
Marscher (1983), and Band & Grindlay (1985). In particular, we assume
familiarity with the excellent overview article by Moffet (1975).

Based on flux-density levels, timing properties, and correlations, SSC
scattering is a good candidate to explain Sgr A*’s X-ray variability. SSC
becomes relevant under the same circumstances—high luminosity in a source of
small size (Moffet, 1975)—as synchrotron self-absorption, and therefore we
discuss the co-evolution of SSC scattering and synchrotron opacity. A simple
model of injection, compression, and expansion can account for many aspects of
the time-variable SED of Sgr A*. A semi-analytic approach allows us to co-fit
the structure functions at all wavelengths considered in this study.

###  4.1 The NIR–X-ray Correlation

As mentioned before, several studies have found the flux-density peak levels
of synchronous NIR and X-ray flares to be consistent with the prediction of
simple SSC models. Furthermore, Witzel et al. (2018) showed evidence that the
NIR variability is at least partially caused by a variable, exponential
synchrotron cooling cutoff close to the NIR; i.e., the synchrotron spectrum of
flux densities is of the form:

$S_{\nu}(t)=I_{\nu}\exp\left[-\left(\frac{\nu}{\nu_{2}(t)}\right)^{\frac{1}{2}}\right]\Delta\Omega\;\;,$
(10)

where the source subtends a solid angle of $\Delta\Omega=\pi(L/D_{s})^{2}$
with $L$ the source diameter and $D_{s}$ the distance to the source. $I_{\nu}$
is the synchrotron intensity at frequency $\nu$, and $\nu_{2}(t)$ is the
variable cutoff frequency due to synchrotron cooling. Building on these ideas,
the analytical equation for SSC radiation (Marscher, 1983) also depends on
$\nu_{2}(t)$ but very weakly because it is the argument of a logarithmic term:

$S(E_{\rm{keV}})=d(-\alpha)\ln{(\nu_{2}/\nu_{m})}\theta^{(4\alpha-6)}\nu_{m}^{(3\alpha-5)}S_{m}^{(4-2\alpha)}E_{\rm{keV}}^{\alpha}\;\;,$
(11)

where $\alpha$ is the optically thin spectral index,444We define $\alpha$ such
that $S_{\nu}\propto\nu^{\alpha}$. Marscher (1983) used the opposite sign
convention., $d(-\alpha)$ are coefficients given by Marscher (1983),
$\theta=L/D_{s}$ is the angular size of the source, $\nu_{m}$ is the self-
absorption turnover frequency, $S_{m}$ is the self-absorption flux density,
and $E=h\nu$ is the photon energy in keV with $h\approx 4.136\times 10^{-9}$
keV/GHz the Planck constant. Equation 11 is valid for

$5.5\times 10^{-9}\gamma_{1}^{2}\nu_{m}\lesssim E_{\rm{keV}}\lesssim
0.2[b(-\alpha)]^{-1}\theta^{-4}\nu_{m}^{-5}S_{m}^{2}$ (12)

with $b(-\alpha)$ coefficients given by Marscher (1983) and $\gamma_{1}$ the
minimum Lorentz factor of the electrons. For energies in this range, the SSC
flux density is a power law with the same spectral index as the NIR. Figure 9
shows a synchrotron–SSC SED for realistic Sgr A* parameters. (Section 4.2
explains how to calculate $I_{\nu}$ from these parameters.)

The logarithmic dependence of $S(E_{\rm{keV}})$ on the cutoff $\nu_{2}(t)$
implies that if $\nu_{2}$ is the origin of the fast variability in the NIR,
this variability power is suppressed in the X-ray light curves. However, it is
not obvious that Equation 11 indeed suggests a clear separation of variability
power: the X-ray flux density is highly variable, and this variability must be
related to changes of the synchrotron source and spectrum itself, i.e., to
changes of $\nu_{m}$, $S_{m}$, and $\theta$. This in turn means that $I_{\nu}$
and $\Delta\Omega$ in Equation 10 are time-dependent with similarly slow
variability as the SSC flux density. Thus, following the idea presented in
Section 3.5 and if SSC is responsible for the X-ray emission, the NIR flux
density $S_{\nu_{\rm{NIR}}}(t)$ is the product of two correlated processes,
$I_{\nu_{\rm{NIR}}}(t)\Delta\Omega(t)$ and
$\exp{\left[-\left(\frac{\nu_{\rm{NIR}}}{\nu_{2}(t)}\right)^{{1}/{2}}\right]}$,
the former being the low-pass-filtered counterpart of the latter. This
situation—the slow process (postulated from the X-rays) feeding back into the
fast process (the NIR variability, described by the product)—is slightly more
complex than the simple case considered in Section 3.5. However, for the right
model for $I_{\nu}(t)$ and proper model parameters, the fast process will
dominate this product, and the result will indeed be similar to the light
curves of Figure 7.

Figure 9: Model radio-to-X-ray SED for Sgr A*. Orange shows the synchrotron
component, blue the SSC component, and red their sum. The SED is shown as if
interstellar extinction were zero. Model parameters (see Table 2 for
definitions) are $n_{e}=1.066\times 10^{8}$ cm-3, $B=10$ G, $L_{0}=2$ $R_{S}$,
$\nu_{2}=20$ THz, $E_{\rm{min}}=0.02$ GeV, and $\gamma=2.4$. The vertical
dashed lines mark the frequencies 230 GHz, 345 GHz, and $\nu_{\rm{NIR}}$. The
blue rectangle marks the 2–8 keV X-ray band and shows typical peak flux
densities during X-ray flares. The increasing and decreasing flanks of the SSC
spectrum (where the spectrum is not described by a simple power law according
to Equation 11) are illustrated by arbitrary exponential cutoffs. The precise
shape of the SSC spectrum in these flanks depends on many assumptions
including the density profile of the source and is beyond the scope of this
study. The spectrum shown is not a fit to any actual data. In particular it
combines low-level NIR flux densities with detectable X-ray emission, a
condition that has not been observed.

According to Equation 11, $S(E_{\rm{keV}})$ is not a direct function of the
optically thin flux density at the $K$-band frequency $\nu_{\rm{NIR}}=136269$
GHz,

$S_{\rm{thin}}=I_{\nu_{\rm{NIR}}}(t)\Delta\Omega(t)\;\;;$ (13)

the self-absorption turnover flux density and turnover frequency also matter.
In order to generate light curves in the submm that show the observed
correlations and delays relative to the NIR, we need to discuss possible
scenarios of co-evolution of $S_{\rm{thin}}$ and the synchrotron opacity. In
particular, it is not correct to assume $S_{m}\propto S_{\rm{thin}}$ in the
context of Equation 11 and to treat $\nu_{m}$ and $\theta$ as independent
parameters (as often done in the literature) and then to argue that
$S(E_{\rm{keV}})$ depends on the optically thin flux density in a highly non-
linear way due to the term $S_{m}^{(4-2\alpha)}$.

###  4.2 Evolution of the Synchrotron Source and Synchrotron Self-absorption

Following the formalism and notation of Moffet (1975), the solution of the
radiative transfer equation
$dI_{\nu}=(\epsilon_{\nu}-\kappa_{\nu}I_{\nu})\;dz$ through a homogeneous slab
of material with constant and isotropic emissivity $\epsilon_{\nu}$,
absorption $\kappa_{\nu}$, and thickness $z_{t}$ is

$I_{\nu}(z_{t})=\frac{\epsilon_{\nu}}{\kappa_{\nu}}[1-\exp{(-\kappa_{\nu}z_{t})}]~{}~{}.$
(14)

For a synchrotron source of homogeneous electron density in a tangled magnetic
field, assuming a power-law electron energy distribution and an isotropic
pitch angle distribution,

$\displaystyle\epsilon_{\nu}$
$\displaystyle=\frac{1}{2}C_{2}n_{0}B^{(\gamma_{e}+1)/2}(C_{1}/\nu)^{(\gamma_{e}-1)/2}$
$\displaystyle~{}~{}\cdot
G^{\prime}(\nu/\nu_{1},\nu/\nu_{2},\gamma_{e})~{}~{},$ (15)

and

$\displaystyle\kappa_{\nu}$
$\displaystyle=\frac{1}{2}c^{2}C_{2}C_{1}^{\gamma_{e}/2}n_{0}B^{(\gamma_{e}+2)/2}\nu^{-(\gamma_{e}+4)/2}(\gamma_{e}+2)$
$\displaystyle~{}~{}\cdot
G^{\prime}(\nu/\nu_{1},\nu/\nu_{2},\gamma_{e}+1)~{}~{}.$ (16)

Here $c$ is the speed of light, $C_{1}$ and $C_{2}$ are constants given in
Table 3, $\gamma_{e}=1-2\alpha$ is the power-law index of the electron energy
distribution, $n_{0}$ is the scaling factor of the electron energy
distribution, and $B$ is the magnetic flux density in gauss.
$G^{\prime}(\nu/\nu_{1},\nu/\nu_{2},\gamma_{e})$ is defined as

$\displaystyle G^{\prime}(\nu/\nu_{1},\nu/\nu_{2},\gamma_{e})$
$\displaystyle\equiv\frac{1}{2}\pi^{1/2}\frac{\Gamma[(\gamma_{e}+5)/4]}{\Gamma[(\gamma_{e}+7)/4]}$
$\displaystyle~{}~{}\cdot\int_{0}^{\infty}x^{(\gamma_{e}-3)/2}F(x)dx\;\;,$
(17)

with

$F(x)=x\int_{x}^{\infty}K_{5/3}(\eta)d\eta\;\;,$ (18)

$K_{5/3}(\eta)$ the modified Bessel function of the second kind, and
$\Gamma(x)$ the Gamma function. Here we set $\nu/\nu_{1}=\infty$ and
$\nu/\nu_{2}=0$. However, this does not imply that the maximal Lorentz factor
$\gamma_{\rm{max}}$ of the electrons is infinite. Equation 10 includes an
exponential term to approximate the effect of a finite $\gamma_{\rm{max}}$.
This approach allows us later to use a lookup table for
$G^{\prime}(\infty,0,\gamma_{e})$ without having to re-evaluate the integral
for every frequency.

A main result of synchrotron theory is that the SED is peaked and
asymptotically approaches two different power laws: the optically thick
spectrum $I_{\nu}=\epsilon_{\nu}/\kappa_{\nu}$ with slope $+2.5$ on the low-
frequency side of the peak and the optically thin spectrum
$I_{\nu}=\epsilon_{\nu}\cdot l$ with spectral index $\alpha=(1-\gamma_{e})/2$
on the high-frequency side. Typical values observed for radio galaxies are
$\alpha\sim-0.7$. The peak is described by ($\nu_{m}$,$S_{m}$), where $S_{m}$
is the intersection of the two power laws, i.e., an extrapolation that is not
actually reached by the SED. Figure 10 shows some example synchrotron spectra.

Figure 10: Example SEDs derived from Equations 10 and 14 to 18. SEDs shown in
green and blue were derived assuming the same ratio $\eta=\theta/B^{1/4}$ but
different electron densities $n_{e}$. For a given $\eta$, the self-absorption
turnover with the coordinates ($\nu_{m}$, $S_{m}$) lies on one of the dashed
grey lines having a power-law index of 2.5 (the slope of the optically thick
part of the spectrum), independent of $n_{e}$ or the optically thin flux-
density levels. The purple dashed line marks the trajectory of the turnover
under adiabatic expansion, i.e., the red SED could be the result of
adiabatically expanding a synchrotron sphere that starts off with the green
SED. The slope of the turnover trajectory related to adiabatic expansion can
be significantly flatter than the iso-$\eta$ trajectories (the exact slope
being given by Equation 31).

Equation 11 for the SSC flux density makes the assumption of spherical
symmetry. As mentioned, Equation 14 is formally derived for a slab of
material, i.e., if $\Delta\Omega=\pi/4\cdot(L/d)^{2}$, $S_{\nu}$ of Equation
10 describes a cylinder homogeneously filled with electrons. However, Band &
Grindlay (1985) demonstrated that the difference between spectra derived with
a slab approximation and with spherical symmetry is very small. (Of course,
$n_{0}$ needs to be multiplied by $1.5$ in order to account for the ratio
between the volume of a cylinder and the volume of a sphere.) Therefore we
will proceed with the slab approximation and use the correction factor of
$1.5$ when calculating electron densities from flux densities.

Moffet (1975) derived a central equation to link the self-absorption turnover
of the synchrotron spectrum described by the observable quantities $\nu_{m}$
and $S_{m}$ with the physical parameters $\theta$ and $B$:

$\nu_{m}^{5/2}=c^{2}C_{1}^{1/2}F^{\prime}S_{m}B^{1/2}/(\pi/4\cdot\theta^{2})$
(19)

with

$F^{\prime}=\frac{(\gamma_{e}+2)\cdot
G^{\prime}(\nu/\nu_{1},\nu/\nu_{1},\gamma_{e}+1)}{G^{\prime}(\nu/\nu_{1},\nu/\nu_{1},\gamma_{e})}\;\;.$
(20)

This equation is similar to an equation of state and allows us to describe how
the observable quantities change under evolution of the physical quantities.
Because $\nu_{m}$ and $S_{m}$ are not observable in our case—the self
absorption turnover of the compact component responsible for the X-ray and NIR
variability is veiled by emission from the constant radio component—we express
$S_{m}$ in terms of the optically thin synchrotron flux density
$S_{\rm{thin}}$:

$\displaystyle S_{\rm{thin}}$
$\displaystyle=\pi/4\cdot\epsilon_{\nu_{\rm{NIR}}}\cdot
L\cdot\theta^{3}\;\;,~{}{\rm and}$ (21) $\displaystyle\nu_{m}^{5/2-\alpha}$
$\displaystyle=c^{2}C_{1}^{1/2}F^{\prime}S_{\rm{thin}}\cdot(\nu_{\rm{NIR}})^{-\alpha}B^{1/2}/(\pi/4\cdot\theta^{2})\;\;.$
(22)

Similarly, the SSC flux density defined in Equation 11

$\displaystyle S(E_{\rm{keV}})$
$\displaystyle\propto\nu_{m}^{(3\alpha-5)}S_{m}^{(4-2\alpha)}\theta^{(4\alpha-6)}$
(23)
$\displaystyle\propto\nu_{m}^{(3\alpha-5)-2\alpha(\alpha+2)}S_{\rm{thin}}^{(4-2\alpha)}\theta^{(4\alpha-3)}\;\;.$
(24)

We are now in position to discuss illustrative scenarios of source evolution
and their consequences for the dependence of the X-ray and submm flux
densities on $S_{\rm{thin}}$. While none of these is directly applicable to
the final model, they are useful to show the dependence of observables on
physical quantities. In the following, we define $\eta\equiv\theta/B^{1/4}$.

####  4.2.1 Case 1

Let us assume $\Delta\eta=0$, $\Delta S_{\rm{thin}}\neq 0$. In this case,
$S_{m}\propto(\nu_{m})^{2.5}$, i.e., the turnover moves along the iso-$\eta$
lines shown in Figure 10. These iso-$\eta$ lines have the same slope as the
optically thick part of the spectrum, and the flux densities at frequencies
below the peak don’t change. From Equation 22,

$\nu_{m}\propto(S_{\rm{thin}})^{\frac{2}{5-2\alpha}}\;\;,$ (25)

and with $\Delta\theta=0$ (i.e., no changes other than in $n_{0}$) and
assuming $\alpha=-1$,

$S(E_{\rm{keV}})\propto(S_{\rm{thin}})^{2}\;\;.$ (26)

This occurs because the increase in $\nu_{m}$ partially counters that in
$S_{m}$. This gives the dependence of $S(E_{\rm{keV}})$ on $S_{\rm{thin}}$ a
significantly lower exponent than the value 6 apparent from Equation 11.

####  4.2.2 Case 2

Next consider the case $\Delta\nu_{m}=0$, $\Delta S_{\rm{thin}}\neq 0$, i.e.,
a constant turnover frequency under changing $S_{\rm{thin}}$. In this case,
$S_{m}\propto S_{\rm{thin}}$, $S_{\rm{thick}}\propto S_{\rm{thin}}$, and

$\eta\propto(S_{\rm{thin}})^{1/2}\;\;.$ (27)

If the change in $\eta$ is mainly a change in $\theta$—a corresponding change
in $B$ would have to be higher by the fourth power—$\Delta B=0$ and
$\alpha=-1$ imply

$S(E_{\rm{keV}})\propto S_{\rm{thin}}\;\;.$ (28)

####  4.2.3 Case 3

As mentioned in the Introduction, a case of particular interest is source
evolution through adiabatic expansion. The fundamental assumptions have been
stated by Shklovskii (1960) in the context of supernova-remnant evolution:

$B\propto\theta^{-2}$ (29)

and

$E\propto\theta^{-1}$ (30)

as $\theta$ changes with time. In contrast to synchrotron cooling, cooling by
adiabatic expansion applies to electrons of all energies at the same rate set
by the expansion speed. van der Laan (1966) showed that in this case,

$S_{m}(\nu_{m})\propto\left(\nu_{m}\right)^{-\frac{7\gamma_{e}+3}{4\gamma_{e}+6}}\;\;.$
(31)

For reasonable values of $\gamma_{e}$, $S_{m}$ has a flatter index than the
iso-$\eta$ lines (Figure 10). Other relations are:

$\displaystyle\nu_{m}$
$\displaystyle\propto(S_{\rm{thin}})^{\frac{2}{5-2\alpha}}\cdot\eta^{\frac{-4}{5-2\alpha}}\;\;,$
(32) $\displaystyle B$
$\displaystyle\propto(S_{\rm{thin}})^{\frac{2}{3-4\alpha}}\;\;,\;\;{\rm and}$
(33) $\displaystyle\theta$
$\displaystyle\propto(S_{\rm{thin}})^{-\frac{1}{3-4\alpha}}\;\;.$ (34)

For the X-ray flux density:

$\displaystyle S(E_{\rm{keV}})$
$\displaystyle\propto\eta^{-4(5-3\alpha+2\alpha(\alpha-2))/(5-2\alpha)}$
$\displaystyle\cdot(S_{\rm{thin}})^{(4-2\alpha)+2(3\alpha-5-2\alpha(\alpha+2))/(5-2\alpha)}$
$\displaystyle\cdot\theta^{(4\alpha-6)}$ (35)

For $\alpha=-1$,

$S(E_{\rm{keV}})\propto(S_{\rm{thin}})^{1.71}\;\;.$ (37)

While $S_{m}$ decreases with $S_{\rm{thin}}$, $\nu_{m}$ decreases as well.
Therefore the optically thick flux densities show a temporary increase some
time after the time of maximum flux density at optically thin frequencies.

In all three cases discussed here, under the physical constraints expressed in
Equation 22, the actual dependence of the SSC flux density on $S_{\rm{thin}}$
is only weakly non-linear with an exponent of 1–2 and does not come close to
the apparent exponent of $(4-2\alpha)$ in Equation 11.

None of the discussed cases of source evolution can by itself reproduce the
submm to NIR phenomenology. Case 1 does not predict NIR-correlated submm
variability at all. Case 2, a scenario of non-adiabatic compression and
expansion, predicts a direct proportionality of submm flux-density changes to
$S_{\rm{thin}}$, which is not observed. Case 3 describes only the cooling of
the synchrotron source by expansion, i.e., decaying $S_{\rm{thin}}$, and the
resulting propagation of delayed peaks towards longer observing wavelengths.
The opposite, increasing $S_{\rm{thin}}$ under adiabatic compression, would be
difficult to understand physically. Furthermore, adiabatic compression would
result in leading maxima at submm wavelengths, also not observed.

###  4.3 A Simple Source Model

To generate submm light curves that correspond qualitatively and
quantitatively to the observed data, we here propose a simple three-step
process of electron injection, compression of magnetic field lines, and
expansion. For simplicity, we model this three-step process as a cyclic
process in a single zone. Analysis of this process makes the assumption that
if $S_{\rm{thin}}$ is rising fast, the variability is injection dominated, and
the opacity is developing according to Case 1. If $S_{\rm{thin}}$ is falling
quickly, the source is adiabatically expanding according to Case 3. Just
before the time of peak flux density, the source is subjected to an episode of
(non-adiabatic) compression and particle escape, i.e., a compression of the
magnetic field lines mainly without electron heating (similar to Case 2 but
with changing $B$). Here, particle escape is necessary because close to the
peak, $S_{\rm{thin}}$ changes little, but in our model $B$ increases. Because
$\epsilon_{\nu}$ depends on $B$ and on the product $n_{0}\theta^{3}$ (which is
constant under particle conservation), particles must escape.

Our three-step source evolution is illustrated in Figure 11, including
transitions between injection and compression and compression and expansion.
The model qualitatively predicts submm light curves correlated with the NIR
and shows a range of possible delays between the two bands. It is, indeed, the
simplest source evolution model that can reproduce the observed phenomenology
without generating artificial symmetries in the light curves or too strict
correlations.

Figure 11: Illustration of the evolution of our simple model for the compact component in Sgr A*, denoted by the dashed circle. Upper panels 1–6 are a cartoon of a possible sequence of injection, compression, and expansion. The steps are: 1) and 2) injection of electrons with a non-thermal energy distribution into a spherical region with a uniform magnetic flux density giving detectable submm and NIR emission. 3) and 4) further injection, compression, and increase of magnetic flux under electron escape giving increased self-Compton scattering efficiency and detectable X-ray emission. 5) minimum size with highest magnetic flux density under ongoing injection followed shortly after by maximum in electron density and in NIR and X-ray emission. 6) adiabatic expansion, no injection giving maximum in submm emission. The funnels on the two sides of the synchrotron sphere are idealized and not implying actual symmetry of the inflow or outflow of particles. In particular, particle escape under compression might well be isotropic. The bottom panel shows the corresponding cycle of the self-absorption turnover position in the $S_{m}$–$\nu_{m}$ plane. Steps are labeled as above, and the diagonal dashed lines are lines of constant $\eta$.  Table 3: Constants Constant | Value | Unit
---|---|---
$C_{1}$ | $16.08$ | MHz $\mu$G-1 GeV-2
$C_{2}$ | $1.16540\times 10^{-26}$ | GeV $\mu$G-1 sr-1
$C_{3}$ | $1.19732\times 10^{-10}$ | yr-1 $\mu$G-2 GeV-1
$C_{4}$ | $1.358688\times 10^{-10}$ | mas-1
$C_{5}$ | $6.241506\times 10^{-17}$ | GeV Jy-1
$C_{6}$ | $0.0098087$ | mas $(1.2\times 10^{10}~{}{\rm m})^{-1}$
$C_{7}$ | $525949$ | minutes/yr
${R}_{S}$ | $1.2\times 10^{10}$ | m
$D_{s}$ | $2.523\times 10^{20}$ | m
$k_{X}$ | 3.39 | (factor)
$A_{K}$ | 2.46 | magnitudes
$A_{M}$ | 1.00 | magnitudes

###  4.4 Model Implementation

We implemented the Section 4.3 synchrotron–SSC mechanism in a semi-empirical
model. The model is semi-empirical because

  * •

the opacity evolution model is derived from submm–NIR phenomenology;

  * •

it is based on two generic, log-normal red-noise processes drawn from the same
random numbers according to Section 3.5 for which we can provide only
empirical reasoning.

One of the red-noise processes is the slow process drawn from a PSD with the
parameters $\gamma_{\rm{slow}}$ and $f_{b,\rm{slow}}$:

$S_{\rm{thin}}(t)=\exp[\sigma_{\rm{slow}}y_{\rm{slow}}(t)+\mu_{\rm{slow}}]\;\rm{Jy}\;\;,$
(38)

with $y_{\rm{slow}}(t)$ the slow Gaussian process and $\mu_{\rm{slow}}$ and
$\sigma_{\rm{slow}}$ the log-normal parameters. This process represents the
variability of the optically thin part of the synchrotron spectrum and the
timescales set by the injection process and expansion cooling. We chose a log-
normal PDF because the synchrotron equations are power laws, and a log-normal
process remains log-normal under multiplication and exponentiation. The other
red-noise process is the fast process drawn from a PSD with the parameters
$\gamma_{\rm{fast}}$ and $f_{b,\rm{fast}}$:

$\nu_{2}(t)=\nu_{\rm{NIR}}\cdot\exp[\sigma_{\rm{fast}}y_{\rm{fast}}(t)+\mu_{\rm{fast}}]+\max(\nu_{m},\nu_{\rm{min}})\;\;,$
(39)

with $y_{\rm{fast}}(t)$ the fast Gaussian process and $\mu_{\rm{fast}}$ and
$\sigma_{\rm{fast}}$ its log-normal parameters. This process represents the
variable location of the synchrotron cooling cutoff that is the result of the
two competing processes at the high frequency tail of the synchrotron
spectrum: the tail of the injection spectrum and the cooling through
synchrotron emission. In this case, we used a log-normal PDF as well to ensure
$\nu_{2}(t)>0$. However, it is a three-parameter log-normal PDF to account for
the fact that the cooling cutoff frequency should be larger than
$\max(\nu_{m},\nu_{\rm{min}})$ at all times, i.e., larger than both the self-
absorption turnover frequency and the transition frequency from synchrotron
cooling to expansion cooling. The latter can be derived from the equation for
the critical frequency for a given electron energy:

$\nu_{c}=C_{1}BE^{2}$ (40)

and the equation for the time after which an electron of initial energy
$E_{0}$ has cooled to $1/e$ of its initial energy,

$t_{1/e}=(C_{3}B^{2}E_{0})^{-1}(e-1)\cdot C_{7}\;\;.$ (41)

We then define $\nu_{\rm{min}}$ as the frequency in GHz where the $1/e$
cooling time is equal to the correlation time scale of the slow process:

$\nu_{\rm{min}}=10^{-3}C_{1}C_{7}^{2}\cdot\frac{f_{b,\rm{slow}}^{2}\;(e-1)^{2}}{(C_{3})^{2}B^{3}}\;\;.$
(42)

Figure 2 of Moffet (1975) shows the spectrum of an individual electron. The
definition above guarantees that the synchrotron cooling break $\nu_{2}$ can
occur only at frequencies that have enough time to cool through synchrotron
emission during an episode of injection, compression, and expansion (with a
typical duration of $1/f_{b,\rm{slow}}$). We define:

$\beta[\eta(t)]\equiv c^{2}\cdot(C_{1})^{1/2}\cdot[\eta(t)\cdot
C_{4}]^{-2}\cdot F^{\prime}\times 10^{-12}\cdot C_{5}\;\;,$ (43)

with

$\eta(t)=\frac{\theta(t)}{B(t)^{1/4}}\;\;$ (44)

and $F^{\prime}$ from Equation 20. We can now derive expressions for the
combined quantities

$\kappa_{\nu}L=\nu^{-(\gamma_{e}+4)/2}\cdot\beta[\eta(t)]\cdot
S_{\rm{thin}}(t)\cdot\nu_{\rm{thin}}^{(\gamma_{e}-1)/2}$ (45)

and

$\frac{\epsilon_{\nu}}{\kappa_{\nu}}\cdot\Delta\Omega=\frac{\nu^{5/2}}{\beta[\eta(t)]}\;\;.$
(46)

With these and $\nu_{2}(t)$, Equations 14 and 10 give the time-variable
synchrotron flux density at each frequency. Similarly, $\eta(t)$ and
$S_{\rm{thin}}(t)$ give

$\displaystyle\nu_{m}(t)$ $\displaystyle=\,0.001\cdot$ (47)
$\displaystyle\left[C_{5}\cdot\beta[\eta(t)]\cdot
S_{\rm{thin}}(t)\left(\nu_{\rm{NIR}}\cdot
1000\right)^{\frac{\gamma_{e}-1}{2}}\right]^{{2}/({\gamma_{e}+4})}$

and

$S_{m}(t)=S_{\rm{thin}}(t)\cdot\left[\nu_{m}(t)/\nu_{\rm{NIR}}\right]^{-\frac{\gamma_{e}-1}{2}}\;\;,$
(48)

and $\theta(t)$ and Equation 11 give the power-law section of the SSC SED.

We parameterized the cyclic source evolution model by linear functions with
variable slopes in the logarithmic $(S_{\rm{thin}},\theta)$ and
$(S_{\rm{thin}},B)$ planes (i.e., power-laws with variable indices):

$\displaystyle B(t)$
$\displaystyle=B_{0}\left[\frac{S_{\rm{thin}}(t)}{S_{\rm{thin},0}}\right]^{u(t)}$
(49) $\displaystyle u(t)$
$\displaystyle=\frac{1}{2}\left(\frac{2}{2\gamma_{\rm{e}^{-}}+1}\right)$ (50)
$\displaystyle\cdot\left\\{\tanh\left[\frac{\dot{S}_{\rm{thin}}(t)}{2S_{\rm{thin}}(t)}\cdot\rm{min}\right]+1\right\\}$
(51) $\displaystyle L(t)$
$\displaystyle=L_{0}\cdot\left[\frac{S_{\rm{thin}}(t)}{S_{\rm{thin},0}}\right]^{-v(t)}$
(52) $\displaystyle v(t)$
$\displaystyle=\frac{1}{2}\left(\frac{1}{2\gamma_{\rm{e}^{-}}+1}\right)$ (53)
$\displaystyle\cdot\left\\{\tanh\left[\frac{\dot{S}_{\rm{thin}}(t)}{2S_{\rm{thin}}(t)}\cdot\rm{min}\right]+1\right\\}\;\;.$
(54)

(The super-dot notation means time derivative in units of minutes-1.) In
Equation 49 for $B$, the power-law index $u$ was constrained to the interval
of $[0,\frac{2}{2\gamma_{\rm{e}^{-}}+1}]$. On average it is
$\frac{1}{2\gamma_{\rm{e}^{-}}+1}$ with deviations to either a larger or
smaller index depending on the fractional derivative
$\frac{\dot{S}_{\rm{thin}}(t)}{S_{\rm{thin}}(t)}$ in units of $\rm{min}^{-1}$.
Similarly, the index $v$ in Equation 52 for $L$ was constrained to
$[0,\frac{1}{2\gamma_{\rm{e}^{-}}+1}]$, averaging
$\frac{1}{4\gamma_{\rm{e}^{-}}+2}$. $B_{0}$ and $\theta_{0}$ are scaling
parameters with units G and $R_{S}$, respectively. This cyclic model is a
simple empirical scenario based on the cases in Section 4.2. The extreme Cases
1 and 3 are asymptotically approached for very rapid changes. Depending on
$\dot{S}_{\rm{thin}}(t)$, our implementation allows for mixing of the
different cases (e.g., slight compression during injection, etc.) and smooth
transitions between the three steps of the source evolution. Figure 12 shows
trajectories in the $(S_{\rm{thin}},\theta_{0})$ plane, the
$(S_{\rm{thin}},B_{0})$ plane, and the $(\nu_{m}$,$S_{m})$ plane calculated
according to Equations 49 to 54.

Figure 12: Cyclic evolution of model parameters. Top panel shows source size
$L$ versus optically thin synchrotron flux density (neglecting interstellar
extinction), middle panel shows magnetic flux density versus the same, and
bottom panel shows self-absorption turnover in the $(\nu_{m}$,$S_{m})$ plane.
All parameter time sequences (blue loops) were derived according to Equations
49 to 54. The grey dashed lines indicate the trajectories for pure injection
($\Delta\eta=0$) and pure adiabatic expansion as labeled. The circular arrows
indicate the direction of time evolution in each panel.

The dependence of $S(E_{\rm{keV}})$ on $S_{\rm{thin}}$ for the peak flux
densities (and the flux-density minima) comes from setting
$\dot{S}_{\rm{thin}}(t)=0$. Then

$\displaystyle\nu_{m}$
$\displaystyle\propto(S_{\rm{thin}})^{\frac{2}{5-2\alpha}}\cdot\eta^{\frac{-4}{5-2\alpha}}\;\;,$
(55) $\displaystyle B$
$\displaystyle\propto(S_{\rm{thin}})^{\frac{1}{3-4\alpha}}\;\;,$ (56)
$\displaystyle\theta$
$\displaystyle\propto(S_{\rm{thin}})^{-\frac{1}{(6-8\alpha)}}\;\;,$ (57)

and for $\alpha=-1$,

$\displaystyle S(E_{\rm{keV}})$
$\displaystyle\propto(S_{\rm{thin}})^{1.86}\;\;.$ (58)

With this model implementation, we are able to predict light curves at submm
and NIR light curves according to Equation 10 and X-ray light curves according
to Equation 11 and compare them via a suited distance function to the observed
light curves in all three wavelength regimes. A high-level block diagram for
generating light curves according to this model is shown in Figure 8b. For the
X-rays, to convert predicted flux densities to count rates, we used the
conversion from the pileup-free count rate to the absorbed energy flux in the
2–8 keV band $c_{\rm{conv}}=2.81\times 10^{-11}$ ergs cm-2 s-1 cps-1 (Yuan &
Wang, 2016). To account for interstellar extinction, we multiplied by 3.39,
the effective correction value derived from the luminosity calculated from
$c_{\rm{conv}}$ and the unabsorbed luminosity (Yuan & Wang, 2016)) and used
the relation

$\Lambda(t)=\frac{S(E_{\rm{keV}})}{\rm{Jy}}\frac{\nu_{8\rm{keV}}^{\alpha}-\nu_{2\rm{keV}}^{\alpha}}{\nu_{5\rm{keV}}^{\alpha}(1+\alpha)\cdot
2.81\cdot 3.39\times 10^{12}}~{}\rm{cps}\;\;.$ (59)

Observed count rates were then calculated according to Equations 3 to 6, but
we assumed that the steady X-ray background is equal for all three detector
modes of Chandra (i.e., $\chi_{I}=\chi_{G}=\chi_{S}$). We used $e_{I}$ rather
than $e_{G}$ as the reference effective area for Chandra.

The free parameters in our analysis are listed in Table 2. In order to compare
the model predictions to observed NIR flux density, we applied extinction
magnitudes $A_{{K}}$ and $A_{\rm{M}}$ given in Table 3. Because these values
have significant uncertainties, we let the fitting modify the extinction
magnitudes via parameters $\Delta A_{{K}}$, and $\Delta A_{{M}}$. We set
parameter
$S_{\rm{thin},0}=0.1\cdot\exp(\mu_{\rm{slow}}-\sigma_{\rm{slow}}^{2})$ Jy, the
mode of the flux-density distribution of $S_{\rm{thin}}$. This lets the fitter
set $S_{\rm{thin}}(t)$ wherever $\mu_{\rm{slow}}$ and $\sigma_{\rm{slow}}$ say
it should be while not having $B(t)$ (Equation 49) or $\theta(t)$ (Equation
52) take on improbable values. The exponent of $0.5$ in the exponential term
of Equation 10 and the factor 0.5 in the argument of the $\tanh$ in Equations
51 and 54 are fiducial parameters which we did not fit. A mismatch in the
exponent can be absorbed by the parameters $\mu_{\rm{fast}}$ and
$\sigma_{\rm{fast}}$, at least where
$\nu_{2}(t)\gg\max(\nu_{m},\nu_{\rm{min}})$. Where $\nu_{2}(t)$ comes close to
its minimum, we estimate the error to be not greater than $10\%$. However,
close to the minimum of $\nu_{2}(t)$, the NIR flux density is far below the
detection limit, and the variability is entirely dominated by the measurement
noise. The factor of 0.5 in the argument of the $\tanh$ gives reasonable
results, and testing showed that the fitting results of the other parameters
are not strongly influenced by its exact value. To fit for this parameter we
would need detailed information on the covariance at several mm to submm
wavelengths and a suited distance function. Our data are not sufficient to
constrain this factor.

Finally, Equations 51 and 54 make use of the derivative
$\dot{S}_{\rm{thin}}(t)$. We take advantage of the fact that we are generating
the mock light curves via FFT from Fourier coefficients. To compute the
derivative from the same Fourier coefficients, we can use the fundamental
relation from Fourier theory:

$F\left[\frac{{d}}{{d}x}g(x)\right]=2\pi if_{x}G(F_{x})\;\;$ (60)

with $F$ the Fourier transform and $G(f_{x})=F[g(x)]$. However, this equation
cannot simply be applied to the discrete case, and for that we followed the
method of wavenumber modification (Sunaina et al., 2018).

###  4.5 The Distance Function

As in Section 3, to use ABC we need a distance function. Our distance function
used several components to guide the ABC algorithm to a valid set of
posteriors:

  * •

The distance between the predicted and the observed structure functions as
defined in Equation 8. We calculated the structure functions and distances for
the three NIR, two submm, and three X-ray datasets separately. We added these
eight distances multiplied by normalized weights, which were the quadratic
difference between the logarithm of the maximum value and the logarithm of the
minimum value of each of the eight observed structure functions $V_{j}$:

$D_{\rm{sf}}=\frac{\sum_{j}\left[\max{\log(V_{j})}-\min{\log(V_{j})}\right]^{2}\phi_{j}}{\sum_{j}\left[\max{\log(V_{j})}-\min{\log(V_{j})}\right]^{2}}\;\;$
(61)

with $D_{\rm{sf}}$ the total distance calculated from structure functions and
$\phi_{j}$ the structure function distance for each dataset. For each
structure function, we used uniform weights $w_{i}=1$ for all bins except for
the last bins of the NIR and X-ray data (i.e., the bin with $\tau>20$ minutes
for the Keck data and the one with $\tau>128$ minutes for all other NIR and
X-ray data), which had $w_{i}=3$. This approach led to uniformly converging
fits at all time lags and in all bands. (See also the discussion by Witzel et
al. 2018, Appendix B.2.)

  * •

To help to constrain the log-normal PDF of $\nu_{2}(t)$, we determined two
other quantities to describe the empirical distribution of observed $K$-band
flux densities: the skewness of the distribution and the fraction of positive
flux densities. Both values are sensitive to the white noise characteristics,
and we determined them only for the VLT dataset, which is the larger of the
two $K$-band datasets. The fraction of positive flux densities is 0.9. This
value was derived from the most recent analysis of the intrinsic-flux density
distribution from the VLTI/GRAVITY interferometer (Gravity Collaboration et
al. 2020c). VLTI/GRAVITY data show significantly less noise than the data from
single-dish telescopes, and its 1 mas angular resolution makes it virtually
free of source confusion. The VLTI/GRAVITY data show a median flux density of
$1.1\pm 0.3$ mJy dereddened, and the peak of the distribution is ${\sim}0.4$
mJy dereddened. While the empirical distribution of VLTI/GRAVITY flux
densities is always positive, convolving it with a Gaussian noise of
$\sigma_{\rm{NACO}}\approx 0.3$ mJy dereddened will result in a tail of
$\sim$10% negative flux densities. We counted the fraction of $n_{\rm{>0}}/n$
points in each VLT mock light curve above zero and the quadratic difference:

$D_{>0}=(n_{\rm{>0}}/n-0.9)^{2}\;\;.$ (62)

For the skewness,

$\displaystyle\rm{Sk}_{\rm{VLT,obs}}$
$\displaystyle=\frac{\frac{1}{n}\sum_{i=1}^{n}(S_{i}-\overline{S})^{3}}{\left[{\frac{1}{n-1}\sum_{i=1}^{n}(S_{i}-\overline{S})^{2}}\right]^{3/2}}$
$\displaystyle=3.44\;\;$ (63)

with $S_{i}$ the observed flux densities and $\overline{S}$ the mean flux
density. For each VLT mock light curve, we calculated the skewness and the
quadratic difference

$D_{\rm{sk}}=(\rm{Sk}_{\rm{VLT,sim}}-\rm{Sk}_{\rm{VLT,obs}})^{2}\;\;$ (64)

and added the result to the distance function.

  * •

Following Witzel et al. (2018, their discussion of Case 3 in Section 4.4), we
additionally used spectral information from simultaneous $K$\- and $M$-band
data. At an average (observed) flux density of $S_{{K}}=0.15$ mJy, Witzel et
al. found the ratio of (observed) NIR flux densities at $M$ and $K$ bands to
be $\Re_{\rm{obs}}=12\pm 0.5$. We calculated this ratio for each parameter set
from our mock data. Because here we are not using simple log-normal PDFs for
modeling the flux densities (as Witzel et al. did), we cannot derive this
ratio analytically from model parameters. Instead, we simulated for each
parameter set equally sampled, measurement-noise-free $K$\- and $M$-band light
curves of 10 000 minutes duration each and determined $\Re[M/K,S_{{K}}]$ as
the average ratio over a suited flux-density range. The flux-density range of
the Witzel et al. simultaneous datasets is $\sim$0.07–0.27 mJy. However, it is
difficult to determine background flux-density levels for single-dish data.
After comparing the median of the NACO flux-density distribution with the
newest VLTI/GRAVITY study (Gravity Collaboration et al. 2020c), we adopted a
background correction factor based on a flux-density range of [0.0–0.17] mJy.
We then defined the distance

$D_{\rm{sp}}=(\Re[M/K,S_{\rm{K}}]-\Re_{\rm{obs}}+\mathcal{N}(0,0.5^{2}))^{2}\;\;,$
(65)

with $\mathcal{N}(0,0.5^{2})$ a normal random number to account for the
uncertainty in $\Re_{\rm{obs}}$.

The final distance function was

$D_{\rm{tot}}=D_{\rm{sf}}+10\cdot D_{>0}+D_{\rm{sk}}+D_{\rm{sp}}/50\;\;.$ (66)

The factors of 10 and 50 are empirical and ensured that none of the additional
distance terms starts to dominate $D_{\rm{sf}}$. This distance function does
not include metrics to quantify the degree of correlation between wavelengths,
and we restricted the $230$ GHz structure function to time lags $<$50 minutes.
We will come back to these two points in the discussion.

Figure 13: Performance of the ABC sampler. Ordinate is the acceptance rate,
and abscissa is the 45% quantile of all distances ($D_{\rm tot}$, Eq. 66) of
the particular iteration. Towards the smallest distances the sample finds $<$1
in 10 000 light curves that obeys the distance threshold. Figure 14: Results
for the Bayesian structure function fit for the synchrotron–SSC model.
Contours show the joint posterior probability density for each parameter pair,
and panels along the right upper edge show histograms of the marginalized
posterior of each parameter defined in Table 2. For each histogram, the dashed
lines mark the 16%, 50%, and 84% quantiles. Instrumental $\sigma$ parameters
are as observed, not extinction-corrected.

##  5 Results

We have implemented the described model and distance function in our C++ code.
The analytic nature of the model allows us to calculate many realizations as
large as the observed dataset in reasonable time. However, because of the
iterative nature of the ABC and the need for a sufficiently large particle
system, we had to run our code on the VLBI compute cluster of the Max-Planck-
Institut für Radioastronomie. The model and the distance function are the
result of more than fifty test runs that all had to be executed over several
days on 200–400 cores (10–20 nodes). Over the course of two years, we estimate
the total CPU time (including all runs necessary for testing, developing, and
implementing the model) to be ${\sim}10^{6}$ hours. To analyze the results, we
also implemented a Python version of our model. This version allows us to
quickly visualize and diagnose the ABC results, including a tool for
generating animations of the timing of the SED, the light curves, and the
derived physical parameters.

Figure 15: Structure functions (Eq. 7) for the eight datasets. Panels from top
to bottom show NIR, X-ray, and submm, all as observed, not corrected for
extinction or pileup. Solid lines show the observed SFs (as presented in
Figure 4), and corresponding dashed curves show the median of 10 000 SFs for
the respective mock datasets. The shaded envelopes denote the model $68\%$
credible intervals for each time lag. The vertical dashed line in the top
panel marks the derived correlation timescale $\tau_{b}=1/f_{b}$. Pairs of
horizontal short-dashed lines, color-coded for each of the NIR instruments,
mark their two respective white noise levels. The lower line of each pair
indicates the measurement noise (Witzel et al. 2018, Eq. 3) and the upper line
the intrinsic red noise of the Sgr A* variability when sampled at timescales
${\gg}\tau_{b}$ combined with measurement noise (Witzel et al. 2018, Eq. 4).
(The upper lines for NaCo and NIRC2 are nearly indistinguishable.) The details
of generating the SFs, including the choice of time-lag ranges, are described
in Section 4.5. The slope of each SF relates to the slope of the underlying
PSD but also depends on the underlying white noise level. Figure 16: Credible
contours (68%, 95%, 99%) for the parameters $\gamma_{1}$ and $f_{b}$. Green
and orange show the PSD parameters of our synchrotron–SSC model (Section 4 and
Table 2), green showing the slow process (identified with X-rays) and orange
the fast process (identified with NIR). Blue contours are for our X-ray-only
analysis (Section 3.5), and red contours are for the NIR analysis of (Witzel
et al., 2018, their Case 3). All posteriors have been marginalized over all
other parameters. The horizontal grey dotted line marks the frequency that
corresponds to a periodic timescale of 50 minutes, which is the orbital
timescale of the astrometric motion measured by Gravity Collaboration et al.
(2018b).

Table 2 lists the priors of all parameters. Priors for all physical parameters
except $B_{0}$ and $L_{0}$ were flat and wide. For $B_{0}$ and $L_{0}$, we
used log-flat priors. Parameters for the log-normal red noise processes were
also flat and wide. We labeled these parameters “slow” and “fast,” but we used
the exact same priors for both PSDs. In other words, we did not force the PSDs
to take on the filter relation discussed in Section 3.5. For parameters
specifying instrumental characteristics, we used narrow Gaussian priors.

The final ABC run presented here is the result of 33 iterations, 32 with a
particle size $n=500$ and the last with $n=5000$. (The model computations had
to use fewer particles than the X-ray-only model because submm and NIR light
curves had to be computed in addition to X-rays.) Figure 13 shows the
evolution of the distance and the acceptance rate. We reached
$D_{\rm{tot}}<0.09$ at an acceptance rate of $5\times 10^{-5}$. (This means
the final iteration generated 204 081 633 light curves in order find 10 000 it
could accept.) The local maximum in the acceptance near $D_{\rm{tot}}<0.2$ is
a hint that at this stage, the algorithm still “learned” relevant information.

The posterior estimates of all parameters converged to well-constrained,
peaked distributions. Figure 14 shows the posteriors and pairwise correlations
of all parameters. The medians and 1$\sigma$ credible intervals are listed in
Table 2. Figure 15 shows the structure functions with the corresponding
1$\sigma$ envelopes drawn from the posterior. The fit describes the observed
structure functions well with only the structure function of the Spitzer/IRAC
data lying partly outside the 1$\sigma$ envelope. This discrepancy is not
statistically significant, and part of it could be caused by our exact choice
of synchrotron cooling cutoff (Equation 10). The true shape of this cutoff
depends on fine details of the process, e.g., the injection spectrum and the
exact cooling mechanism. These depend on the geometry of the magnetic field
lines and other specifics. The shape we chose predicts a slightly higher
variability in the middle range of time lags at 4.5 $\micron$. In the absence
of more specific information on the coevolution of the spectral indices in the
IR bands, it is not possible to determine a more realistic scenario.

The physical parameters of the model are tightly constrained:
$\gamma_{e}=2.95\pm 0.2$, $B_{0}=8.6^{+1.8}_{-1.4}$ G, and $L_{0}=2.7\pm
0.3~{}R_{S}$. Figure 16 shows the posteriors of the PSD parameters. The PSD
slope of the slow process comes out a little steeper in the combined model
than when analyzing the X-ray data alone but not significantly so. A bigger
difference is between the fast and slow process timescales and the timescale
from the earlier NIR analysis. However, both timescales ($135^{+26}_{-20}$
minutes and $82^{+28}_{-16}$ minutes, respectively) are consistent with the
earliest report of the NIR break timescale of $154^{+124}_{-87}$ by Meyer et
al. (2009).

In order to determine posteriors of observable and physical quantities, we
used the model to construct 1000 mock light curves of 700 minutes duration in
all bands. The results are given in Table 2. The mock light curves were then
“observed” at uniform intervals, unlike the real data. This may give the
observed data higher variance than predicted by the mock data. Figure 17
compares the posterior of $(\log{S_{\rm{NIR}}},\log{S_{{X}}})$ to observed
flux-density pairs. Figure 18 shows calculated NIR and submm spectra indices.
The model NIR spectral index becomes flatter towards brighter flux densities
with a maximum around $\alpha_{\rm{NIR}}=-1$, while the submm spectral index
is a wide scatter cloud terminated sharply at $\alpha_{\rm{submm}}=+2.5$.
These are as expected for a synchrotron–SSC source.

Figure 17: Probability distributions of simultaneous NIR and X-ray
measurements. Grey points show distributions from 1000 mock multi-wavelength
light curves (Section 5), and lines show their 68% and 95% credible contours.
Red points show real measurements from simultaneous observations. All values,
X-ray and NIR, mock and real, are intrinsic with the real data corrected for
extinction by the canonical factors given in Table 3. Upper panel shows
$[S_{{K}},S(E_{\rm{keV}})]$ with data from Fazio et al. (2018). The large
orange data point shows the peak value of the NIR flare and its X-ray
counterpart observed by those authors (their Table 3). The even more extreme
event observed by Do et al. 2019 had no simultaneous X-ray observations. At a
corrected $S_{2.2~{}\micron}=60$ mJy, that observation would be just to the
right of the plot area. The bottom panel shows $[S_{{M}},S(E_{\rm{keV}})]$
with observed data from Boyce et al. (2019). To be consistent with the
presentation of the X-ray data in the respective publications, the X-ray count
rates are shown in 600 s bins for the upper panel and 300 s bins for the
lower. Even in those bin widths, X-ray counts can be zero. In order to show a
continuous posterior in log space, we added $\chi_{S}=0.001$ cps to the count
rate and a Gaussian random variable of $\sigma=0.18$ to the values of the
logarithm to dither both the mock and observed X-ray data. Figure 18:
Posterior distributions calculated from mock light curves (Section 5). Points
in the upper panel show NIR spectral index versus flux density
$S_{\rm{2.2~{}\micron}}$, and those in the lower panel shows the submm
spectral index versus $S_{\rm{230~{}GHz}}$. Observations have to be corrected
for extinction before comparison with quantities shown here.

Figure 19 shows the posterior distributions of physical and observable
parameters based on the models and mock light curves. The delay (panel a) was
computed from the maximum of the cross-correlation function of the 230 GHz
flux densities with $S_{\rm{thin}}$ for each 700 minute light curve. We
accepted only delays with peaks at least a factor of 5 higher than the
standard deviation of the cross-correlation function. The high-energy electron
density (panel d) was calculated according to:

$\displaystyle n_{0}(t)$ $\displaystyle=\frac{12\cdot
S_{\rm{thin}}(t)}{\pi\theta^{3}\cdot
d}\frac{(\nu_{\rm{NIR}}/C_{1})^{\frac{\gamma_{e}-1}{2}}}{C_{2}\cdot
B^{\frac{\gamma+1}{2}}G^{\prime}(\infty,0,\gamma_{e})}\;\;,$ (67)
$\displaystyle n_{e}(t)$
$\displaystyle=\frac{n_{0}(t)}{(1-\gamma_{e})}\cdot\left\\{\left[\frac{\nu_{2}(t)}{C_{1}\cdot
B}\right]^{(1-\gamma_{e})}-E_{\rm{min}}^{(1-\gamma_{e})}\right\\}$
$\displaystyle\mathrm{for}\;\;\gamma_{e}>1\;\;.$ (68)

All other quantities were calculated according to the equations in Section
4.4.

Figure 19: Posterior distributions of parameters for the synchrotron–SSC model
(Section 5) calculated from the mock light curves. These posteriors represent
both the variance propagated from the posterior as well as the variance from
the variability. The large skewness of the NIR flux-density distribution
biases the parameters towards small NIR fluxes. Panels show: a) time lag from
NIR to submm (positive means NIR leads), b) magnetic flux density, c) source
size in units of $R_{S}$, d) density of high-energy electrons in the source
volume, e) frequency of maximum synchrotron flux density, f) frequency of
cooling cutoff of synchrotron spectrum, g) intrinsic 4.5 to 2.2 µm spectral
index, and h) 230 to 345 GHz spectral index.

The synchrotron–SSC model gives reasonable values for physical and observable
parameters (Figure 19). The submm to NIR time delay (panel a) is a positively
skewed distribution in the interval of [0,100] minutes and peaks at ${\sim}20$
minutes. The magnetic flux density $B$ (panel b) varies between 6 and 30 G and
peaks around 13 G. The source diameter $L(t)$ (panel c) is a symmetric
distribution around $L=2.2~{}R_{S}$ with values between 1.3 and 3.4 $R_{S}$.
The electron density (panel d) is distributed around $n_{e}=4\times 10^{-7}$
$\rm{cm}^{-3}$ and varies by somewhat more than an order of magnitude in both
directions. The spectrum shows a self-absorption turnover (panel e) $\nu_{m}$
most often near 280 GHz but ranging between 100 GHz and 750 GHz. The cooling
cutoff (panel f) $\nu_{2}\sim 3\times 10^{4}$ GHz and rarely goes below
$5\times 10^{3}$ GHz or exceeds $10^{6}$ GHz. The distributions of
$\alpha_{\rm{NIR}}$ (panel g) and $\alpha_{\rm{submm}}$ (panel h) are
complementary to Figure 18. Table 2 lists medians and 1-$\sigma$ credible
levels for each of these quantities.

Even two uncorrelated light curves will sometimes show peaks near the same
time. The “false alarm probability” gives the probability of finding such
spurious correlations. at any time lag due to the white noise and the
autocorrelation (the red noise) in each of the two bands. Figure 3 shows the
95% false alarm probability for the submm–NIR light curves. This was derived
by calculating the DCC of the 230 GHz and 4.5 µm mock data in the same way as
for the observed data but using pairs of 700 minute light curves that do not
belong to each other. By modeling five datasets with the cadence of the
observed data, we generated 200 DCC functions from our 1000 light curves. At
each time lag we then determined the 95% quantile.

Figure 20 shows the posterior of complementary cumulative distribution
functions (CCDF) of $K$-band flux densities in comparison to the empirical
CCDF from NIRC2 and NACO data. For all fluxes, the model is 3$\sigma$
consistent with the observed data, and up to about 20 mJy (derredened) it is
2$\sigma$ consistent. However, the model does not predict the brightest flux
densities to be as frequent as we observe them.

Figure 20: Cumulative distribution functions (CCDFs) of Sgr A* 2.12 µm flux
densities (as observed, not corrected for extinction). The black upper line
shows the CDF observed by the VLT and Keck with the dashed section indicating
the CDF needed to match the single largest flux-density ever observed (Do et
al. 2019). The lower dashed blue line shows the median CDFs from the mock
light curves. Shaded areas show $68\%$, $95\%$, and 99% credible intervals
(Table 2).

Finally Figure 21 shows a snapshot of our animation. The animation is the most
concise way to illustrate the properties of our model. However, presenting
model light curves this way is based on picking a particular parameter set
from the posterior. It is the nature of the Bayesian approach that this
posterior includes particles that might not generate one or another aspect of
the observed data because, after all, the observed data could be an unlikely
realization of the underlying process. In order to allow the reader to become
more familiar with the phenomenology of the model and the presented posterior,
we are publishing the particle system (the “chain” of parameters), the
observed data, and our Python implementation of the model and of the animation
code. The details of the code and the data repository are described in
Appendix A.

Figure 21: Snapshot of the SED animation for one mock light curve (from
particle 1655) generated with our Python package (Appendix A). Panels a) and
b) show the observed SED at minute 732, a) in flux density and b) in
luminosity units. Orange lines show the synchrotron emission, blue the SSC
emission, and red their sum. The dashed line marks $\nu=230$ GHz, the short
grey solid line the typical range of NIR variability when the source flux is
above the detection limit, and the grey rectangle the typical range of peak
2–8 keV flux densities. Panel c) shows spectral index versus flux density for
minutes 683 to 732 of the mock light curve with the earlier points being more
transparent. Green shows NIR, and orange shows submm. Panel d) shows X-ray
flux “as measured” (cps) versus $S_{2.2~{}\micron}$, again for minutes 683 to
732 with the earlier points being more transparent. Panels e)–i) show time
series of various parameters as labeled. Left and right ordinate labels are
color-coded to match the quantities plotted. The vertical line marks minute
732. Points to the right of that are shown as intrinsic values, while in
panels e)–g), points to the left are shown “as measured,” i.e., with
measurement noise added. (This is not applicable to panels h and i.) Red
points represent a “measurement” at minute 732\. The X-ray count rate is
averaged over 600 s bins.

##  6 Discussion

The ABC algorithm results in constrained distributions of all parameters and
describes the structure functions of nearly all datasets within the 1$\sigma$
envelopes. The only exception is Spitzer/IRAC, where parts of the structure
functions show slightly less variability than the model at timescales $\sim$40
minutes. This deviation is not significant. The values of all parameters are
reasonable and consistent with previous analyses (e.g., Eckart et al. 2012)
and prior knowledge: magnetic flux density of $5<B<30$ G, source sizes
${<}4~{}R_{S}$, the upper limit derived from 1.3 mm VLBI observations
(Doeleman et al. 2008; Lu et al. 2018), and $n_{e}=4\times 10^{7}$ cm-3,
comparable to the solar corona. The effective collection area of each of the
Chandra modes and the white noise levels in the NIR and submm are consistent
with independently determined values. The quiescent count rate
$\chi=1.1^{+0.8}_{-0.5}\times 10^{-3}$ cps is identical with the background
count rate of comparison apertures close to Sgr A* (Yuan & Wang, 2016).
However, Yuan & Wang found the sum of background count rate, truly quiescent
level, and undetected variability $\chi_{\rm{I}}\approx 4.8\times 10^{-3}$
cps, $10\%$ of which they attributed to undetected, weak flares. In contrast,
in our analysis all X-ray photons other than the expected background rate can
be attributed to the SSC process. At times of low activity, this appears
quiescent.

The initial idea of two correlated processes related by a low-pass filter is
consistent with our final values for $(\gamma_{\rm{slow}},f_{b,\rm{slow}})$
and $(\gamma_{\rm{fast}},f_{b,\rm{fast}})$ as shown by Figure 16. While the
contours of our pre-analysis with generic log-normal flux-density
distributions overlap, the corresponding contours of the processes in the
synchrotron–SSC model do not. This is not the result of the priors, which are
equal for both PSDs. The break frequencies of both processes are
systematically higher than in the pre-analysis because of their dependence on
the underlying flux-density model, which is given by the radiative model. In
particular, the fast process is marginally consistent with a timescale of 50
minutes, which is the orbital timescale of the astrometric motion measured by
VLTI/GRAVITY (Gravity Collaboration et al. 2018b).

The consequence of the low-pass filter is the decorrelation of the two
processes towards higher frequencies. We can provide some qualitative
reasoning for this behavior. The variability is caused by injection of non-
thermal electrons whose energy distribution is truncated at some maximum
energy $E_{\rm{max}}$. This cuts off the synchrotron spectrum at a frequency
close to the NIR. The injected electrons are subject to both expansion cooling
and synchrotron cooling. The former affects the entire spectrum, while the
latter (in the time frame of injection and expansion cooling) affects only the
highest frequencies. The timescale of expansion cooling depends on the rather
low expansion velocity, but the synchrotron cooling timescale at frequencies
above the NIR is mere seconds to minutes. Therefore rapid fluctuations of the
injection process will be tracked by NIR emission, but the submm emission will
track only the long-term average. It should be possible to derive a time-
dependent model of the injection process from our semi-empirical variability
model, but this is beyond the scope of our analysis.

The distance function used for our ABC implementation does not use any measure
of correlation between the bands. Instead, the correlation is naturally built
into the model. The exact phenomenology of the correlation, however, depends
on the parameters that the ABC algorithm finds in the attempt to describe the
structure functions of all bands. Figure 17 shows a scatter plot of the
posterior NIR and X-ray flux-density pairs. The observations are in good
agreement with the posterior with the exception of a rare, very high $K$-band
flux density at a low X-ray level (Fazio et al., 2018, their Figure 3).
Estimating the probability of this event is difficult, given the different
cadences of the mock and real light curves, but the event is a hint that the
real source is more complicated than our purely cyclic model.

Boyce et al. (2019) offered another way to investigate the NIR to X-ray
correlation. Their cross-correlation analysis of the simultaneous Chandra and
Spitzer data found a tendency for the X-ray peaks to precede the NIR. However,
Boyce et al. could not claim a statistically significant delay given the large
uncertainties. Boyce et al. also presented an overview of earlier results,
many of which showed the X-ray to follow the NIR but again without high
significance. Our 700 minute mock light curves show no significant delay
between the X-ray and NIR bands.

Our model predicts correlation between the NIR and submm light curves as well.
We quantify this in terms of the cross-correlation between $S_{\rm{thin}}$ and
the corresponding 230 GHz mock light curves as shown in Figure 19a. Observed
delays range from almost synchronous variability (as observed by Fazio et al.
2018 and in this paper in the case of the 345 GHz SMA data from 2014 Jun 17)
to delays up to 90 minutes (as found by Eckart et al. 2008a). The model
distribution peaks at ${\sim}22$ minutes with a FWHM of ${\sim}10$ min, which
corresponds nicely to the result we obtained from our ${\sim}48$ hours of
synchronous NIR and 230 GHz data. The cross-correlation in Figure 3 shows a
pronounced peak at ${\sim}27$ min. The observed light curves have a maximum
correlation $\rm{R}=0.21$. However, the observed NIR flux densities, if our
model is correct, are not identical with the slow process $S_{\rm{thin}}(t)$
but include the faster varying effects of the cooling cutoff. In order to
approximate the quantity $S_{\rm{thin}}(t)$, we also present the cross-
correlation of the 230 GHz data with the logarithm of the low-pass-filtered
NIR flux densities. This results in a significant increase of the correlation
to $\rm{R}=0.43$. The 95% false alarm probability levels calculated from our
model and the posterior in Section 2 indicate that the data are consistent
with our model assumptions555However, we have not tested the significance of
the correlation independently (with, e.g., a Granger causality test). The
cross-correlation of Figure 3 cannot be used as independent evidence for the
existence of this type of correlation. All we can conclude here is that the
Spitzer–SMA synchronous dataset is consistent with a delay that is predicted
by our model..

Several papers have discussed expanding plasmon flare evolution models for Sgr
A*. The earliest, by Yusef-Zadeh et al. (2006a), described the temporal
evolution of flares in the cm wavelength regime. Eckart et al. (2006, 2008a,
2012) explained a possible correlation and delay between NIR and a submm
flares. Marrone et al. (2008) presented an example of simultaneous IR and
submm data with a delay (IR leading, submm following) of $20\pm 5$ minutes,
consistent with the data presented here and with our model and also discussed
adiabatic expansion. All these expanding plasmon models result in typical
expansion velocities of ${\sim}0.01c$, with $c$ the speed of light. Peak
expansion velocities derived from our time dependent model reach ${\sim}0.01c$
and are consistent with the earlier estimates.

In contrast to Yusef-Zadeh et al. (2006a), Falcke et al. (2009) interpreted
20–40 minute delays in the 20 to 40 GHz regime in the context of the frequency
dependence of VLBI sizes and saw evidence for a relativistic outflow.
Brinkerink et al. (2015) even followed the realtime progression of variability
maxima from 100 to 19 GHz and similarly derived relativistic outflow
velocities of up to ${\sim}0.77c$, i.e., a jet. Whether the similar delays
found at radio frequencies and in the submm to NIR regime indicate a direct
relation between the variability of those regimes or are merely coincidental
cannot be decided here. Also, we cannot determine whether the compact
component responsible for the fast, high-frequency variability is located in a
jet or in the accretion disk.

Neilsen et al. (2015) pointed out that the relation between NIR and X-ray
fluxes is expected to be mildly non-linear. Assuming power-law distributions
for the flux-density distributions in the NIR and the X-ray, they derived the
exponent of the power-law dependence of X-ray flux densities on NIR flux
densities (their Equations 8–11). In the case of a power-law index of
${\sim}2$ for the X-ray distribution and ${\sim}4$ for the NIR (Witzel et al.
2012; Neilsen et al. 2015), the dependence should be $S(E_{\rm{keV}})\propto
S_{\rm{NIR}}^{3}$. Our simple model predicts an exponent between 1.7 and 2.0
(Section 4.2). The higher non-linearity in the relation of observed NIR flux
densities to simultaneously measured X-ray flux densities can, however, be
explained by the exponential cooling cutoff (Equation 10).

The model predicts the self-absorption turnover to vary between 150 GHz and
700 GHz, yielding a wide range of spectral indices at these wavelengths. On
the other hand, in the NIR the spectal index is typically steeper than the
optically thin spectral index $\alpha=(1-\gamma_{e})/2$ due to the cooling
cutoff. Both regimes, one used by mm VLBI, the other by VLTI/GRAVITY, are
limited in the sense that they are subject to rapid changes in the state of
the radiative transfer: the submm by its regular changes between the optically
thick and thin branch of the spectrum which potentially causes rapid changes
in the source structure and the NIR by fast synchrotron cooling, which makes
it hard to follow the compact component of Sgr A* in its entire evolution.

The median brightness temperature at the self-absorption turnover
$T_{b}=1.8\times 10^{10}$ K. It varies between 80% and 90% of the
equipartition brightness temperature $T_{\rm{eq}}$. Assuming a pure electron
plasma and a Doppler factor of 1, this value suggests that the source is
energetically dominated by the magnetic flux.

The model spectral index in the NIR has an upper limit of
$\alpha_{\rm{NIR}}\approx-0.8$ with many bright flares showing spectral
indices as steep as $\alpha_{\rm{NIR}}\approx-1.2$ (Figure 20). The upper
limit of $\alpha_{\rm{NIR}}$ is defined by the spectral index
$\alpha=(1-\gamma_{e})/2$ of the optically thin spectrum without cooling
cutoff. The NIR spectral index reaches the optically thin limit when the
cooling cutoff is at high energies. In our model fit, $\gamma_{e}$ is
determined by the variance of the submm variability relative to the variance
of the NIR variability. $\alpha\approx-1.0$ is steeper than the canonical NIR
spectral index for bright flares of $\alpha_{\rm{NIR}}\approx-0.65$
(Eisenhauer et al. 2005; Ghez et al. 2005; Gillessen et al. 2006; Krabbe et
al. 2006; Hornstein et al. 2007; Bremer et al. 2011; Witzel et al. 2014).
However, spectral index analyses in the NIR have two main sources of
systematic uncertainty: the extinction correction and background flux levels.
While the extinction correction is precise enough and cannot explain this
discrepancy, unaccounted background flux makes a big difference. However, Do
et al. (2019) reported $\alpha_{\rm{NIR}}\approx-0.4$ (with canonical
extinction correction) at such high flux-density levels that background
contamination does not play a role. Such a flat spectral index could be the
result of the SSC scattering contributing to the NIR. The lower limit of the
power-law section of the SSC spectrum is given in Equation 12. As seen in
Figure 9, for source parameters discussed here, the SSC spectrum sets in
around NIR frequencies and might become dominant once in awhile.

Figure 22: Variable SED in comparison to overall SED of Sgr A*. The dashed
lines show snapshots of the SED from our model for three different
(dereddened) NIR flux-density levels. SED values are taken from Zylka et al.
(1992); Telesco et al. (1996); Falcke et al. (1998); Cotera et al. (1999); An
et al. (2005); Dodds-Eden et al. (2009); Schödel et al. (2011); Dodds-Eden et
al. (2011); Bower et al. (2015); Liu et al. (2016); Stone et al. (2016); von
Fellenberg et al. (2018); Witzel et al. (2018); Bower et al. (2019); Gravity
Collaboration et al. (2020c). The (dereddened) NIR points represent different
attempts to quantify the mode (peak) of the NIR flux-density distribution
(compare discussion in Witzel et al. 2018). At the Herschel bands (100, 160,
and 250 $\micron$) two numbers are presented: the amplitude of the detected
variable signal (blue lower limit, grey points) that serve as an lower limit
of the steady flux density, and for 100 and 160 $\micron$ upper limits of the
steady flux density derived from the assumption of 15% variability (von
Fellenberg et al. 2018).

Another, potentially related, problem is the submm variability variance: at
230 GHz the variance at longest timescales cannot be reproduced by the model.
This is the reason we did not include the last bin of the 230 GHz structure
function in our distance function. Furthermore, Figure 1 shows significant
offsets in average flux density between the individual SMA and ALMA epochs. If
real, these offsets would point to a long term variability component that is
not accounted for in our model. Several explanations are possible:

  * •

We modeled the synchrotron source as a sphere homogeneously filled with
electrons, resulting in an optically thick spectral index
$\alpha_{\rm{thick}}=+2.5$. In the case of a different electron density
profile, this slope can be significantly different with consequences for the
position of the self-absorption turnover and the relative flux densities
between the 230 GHz, 345 GHz, and the NIR.

  * •

A separate mechanism, perhaps a second electron population, is required to
produce the non-varying emission, which accounts for most of the submm and
longer-wavelength flux most of the time. This second (or rather primary)
component is not needed to account for most of the observed submm variability,
but it might contribute, especially on long timescales. This would reduce the
variability associated with the compact source, resulting in a flatter
optically thin spectral index.

  * •

Instead of one zone undergoing cycles, there might be distinct zones coming
and going with different initial conditions, explaining the distinct levels of
the submm epochs.

Figure 22 compares the SED of the compact component modeled here with the so-
called steady SED of (selected) literature values. A minimum requirement for
our model should be that it does not violate the steady SED in the submm or
radio. For a source that is variable in all parts of the spectrum, it is a
matter of definition what can be considered the steady flux. While in the
radio the variability amplitude is a small fraction of the mean flux, this
problem is prominent NIR, where the source is entirely dominated by variable
flux, and no truly steady point source has been identified. Additionally,
observations of the mean flux (especially in the presence of variability) are
limited by their duration, their resolution, extended flux surrounding the
Galactic center, and the atmosphere. Despite the difficulties, we see a
certain degree of scatter at radio to submm flux densities which most likely
are caused by intrinsic variability. At frequencies above submm, the only
detections of Sgr A* are a far-infrared measurements with Herschel (Stone et
al. 2016; von Fellenberg et al. 2018) and NIR and X-ray measurements as
presented here. Our model can describe the variability amplitudes at NIR and
FIR frequencies well and even follows the spectral index suggested by the
observed data in the FIR while not violating any upper limits or the radio or
submm part of the SED. In the submm, our model significantly contributes to
the overall variability but cannot explain the entire flux density. At lower
radio frequencies, the variable contribution decreases quickly, and the SED is
dominated by one or more other source components, which might contribute to
the excess submm variability at longer timescales that we cannot account for.

In synchrotron–SSC models, the X-rays have the same spectral index as the
optically thin part of the synchrotron spectrum. In their careful analysis of
the X-ray spectral slope, Ponti et al. (2017) found a photon index
$\Gamma=2.27\pm 0.12$, which corresponds to a spectral index of
$\alpha_{X}=1.27\pm 0.12$, about 1$\sigma$ consistent with our $\alpha=0.98\pm
0.1$. Haggard et al. (2019) analyzed the two brightest X-ray flares detected
so far and find $\Gamma=2.06\pm 0.14$ and $\Gamma=2.03\pm 0.27$, respectively,
again consistent with our results.

Finally, our model gives the brightest observed NIR flare so far (Do et al.,
2019) only ${\sim}10^{-3}$ probability. In other words, the brightest observed
NIR states of Sgr A* are somewhat underrepresented in our mock light curves.
Dodds-Eden et al. (2011) and Gravity Collaboration et al. (2020c) interpreted
the NIR flux-density distribution as the result of two states with variability
from two distinct physical processes. Whether this is the case or instead very
bright events are only extreme cases of the process described here cannot be
determined. In order to generate extreme events more regularly, it would
suffice to replace log-normality of the underlying fast and slow processes
with a suited distribution. However, proper inclusion of extreme values in the
statistics is difficult because the probability of rare events is very
uncertain.

##  7 Summary

This paper has presented a comprehensive analysis of submm, NIR, and X-ray
light curves of Sgr A* including two new epochs of ALMA data and four new
epochs of SMA data. These include about 2 days of simultaneously observed NIR
and submm light curves, which show definite correlation between submm and NIR
variability and are consistent with a median delay of $\sim$20 minutes of the
submm with respect to the NIR.

This paper is also the first analysis of the X-ray PSD. The X-ray and NIR
variations are correlated but with correlation decreasing at short variability
timescales.

A simple but physically consistent model explains most features of Sgr A*’s
variability. The model was constructed to match the variability amplitudes and
timing properties (structure functions) of the NIR and X-ray emission
considered separately. Given those, the model predicts the submm variability
properties, the time lag between NIR and submm variations, the submm and NIR
spectral indices, and the observed correlation between NIR and X-ray
variability. The model consists of a single zone in an external magnetic field
with random injections of high energy electrons and cyclic expansion and
contraction. The radio emission is optically thick synchrotron, the NIR is
optically thin synchrotron, and X-rays arise from synchrotron self-Compton
emission. The variability comes from varying density of high-energy electrons
and source size, which affect the amplitude of synchrotron emission, the
optical depth, the synchrotron cutoff frequency, and the self-Compton
amplitude. The minimum required magnetic field is $\sim$8.5 G, the maximum
source size is $\sim$2.7 $R_{S}$, and the high-energy electrons have a power-
law index $\sim$2.95. The variability processes have timescales of roughly 82
and 135 minutes and are ruled by the tradeoff between electron injection and
expansion cooling in the submm and X-rays and by the tradeoff between
injection and synchrotron cooling in the NIR. General relativity plays no role
in the model.

The predicted changes in source parameters are moderate, typically about a
factor two in source size and magnetic flux density and an order of magnitude
in electron density. The submm radiation of this compact component of Sgr A*
changes from optically thick to thin and back at the variability timescales.
This likely creates complex changes in the intrinsic source structure during
VLBI observations.

Postulates of the model include the PSD of the underlying trigger processes
and that rises and falls are symmetric in time. The model allows for but does
not require different correlation times at different frequencies and time lags
between frequencies. Where the model may fall short is that the most extreme
observations are improbable (but not impossible), and the observed NIR
spectral index at high flux densities is flatter than the model predicts. That
last could probably be remedied by allowing SSC emission at NIR frequencies.

All in all, the model’s shortcomings are related to rare observations such as
the brightest flux-density states and therefore are documented insufficiently
in our sample to warrant a more complex model. For the vast majority of
observations, the model presented here is a sufficient representation. At the
least, it provides a baseline to help identify extraordinary variability,
e.g., states dominated by effects of general relativity in the future.
Relativistic effects alone cannot be the origin of the variability because
they cannot explain the rich phenomenology of correlations between the
wavelengths. Our model in its current form does not require any relativistic
boosting or other effects to describe the variable flux densities. Whether and
how this is commensurable with the findings of apparent circular motion close
to the last stable orbit by Gravity Collaboration et al. (2018b) has to be
investigated.

Physical conclusions from this work include:

  * •

There is a compact component with size 2–$3R_{S}$ that dominates the high
frequency regime.

  * •

There is more evidence now for SSC generating the X-rays (and maybe
contributing to the NIR).

  * •

The NIR is linked to the submm variability, and expansion is a strong
candidate process to explain the correlation.

Some things the model does not explain include:

  * •

Long-timescale variability in the submm.

  * •

Perhaps the most extreme NIR flux densities.

This paper has presented a quantitative assessment of what a simple model can
achieve. It can serve as a starting point for further work to test and improve
the model. Theoretical work should investigate acceleration mechanisms for
electrons and whether suitable populations of high-energy electrons can be
generated and on what timescales.

The model presented here was enabled and its parameters determined by the vast
library of light-curve monitoring data accumulated over the past two decades.
Simultaneous light curves at different wavelengths were especially valuable
for testing (though not deriving) the model. Further work should include
better simultaneous measurements of the X-ray and NIR spectral indices, which
should be identical if our model is correct. Monitoring Sgr A*’s variability