The Relationship Between Simulated Sub-Millimeter and Near-Infrared Images of Sagittarius A* from a Magnetically Arrested Black Hole Accretion Flow

Sagittarius A* (Sgr A*), the supermassive black hole at the center of the Milky Way, undergoes large-amplitude near-infrared (NIR) flares that can coincide with the continuous rotation of the NIR emission region. One promising explanation for this observed NIR behavior is a magnetic flux eruption, which occurs in three-dimensional General Relativistic Magneto-Hydrodynamic (3D GRMHD) simulations of magnetically arrested accretion flows. After running two-temperature 3D GRMHD simulations, where the electron temperature is evolved self-consistently along with the gas temperature, it is possible to calculate ray-traced images of the synchotron emission from thermal electrons in the accretion flow. Changes in the gas dominated ($\sigma=b^2/2\rho<1$) regions of the accretion flow during a magnetic flux eruption reproduce the NIR flaring and NIR emission region rotation of Sgr A* with durations consistent with observation. In this paper, we demonstrate that these models also predict that large (1.5x - 2x) size increases of the sub-millimeter (sub-mm) and millimeter (mm) emission region follow most NIR flares by 20 - 50 minutes. These size increases occur across a wide parameter space of black hole spin ($a=0.3,0.5,-0.5,0.9375$) and initial tilt angle between the accretion flow and black hole spin axes $\theta_0$ ($\theta_0=0^{\circ}$, $16^{\circ}$, $30^{\circ}$). We also calculate the sub-mm polarization angle rotation and the shift of the sub-mm spectral index from zero to -0.8 during a prominent NIR flare in our high spin ($a=0.9375$) simulation. We show that, during a magnetic flux eruption, a large ($\sim10r_g$), magnetically dominated $(\sigma>1)$, low density, and high temperature ``bubble'' forms in the accretion flow. The drop in density inside the bubble and additional electron heating in accretion flow between 15$r_g$ - 25$r_g$ leads to a sub-mm size increase in corresponding images.


INTRODUCTION
Sagittarius A* (Sgr A*), the object at the center of the Milky Way, is a bright radio source (Balick & Brown 1974;Ekers et al. 1975;Lo et al. 1975) with variable near-infrared (NIR) emission reaching a factor of 10x -25x (Genzel et al. 2003b;Ghez 2004;GRAVITY Collaboration et al. 2020b) its median value in flux (Schodel et al. 2011;Dodds-Eden 2011;Witzel 2018).Sgr A* also exhibits X-ray flares reaching 50x the quiescent value (Markoff et al. 2001) which often coincide with a NIR flare (Baganoff 2001;Genzel et al. 2003a;Ghez 2004;Eckart 2008;Ponti et al. 2017;GRAVITY Collaboration et al. 2021).Astrometry of stars orbiting Sgr A* (Ghez 2008;Genzel et al. 1997) and the resolved black hole shadow at 230 GHz captured by the Event Horizon Collaboration (EHT) (Event Horizon Telescope Collaboration et al. 2022a) provide consistent and strong evidence that Sgr A* is a supermassive black hole with mass ∼ 4 × 10 6 ⊙.Linear polarization fractions of 10% -40% in the NIR and 230 GHz indicate that the emission region around the black hole is dominated ★ E-mail: arpi.grigorian@colorado.edu by synchrotron radiation (Eckart 2008;Trippe et al. 2007;Wielgus et al. 2022).The lack of sub-mm/mm variability on the same order as the NIR variability suggests that the origin of the flares are electrons accelerated to extreme velocities ( ≳ 10 3 ) (Markoff et al. 2001;Dodds-Eden 2009).For example, the magnetic Rayleigh-Taylor instability stemming from the interface of a dense plasma and a magnetic cavity can drive efficient particle acceleration (Zhdankin et al. 2023).Weaker flares may be also be caused by the lensing of turbulent accretion flow (Chan et al. 2015).However, the origin of the NIR flares, especially stronger ones, remains unclear.mm polarization angle (Marrone et al. 2006).Using the ALMA Observatory, Wielgus et al. (2022) recently found an example of the 230 GHz polarization angle that makes a quasi-circular path in Stokes U, Q (similar to the path of NIR polarization angle) over a period of 74 ± 6 minutes following an X-ray flare.
The GRAVITY observations are well described by a hot-spot at ≃ 6  −10  orbiting around the black hole at a low ( ≤ 30 • ) inclination angle and polodial field (GRAVITY Collaboration et al. 2018;et al. 2020).However, these models do not address possible physical mechanisms that could change the gas and electron temperature in the inner accretion flow.
Three-dimensional General Relativistic Magneto-Hydrodynamic (GRMHD) simulations of a black hole surrounded by an accretion disk or flow offer insights into the physical mechanisms that may be present in Sgr A*.Two key GRMHD black hole accretion disk models often explored in the literature are Standard and Normal Evolution (SANE) disks and Magnetically Arrested Disks (MAD).SANE disks, which do not reach high levels of magnetic flux threading the horizon (Narayan et al. 2012), likely rely on the Magnetorotational Instability (MRI) for angular momentum transport (Hawley & Balbus 1991;Balbus & Hawley 1998).MADs, by contrast, accrete large amounts of magnetic flux that ultimately saturates (Bisnovatyi-Kogan & Ruzmaikin 1974, 1976;Igumenshchev et al. 2003;Begelman et al. 2022).
In this paper, we study similar MAD simulations as Dexter et al. (2020b) over a wider parameter space.Specifically, we vary the black hole spin  and consider simulations that begin with an initial tilt angle  0 between the black hole spin axis and the accretion flow spin axis.In Section 2 (Methods), we describe the GRMHD set-up, how images are created using outputted data from our simulations, and how we extract key observables from those images.Then, we focus on an example from the high spin ( = 0.9375), prograde ( 0 = 0 • ) simulation where the dissipation of magnetic flux near the horizon creates a magnetically dominated, low density "bubble" (hot spot) within 10   of the accretion flow.The bubble itself is a high  =  2 /2 > 1 region, where  2 /2 is the magnetic pressure and  is the fluid density.However, the bubble formation and orbit around the black hole is coincident with significant electron temperature increases in regions where the plasma parameter  < 1.When the bubble eventually reaches larger radii (15   -30   ), the electron temperature also increases in this region.We show these changes in the  < 1 regions of the inner accretion flow not only power a NIR flare and quasi-circular NIR centroid rotation (similar to results from Dexter et al. (2020b)), but also cause a significant increase in both the sub-mm and mm emission regions about 1 hour after the NIR flare peak.This example produces the rotation of the sub-mm polarization angle with a period of ∼2 hours, consistent with observations (Marrone et al. 2006;Wielgus et al. 2022).The luminosity spectrum evolves during the course of the NIR flare, and the peak of the sub-mm size increase is simultaneous with a decrease in the synchotron spectrum peak frequency   and, as a result, a decrease in the spectral index near 230 GHz (1.3 mm).
In Section 7.1, we demonstrate that the relationship between NIR flares and the sub-mm/mm size increases is a feature of all simulations in our parameters space.The relationship between the NIR total flux and the sub-mm/mm size increases is stronger and more consistent that the relationship between the NIR total flux and the sub-mm total flux increases.We also present the NIR centroid behavior during our simulations.Finally, we discuss limitations of the GRMHD model, the dependence of size increases on inclination angle, and the implications our results have on illuminating the physical mechanisms behind NIR flares in Sgr A*.

The GRMHD MAD Model
We carried out 3D GRMHD simulations for a range of black hole spin  and initial tilt angle of the accretion flow  0 (See Table 1 for the full parameter space).The dimensionless black hole spin is defined as  ≡ /(  2 /), where  is the angular momentum and  is the black hole mass.We refer to simulations with || > 0.5 as high-spin; simulations with || ≤ 0.5 are low-spin.The angle  0 is measured between the spin axis of the black hole and the initial angular momentum vector of the accretion flow.Simulations with  0 = 0 • are prograde (the black hole spin and the accretion flow spin are aligned).We use the convention of negative spin values for simulations where the black hole spin and the accretion flow spin are completely anti-aligned.
Simulations were carried out using the public harmpi code (Tchekhovskoy 2019), a 3D version of HARM (Gammie et al. 2003;Noble et al. 2006).All simulations studied here were initialized from a hydrodynamic equilibrium torus with an inner radius of 12  and a pressure maximum at 25  .The (arbitrary) initial magnetic field was composed of a single poloidal loop so that max(   )/max(   ) = 100, where   ,   are the gas and magnetic pressures, respectively (See Tchekhovskoy et al. (2011)).In order to study MAD conditions, the initial magnetic field was concentrated further out in the torus so that magnetic flux could accumulate and saturate on the black hole.For more details, see Dexter et al. (2020b).
We use Kerr-Schild coordinates and a grid resolution of 320x256x160 in the , , and  directions, respectively.The  values are spaced log-normally to resolve the important physics occurring near the horizon.The  coordinate is irregular, chosen to concentrate resolution near the midplane at small radius.The metric in all simulations is fixed as the stationary Kerr metric determined by the black hole mass and spin.We used a temporal cadence of 10 GM/c 3 when plotting the fluid variables and observables calculated from ray-traced images as functions of time.
We also include a scheme in our simulations that self-consistently evolves the electron internal energy density for four different electron heating prescriptions along with a single, MHD fluid (Ressler et al. 2015).For this work, we choose the electron energy density calculated from the heating prescription developed by Werner et al. (2018) from their study of particle acceleration from magnetic reconnection using 2D particle-in-cell (PIC) simulations (W18).They found an empirical relationship between the electron heating fraction   and the ion magnetization   ≡  2 /4     2 , where  is the upstream magnetic field, and   is the ion density upstream of the reconnection region (see their Equation 3).The electron heating fraction ranges from 1/4 for low magnetization to 1/2 for high magnetization.W18 was able to reproduce a variety of observed signatures for Sagittarius A* in previous studies (Dexter et al. 2020b;Dexter et al. 2020a).By applying W18 throughout the entire coordinate space of the simulation, our model assumes that 100% of the dissipated heat is from reconnection.
We evolve our simulations out to times ranging from 20,000 GM/c 3 to 100,000 GM/c 3 (See Table 1 for the exact duration for each simulation in our parameter space).The durations are long enough for the fluid in inner (< 50  ) region to reach equilibrium.For all the  = 0.9375 simulations and the  = 0.5 simulation with  0 = 30 • , we ignored the first 10,000 GM/c 3 to allow the simulation to evolve from the initial configuration of the fluid to a steady state.However, the remaining low spin  ≤ 0.5 simulations took longer than 10,000 GM/c 3 to start exhibiting consistent magnetic flux eruptions.
Our goal is to compare simulations when they are clearly in a MAD state.Therefore, the initial period of time that we ignored for these simulations was longer and ranged from 16,000 GM/c 3 to 26,000 GM/c 3 (See Table 1).However, because the simulation durations were long enough, for each simulation there was a significant (> 16,000 GM/c 3 ) continuous period time containing magnetic flux eruptions that we could compare across our parameter space.

Synthetic black hole images
We created monochromatic images for each simulation in Table 1 using the the ray-tracing radiative transfer code grtrans (Dexter et al. 2009;Dexter 2016).The code calculates all four Stokes parameters for each image: total intensity I, linear polarization Q and U, and circular polarization V.The grtrans code includes all relativistic effects in a Kerr space time.
We calculate the synchrotron emission with grtrans assuming a thermal distribution of electrons.We do not include emission from non-thermal electrons (see the Discussion for how non-thermal electrons may effect the sub-mm, mm, and NIR emission.)We set the inclination angle , the angle between the black hole spin axis and the line of sight, to  = 25 • (unless otherwise noted).Since the angular momentum of the black hole is along positive  (spinning counterclockwise) in our simulation convention,  < 90 • means that the black hole, and, generally, the accretion flow, are spinning counter clockwise on the sky.
When calculating images, we omit regions of the accretion flow which have a  value greater than -cut.Unless specified otherwise, we used -cut = 1.For our fiducial example, we also produced images with -cut = 10 (See Section 8.1 in the Discussion).
We made images at 3.5 mm (86 GHz), 1.3 mm (230 GHz), and 2.2m (136 THz) for all simulations in this work.We refer to all 3.5 mm images as mm (millimeter), 1.3 mm images as sub-mm (sub-millimeter), and 2.2m images as NIR (near-infrared).Most images in this study use the following fields of view as a function of frequency: 150 2 as 2 for the NIR (2.2m) images, 300 2 as 2 for the sub-mm (1.3 mm) images, and 400 2 as 2 for the mm (3.5 mm) images.The black hole is always centered in the field of view.With the exception of images used to calculate the luminosity spectrum, images are 300x300 pixels in size.
In order to report distances in terms of subtended angles on the sky (as) and times in units of minutes/hours/days, we use 8 kpc for the distance to the Galactic Center and  = 4 × 10 6  ⊙ for the mass of the black hole.For reference, 1,000 GM/c 3 is approximately 5.6 hrs.
In GRMHD simulations, the mass accretion rate across the horizon is calculated self-consistently.However, in MAD simulations, matter is accreted in clumps and this "true " mass accretion rate is noisy.Instead of using the "true" mass accretion rate to scale our images, we choose a constant value for the mass accretion rate mdot every 500 frames so that the average sub-mm total flux is near the observed value of ∼3 Jy for Sgr A*.When calculating images at higher -cut = 10, we use a different value of mdot so that the sub-mm images for each -cut meet this observational total flux constraint.
The only exception to the rule above for scaling images is set of -cut = 1 NIR images from the high-spin, prograde simulation between 10,000 GM/c 3 and 57,190 GM/c 3 .For this case, we first found a value of mdot that recreates the observed Sgr A* sub-mm total flux for the first 500 frames and calculated all the NIR images during this period of time with the same choice of mdot.Then, we calculated a "smoothed" curve for the mass accretion rate ⟨ ⟩ as a function of time by taking the actual, intrinsically noisy mass accretion rate across the black hole horizon  and smoothing it with a Gaussian kernel.We choose a standard deviation of the Gaussian kernel that produces a ⟨ ⟩ which is monotonically decreasing so that any peaks in  do not erroneously create peaks in the NIR total flux curve (we also use this procedure to scale the dimensionless magnetic flux on the black hole horizon, see Section 2.7).We then scaled the NIR images by ⟨ ⟩ 2 .The NIR images from the high-spin, prograde simulation after time 57,190 GM/c 3 were calculated using the former method.

Calculating Total Flux, Centroid, and Image Size
We calculate the light curves, sizes, and centroid positions directly from the grtrans images.
The moment of an image  , with image coordinates ,  is defined as (1) The total flux   = ∫   cos()Ω from a far-away source can be rewritten as where  is the distance to the Galactic Center and ,  are coordinates on the sky.In this coordinate system, the black hole is centered in the image.Discretizing the integral yields where ,  is the total height and width of the field of view in units of the gravitational radius   = GM/c 2 and N is the total number of pixels along the  and  directions.
The centroid position on the sky  (, ) is given by  (, ) = ( 10 / 00 ,  10 / 00 ) = ( Ī , Ī ). (4) To find the size, or the spread of the intensity map, for each image, we first solve for the eigenvalues of the following matrix (5) which gives the major (   ) and minor (    ) axes of an ellipse that corresponds to the "spread" of the intensity on the image.We approximate the emission region as Gaussian and find the maximum and minimum Full-Width Half-Max (     ,       , respectively) by using the simple relation

The Cross Correlation Coefficient
Each of the observables defined above at a single frequency, such as the total flux, Ī , Ī ,      ,       , produces a curve that is a discretized function of time with equal time steps: a time series.The cross-correlation of two time series is a convolution of one curve with the other that measures their degree of similarity as well as any characteristic lags between them.We use the following definition of a discrete cross-correlation coefficient CC The cross correlation CC [ (), ()] () is the functional that takes two curves (), () and outputs a function that depends on the separation time delay  between the two curves.The quantities ⟨ , ⟩ and ⟨ , ⟩ each correspond to the mean value to their respective functions , in the interval ( , ).The maximum of CC() shows the strength of the correlation between two curves and the location of the peak value is a characteristic delay time  * .
Cross correlation values greater than 0.5 are strong, between 0.25 and 0.5 are weak, and values less than 0.25 are practically uncorrelated.
We can use the cross correlation function to define the auto correlation for a single time series (): CC [ (), ()] ().The autocorrelation function measures the periodicity of a time series; a strong auto-correlation coefficient at a non zero  * suggests the curve contains a mode that has frequency 1/ * .

The Structure Function
The first order structure function for a finite, continuous time series A () is defined as (Simonetti et al. 1985;Hughes et al. 1992) where  is the total duration for the time series.We then find the discrete structure function where  is the total number of points in our discrete time series  (  ), and apply it to our discrete time series that represent simulated observables.
The structure function ignores any DC offsets present in a time series (Simonetti et al. 1985;Dexter et al. 2022;Witzel 2018).Therefore, we can compare the simulated structure function to observation even if there is a source of steady emission in Sagittarius A* or background emission that is not included in our model.For curves that are well modeled by a Damped Random Walk (DRW), the structure function reaches an asymptotic value  ∞ at late times (MacLeod et al. 2010).If the structure function becomes close to  ∞ at finite time   , this time represents when the series becomes uncorrelated or is transitioning from red noise to white noise.We choose   as the time when the structure function  reaches its first prominent peak before plateauing.

Computing Average Values of Fluid Variables
In this work, most of our plots of fluid variables show averages of that variable over the polar angle  or the azimuthal angle .For some function  (, , ) that represents a scalar field of the fluid, the polar and azimuthally averaged values ⟨  (, )⟩ and ⟨  (, )⟩ are, respectively Sometimes it is more instructive to show the weighted average of a fluid variable with respect to the density .The weighted averages over  and  are given, respectively, as

Calculating the Dimensionless Magnetic Flux on the Event Horizon
As the simulation evolves, magnetized material accretes onto the black hole.This process allows the magnetic flux threading the surface just outside the event horizon to accumulate.The total magnetic flux Φ at over the spherical shell at radius  is defined by where   is the radial component of the three-vector magnetic field and √ − is the determinant of the metric.
Instead of Φ, we use the dimensionless magnetic flux   given by We calculate smoothed mass accretion rate ⟨ ⟩ by taking the actual, intrinsically noisy mass accretion rate across the black hole horizon  and smoothing it with a Gaussian kernel.We set the standard deviation of the Gaussian kernel so ⟨ ⟩ is monotonically decreasing; this step ensures that the variability of  does not translate to the dimensionless magnetic flux curve (this is similar to the procedure used to scale the total flux of some of our NIR images).Therefore, all the fluctuations seen in   are due to changes in the magnetic flux rather than due to variability in the mass accretion rate.

MAGNETIC FLUX ERUPTIONS IN THE MAGNETICALLY ARRESTED ACCRETION FLOW
For all simulations studied, the magnetic flux on the horizon   enter periods of slow build up followed by rapid dissipation.Most drops in   , called magnetic flux eruptions, are also accompanied by significant increases in the NIR total flux: a NIR flare.
Figure 1 shows that   cycles many times during the high-spin, prograde simulation.Most significant drops in magnetic flux are also accompanied by NIR flares.We find the average time between magnetic flux eruptions to be roughly ∼1,600 GM/c 3 (∼9 hours), although the low (0.2) auto-correlation peak of the   curve suggests no consistent period of recurrence for the magnetic flux eruptions.
Figure 2 shows the evolution of several fluid variables during the course of a magnetic flux eruption taken from the high-spin, prograde simulation (all times occur within the time interval outlined by a black box in Figure 1).
Shortly after the magnetic flux peak at 37,230 GM/c 3 , a low density, high temperature "bubble" begins to form in the accretion flow.At 37,290 GM/ 3 , the NIR total flux curve reaches the most prominent peak during this magnetic flux eruption event.The bubble is now ∼ 10  in size and the electron temperature in the bubble is 100x -1000x the median value (second column in Figure 2).The fluid parameter  =  2 /2 has also increased above and below the equatorial plane; only a thin region within ∼ 7  of the accretion flow remains mostly  < 1.
As the bubble rotates around the black hole, it continues to grow.By 37,360 GM/c 3 , the bubble is ∼ 15  in size (third column in Figure 2).At this time, the innermost (≲ 10  ) region of the accretion flow is mostly magnetically dominated.The sub-mm size is also increasing (See Figures 3 and 4).
Around 37,400 GM/c 3 , the low density region becomes elongated; it stretches to radii ∼ 20  .Density plots of the slice along the equatorial plane show that a single, low density region breaks apart into multiple low density regions.By 37,460 GM/c 3 , the slice of the equatorial plane also shows that the low density region within ∼ 10  remains magnetically dominated but the  value of the low density region at larger radii has decreased to less than one.Also, at 37,460 GM/c 3 , the electron temperature has increased to 10x the median value in the 15  -25  region of the accretion flow.The temperature within ∼ 10  remains elevated (last column of Figure 2).The sub-mm size at 37,460 GM/c 3 is at a maximum (See Figure 4 for the sub-mm image corresponding to this time).Throughout all the events described above, the magnetic flux has been monotonically decreasing.The magnetic flux reaches a minimum at 37,540 GM/c 3 .By this time, the outer region of the bubble has reached 30   .The sub-mm size is still elevated from its median value but is decreasing.
The total time between the magnetic flux peak and its local minimum is 310 GM/c 3 , or 1.7 hours.

NIR TOTAL FLUX AND SUB-MM/MM TOTAL FLUX AND SIZE DURING A MAGNETIC FLUX ERUPTION
During the course of our fiducial magnetic flux eruption explored above (Section 3), the NIR flux significantly increases or "flares".It reaches a peak value of 22.2 mJy, or 10x of the median NIR flux, at 37,290 GM/c 3 (Figure 3).After the NIR flux peak, the NIR total flux stays elevated and reaches several lesser peaks over the course of roughly 1 hour before quickly returning to its median value.This NIR flare is also accompanied by significant (∼2x) increases in both the sub-mm and mm       and      .Both the sub-mm and mm size reach a peak value about 1 hour after the NIR total flux peak.The mm      peak lags behind the submm      peak by 10 minutes.Near the peak value of the sub-mm and mm size increases, the NIR total flux curve quickly drops back to its median value.The sub-mm/mm sizes, however, slowly return to their quiescent values over the course of ∼30 minutes.
Before the NIR flare, the sub-mm total flux shows significant variability between 2 Jy to 5 Jy.During the NIR flare, the sub-mm flux climbs again to 5 Jy.Therefore, the peak in sub-mm flux is not unique to periods with NIR flaring.
After the second 5 Jy peak, the sub-mm flux quickly drops below 2 Jy.The minimum in the sub-mm total flux curve is simultaneous with the sub-mm      peak as well as a minor peak in the NIR total flux curve (See dotted vertical blue line in Figure 3).
The mm total flux, both during the quiescent and flaring periods, shows less variability than its sub-mm counterpart and stays near 2.5 Jy.
Figure 3 also shows a minor NIR flare near 36,800 GM/c 3 , which is coincident with a minor magnetic flux eruption event.The NIR flux during the flare only reaches a peak value ∼6 mJy and stays elevated for 80 minutes.Also, unlike the prominent NIR flare, there are no significant sub-mm/mm size increases within 1 hour after the minor NIR flare.The low density "bubble" that forms during the minor flare is smaller in size than its counterpart from the prominent flare, demonstrating that the geometry of the bubble has an impact on the resulting sub-mm/mm size.

THE SUB-MM AND NIR CENTROID MOTION DURING
A MAGNETIC FLUX ERUPTION

The NIR Centroid
Prior to the prominent magnetic flux eruption discussed above from the high-spin, prograde simulation, the NIR centroid curves stay near the median value of (-6.5 as, 0.3 as).The NIR centroid position does not change significantly during the minor NIR flare (Figure 3).However, during the second, prominent NIR flare, the NIR centroid makes two orbits around the black hole over a period of ∼ 3 hours, tracing a quasi-circular path with a peak diameter of 46 as (See the right side of Figure 4 and Figure 3).The NIR centroid returns to its quiescent position at the end of the NIR flare.

Sub-mm Centroid Motion and Linear Polarization
During the prominent NIR example flare, the sub-mm centroid makes a similar, quasi-circular orbit on the sky as the NIR centroid over a period of ∼ 3 hours (Figure 5).The sub-mm centroid position is also less noisy than its NIR counterpart.Considering both centroids follow the counter-clockwise motion of the accretion flow, the submm centroid is ahead of the NIR centroid by ∼ 20 minutes.The effective diameter of the sub-mm centroid path, which reaches an extreme at 28 as, is also smaller than its NIR counterpart.
During the 5 hours prior to the flare, we find that the sub-mm linear polarization integrated over the image behaves like a random walk in Stokes U and Q space with an average linear polarization fraction ( = (|U| 2 + |Q| 2 )/|I| 2 ) of 4.5%.However, right as the submm centroid begins its quasi-circular path, the sub-mm polarization vector changes behavior and slowly rotates in Stokes U and Q.For the first ∼2 hours, shown in Figure 6, the path traced in Stokes U and Q space is a clear loop; the shape of the path in Stokes U and Q looses this character for the final ∼ 1 hour of the sub-mm centroid motion.While the polarization vector angle slowly rotates, the total linear polarization fraction reaches a maximum value of 11%.

THE LUMINOSITY SPECTRUM AND SUB-MM SPECTRAL INDEX DURING A MAGNETIC FLUX ERUPTION
We calculate the spectrum () directly from the intensity maps by finding the total flux   over a wide range of frequencies .Then, we find   = 4 2   , where  is the distance from the galactic center to the observer.
The spectral index is defined as  ≡ − log( ())/ log .We calculate the spectral index near 230 GHz (1.3 mm),   , by finding log( (log())) for five frequencies between 210 GHz and 250 GHz and then fitting the function to a line using numpy.polyfit.
Figure 7 shows the resulting spectra.The spectrum at  = 36,200 GM/c 3 , when the accretion flow is an a quiescent state, has a shape similar to the power spectrum of a single electron and   = -0.07.At the NIR flare peak ( = 37,290 GM/c 3 ), the high frequency tail of  () is significantly elevated with only a slight increase in luminosity for sub-millimeter and millimeter frequencies but   remains relatively unchanged with a value of 0.06.There is also little to no change in the total luminosity for millimeter and radio frequencies.
When the simulated images of the accretion flow undergo sub-mm and mm size increases, the luminosity spectrum continues to evolve.At ( = 37,460 GM/c 3 ), when the sub-mm size reaches its peak value, the spectrum becomes a combination of an elevated high frequency The density  (in code units), the temperature in units of /   2 weighted by the density , and the plasma  parameter 2 / 2 , where  is the fluid pressure and  2 /2 is the magnetic pressure, all averaged over polar angle .We also show the plasma  parameter  2 /2 averaged over azimuthal angle  (see Equations 10 -13 for how we compute averages).Note that the outermost radius shown in the  parameter plots is less than the outermost radius shown for the other fluid parameter plots.A low-density, high temperature bubble forms in the innermost accretion flow and breaks apart before moving to larger radii during the 2,000 GM/c 3 interval outlined by a black box in Figure 1.In the far left column, the accretion flow is in a quiescent state.The NIR total flux is near its median value and the sub-mm size is ∼50 as.The temperature and density is high close to the event horizon at ∼ 2  and then falls off at larger radii.Spiral structure in the density plot is visible.At  = 37,290 GM/c 3 , a large scale (10   ) low density, high temperature, and low plasma  region forms.Also, at this time, the magnetic flux is expelled from the horizon and NIR total flux curve is at a maximum.At  = 37,360 GM/c 3 the low density bubble has increased in size and pushed material out to larger radii.The temperature in this region stays elevated and additional smaller, even higher temperature spots appear throughout the bubble.In the rightmost column, at  = 37,450 GM/c 3 , the sub-mm emission region size is at a maximum.The bubble is dissipating.The temperature in the bubble has started to cool and the density of the bubble has increased.There is a small (10x) increase in temperature in the 15   -25   region.Following the trajectory of the bubble in the last three columns, one can see that it rotates counter-clockwise with the general accretion flow.Prior to the second flare, both the NIR centroid (middle subplot) and the sub-mm/mm size (lower subplot) undergo only small deviations from their median values.However, during the second flare, the NIR centroid changes its behavior and makes a quasi-circular orbit around the black hole.The sub-mm and mm sizes begin to climb.The sub-mm peak size increase at 37,460 GM/c 3 coincides with a final, lesser NIR total flux peak during the second flare before the total flux quickly drops to its median value.
The mm      reaches at maximum only 10 minutes after the sub-mm      peak and outside the flaring period.Then, the sub-mm and mm sizes slowly decrease back to their quiescent values.The entire event, from the beginning of the first flare to the end of the sub-mm/mm size increases, lasts 800 GM/c 3 , or just under 4.5 hours.However, the separation between the peak value of the prominent NIR flare and the sub-mm      peak is only 57 minutes (270 GM/c 3 ).The width of the sub-mm size increase peak is less than 2 hours.We also include the observed sub-mm size from the EHT (Event Horizon Telescope Collaboration et al. 2022a) and the observed mm sizes from ALMA (Issaoun et al. 2021) in the bottom panel of the plot.The light blue region around the blue EHT line represents the 68% credible interval; the light grey region around the observed mm      /    lines represents the 95% confidence interval for each measurement.We discuss the agreement between the sub-mm size and our model, as well as the tension between our calculated mm sizes and the observed values, in Section 8.6 of the Discussion.
tail and a synchrotron spectrum with a lower peak frequency   than at the earlier, quiescent time.The shift in the peak synchrotron frequency is clearer at time  = 37510 GM/c 3 , when the sub-mm size is still 1.6x times the median value but the NIR flare has dimmed.
This shift in the synchrotron thermal spectrum at  = 37,460 GM/c 3 and  = 37,510 GM/c 3 causes   to become negative compared to the quiescent and peak NIR total flux times.We calculate   = -0.84 at the peak sub-mm size and   = -0.49at  = 37,510 GM/c 3 .

NIR Total Flux and Sub-mm Size
By running the high-spin, prograde simulation to 96,310 GM/c 3 , we were able to characterize the long-term behavior of the sub-mm and NIR total flux, centroid motion, and size.The auto-correlation functions for both the NIR total flux and the sub-mm     (Figure 8) show no prominent peaks besides  = 0, demonstrating a lack of periodicity.The structure functions  for the NIR total flux, sub-mm total flux, and the sub-mm      (Figure 9) all reach an asymptotic value   near 1,000 GM/c 3 (5 hours for Sgr A*).Since the NIR total flux curve   is less than the time between  magnetic flux eruptions, the NIR flares are likely independent from one another.Significant increases in the sub-mm      follow many NIR flares during the high-spin, prograde simulation run (Figure 10).This relationship, however, is not unique to our fiducial simulation and persists across our entire parameter space.Figure 11 shows a positive example of a sub-mm size increase following a NIR flare for every simulation studied.The cross-correlation coefficient between the NIR total flux and the sub-mm      for every simulation has a single, significant peak with a value near 0.6, a positive characteristic time delay  * < 1 hour, and peak width of ∼3 hours (Figure 12).In our convention, a positive  * means that the NIR total flux curve's behavior precedes the sub-mm      curve's behaviour.
However, we do find a dependence between the black hole spin and accretion flow orientation on  * : low-spin and low initial tilt angle ( 0 ) simulations produce the shortest  * ∼ 20 minutes while increasing either the spin or  0 increases the value of  * to 40 minutes -50 minutes.
In all simulations studied, the NIR total flux, unlike the the sub-mm      , is intrinsically noisy.This noise reduces the value of the peak cross-correlation coefficient between the NIR total flux and the sub-mm      .Performing a moving average of only the NIR total flux for the high-spin, prograde simulation with a window size between 120 GM/c 3 -240 GM/c 3 (40 minutes -80 minutes) preserves the periods of significant flaring while reducing the NIR total flux noise.When we correlate this averaged NIR total flux curve with the sub-mm      , the peak cross-correlation coefficient between the two curves increases from 0.6 to 0.7 -0.8.

The Probability of a Sub-mm Size Increase
Although the cross-correlation coefficient illustrates the relationship between the NIR total flux and the sub-mm size increases, the curve alone does not assign a probability of a size increase given a NIR total flux flare.To derive those probabilities, we first define an observation window to be a time-ordered set of all observables within a fixed time interval.Each observation window is indexed by its start time.We define the random variable  as a Boolean-valued outcome of the following trial: a success is if there is at least one point within an observation window of the sub-mm      above some threshold value     , where     is normalized by the median sub-mm      , and a failure otherwise.We then define     as the total NIR flux normalized by its median value (     =     /(    )).The Boolean random variable  is then defined to be true if the total flux during an observation window ever exceeds the threshold value     and false otherwise.We treat both  and  random variables as independent and identically distributed.We then calculate the conditional probability (    >  |     >  ) for two fixed duration observation windows: 50 minutes and 100 minutes.
Figure 13 shows these conditional probabilities for the fiducial high-spin, prograde simulation.An increase in NIR flux during an observation window greatly improves chances of detecting a sub-mm size increase during the same window.Even restricting observations to periods of minor (x2 the median value) elevated NIR total flux significantly increases the chances of a positive sub-mm size increase detection.For NIR flares that reach 8x the median NIR total flux, the chance of detecting at least a moderate size increase > 90%.Without considering the NIR total flux, the chance of detecting a moderate size increase is ∼ 20%.

The Sub-mm and NIR Flux
Unlike the consistent behavior of the cross-correlation coefficients between the sub-mm size and the NIR total flux, the cross-correlation coefficient between the sub-mm total flux and NIR total flux varies significantly across the simulation parameter space (See  Luminosity spectra for the high-spin, prograde simulation at key times before, during, and after the example NIR flare.At  = 36,200 GM/c 3 , the accretion flow emission is in quiescence.The spectrum is flat at low frequencies and falls exponentially at high frequencies.During the peak value of the NIR total flux flare ( = 37290 GM/c 3 ), the spectrum develops an elevated high frequency tail but is only slightly elevated at sub-mm and mm frequencies;   remains near zero.However, during the peak sub-mm size,   = -0.84.A time steps later, at 37,510 GM/c 3 , the NIR flux is no longer elevated but the size is sill 1.6x the median value and   remains negative.The sub-mm size increase is coincident with a shift in the synchrotron spectrum to lower frequencies and a more negative spectral index near 230 GHz.
wider than their counterparts between the sub-mm size and the NIR total flux.Although we do find correlations between the sub-mm total flux and NIR total flux, the correlations between the sub-mm size and the NIR total flux is both stronger and more consistent than the former across our parameter space.

The NIR Centroids
The NIR centroid motion varies considerably between flaring states and quiescent states.During quiescent periods, the NIR centroid stays near a median value with several "jumps" around the black hole.However, during flaring periods, the NIR intensity map is dominated by a bright region which moves in a smooth, continuous motion.For some extreme NIR flares, this continuous motion forms a quasicircular pattern on the sky.
This qualitative behavior occurs in all simulations studied.However, the effective diameter of the centroid path during different NIR flares varies considerably both during the course of a single simulation and across the simulation parameter space.The example from the fiducial simulation explored in Section 4 shows the largest (46 as) effective diameter an NIR centroid makes on the sky out of all simulations studied here.Other quasi-circular centroid motions during NIR flares from the same simulation show smaller effective diameters of ∼25 as and periods of 80-120 minutes.With the exception of the  = 0.3,  0 = 16 • simulation, which only produces centroid motion with a diameter of 17 as, the largest centroid motion diameters from other simulations explored range from ∼ 20as − 25as.
In the example from our fiducial simulation, the sub-mm centroid .The auto-correlation functions for the NIR Total Flux (red) and the the sub-mm maximum size increase (light blue) over the time 10,000 GM/c 3 to 96,310 GM/c 3 for the fiducial high-spin, prograde simulation.Besides the expected peak at  = 0, there is no significant secondary peak for a  range of from 0 GM/c 3 to 50,000 GM/c 3 , just over half the simulation time.Therefore, there are no significant periods of recurrence for both the NIR Total Flux and the sub-mm size increases.9) for the total NIR flux     , total sub-mm flux   , and the sub-mm max size increase      plotted a log scale for the fiducial high-spin, prograde simulation.Each curve is normalized by its mean value.The structure function for     and   reach a plateau value at ∼5 hours.The structure function for the sub-mm      shows some variation over after reaching a peak value but still stays near the mean, normalized value of 0.18. .The Total NIR Flux and sub-mm      for the full (10,000 GM/c 3 to 96,310 GM/c 3 ) high spin, prograde simulation.Even at late times, there are many examples of significant size increases following a period of NIR flares.Some smaller amplitude flares are able to produce large size increases while large amplitude flares can produce less extreme size increases.As explored in the Discussion, the actual size increase observed is dependent on the geometry of the low density, high temperature bubble.The size increases generally have a single peak with similar rising and falling behavior that corresponds to the ejection and then re-accumulation of material near the horizon, respectively.However, the NIR flares demonstrate a variety of morphologies.Some flares exhibit a single, extreme peak followed by lesser peaks (e.g. times ∼57,000 GM/c 3 and ∼92,500 GM/c 3 ) while other flares show several medium amplitude peaks in a row (near 23,000 GM/c 3 and 89,000 GM/c 3 ).also traces a quasi-circular orbit on the sky while the NIR centroid makes a similar orbit.However, we find many examples of NIR flares across our parameter space that exhibit the quasi-circular NIR motion without similar sub-mm centroid motion.However, these same NIR flare examples do display significant sub-mm size increases.
One interesting 9 hour period of time during the prograde  = −0.5 simulation contains 10 consecutive NIR flares.During this period, the NIR centroid makes several loops around the sky rather than the usual one or two loops.
We can also use the cross-correlation coefficient to show that the quasi-circular motion NIR centroid seen during NIR flares does not continue during quiescent periods.We treat the  and  centroid components, Ī and Ī , as separate curves and find the cross-correlation coefficient between them, CC [ Ī , Ī ] for short ∼2,000 GM/c 3 (11 hour) time intervals that contain a significant NIR flare and for the entire duration of the simulation.If the centroid traced a perfect ellipse on the sky with frequency ,then the maximum of the cross correlation coefficient CC [ Ī , Ī ]   = 1 at  = √ 2/2.For the short intervals of time, we find examples that reach CC [ Ī , Ī ]   of 0.6 -0.8.However, the values of CC [ Ī , Ī ]   for the full duration of each simulation in our parameter space is much lower at 0.2 -0.4.The upper-end of this interval comes from the  = −0.5 simulation where there is a 9 hour period of NIR flaring that results in many orbits of the NIR centroid around the black hole.

DISCUSSION
In this work, we show that magnetic flux eruptions cause a low density, high temperature magnetized bubble to form in the inner accretion flow.This bubble orbits around the black hole and breaks into smaller regions before reaching large (> 30  ) radii.The low density of the bubble as well as temperature increases in the  < 1 regions of the accretion flow lead to significant (1.5x -2x) increases in the sub-mm and mm emission region.These size increases follow NIR flares by ∼ 1 hour for the high-spin, prograde simulation.The mm size increase trails behind the sub-mm size increase by a few minutes, likely due to the fact that the mm emission region extends to larger radii compared to the sub-mm emission region.The submm size increase coincides with a sudden drop in the sub-mm total flux and can cause the sub-mm spectral index to become more negative.This relationship between NIR flaring and sub-mm/mm size increases occurs over a parameter space of black hole spin and initial tilt angle: parameters currently unconstrained for Sgr A*.
We showed that magnetic flux eruptions produce quasi-circular centroid tracks for simulations over a variety of black hole spin values; the largest diameter track occurred during the high spin case.We identified several examples of NIR flares where the NIR centroid follows a quasi-circular path but the sub-mm centroid does not.The sub-mm and NIR centroid motions are weighed by nearly complimentary regions in their respective emission maps (See Figure 4).There are also many sources that can effect the position of the submm centroid: the location of the shadow caused by the low density "bubble," additional heating outside the "bubble," and the anisotropic  3).Every column is a selected interval of time from a different simulation; the spin  and initial tilt angle  0 of the simulation are written in bold in each top left corner.The light red region highlights the interval of time where the NIR total flux is in a flare state.The light blue, dashed line aligns with the time where the sub-mm      is at a maximum.The black, horizontal solid and dashed lines are the observed median intrinsic (de-scattered) mm (86 GHz) major and minor axis, respectively, from Issaoun et al. (2021).The grey, shaded region corresponds to uncertainty (95% confidence interval) of the mm observed size.The blue, horizontal line represents the observed sub-mm size and the light blue shaded region is the boundary of the 68 % confidence interval of the observed sub-mm size (Event Horizon Telescope Collaboration et al. 2022a).In all the examples shown, the sub-mm size increase peak occurs a few minutes before the mm size increase peak.All simulations produce the observed sub-mm size during quiescent periods, but both the mm size and eccentricity is lower than the observed value.The cross correlation coefficients between the NIR and the sub-mm total flux show both wider peaks and smaller peak values than their counterparts between the NIR total Flux and sub-mm      .There is also greater variability of the cross correlation peak values, the peak lag times, and the peak widths across our parameter space for the NIR and sub-mm total flux versus NIR total flux and sub-mm      .Therefore, a correlation between the NIR and sub-mm total flux exists but is weaker and less consistent than the correlation between the NIR total flux and sub-mm      .
geometry of the accretion flow itself.Therefore, NIR centroid motion does not guarantee similar sub-mm centroid motion.
We also considered a single, high spin ( = 0.9375) SANE simulation with  0 = 16 • .After running the simulation for 25,000 GM/c 3 , we found one example of a small NIR flare (< 1 mJy) and no examples of significant sub-mm size increases.Therefore, the extreme variability in sub-mm size may be a unique feature of magnetic flux eruption events that occur in MAD simulations.
In our work, the accretion rate is fixed over 5,000 GM/c 3 periods of time by matching the calculated, average sub-mm total flux to its observed value (∼3 Jy) for Sgr A* .In the real Sgr A* system, stellar winds from the ∼30 Wolf-Rayet stars that populate the galactic center provide a significant source of accretion material (Ressler et al. 2018(Ressler et al. , 2019(Ressler et al. , 2020)).

Electron Heating, Non-Thermal Electrons, and 𝜎-cut
Self-consistently evolving the electron temperatures requires a choice of electron heating prescription.
In their study of sub-mm emission from a magnetic flux eruption, Jia et al. (2023) applied the  ℎℎ −   electron heating model (Mościbrodzka et al. 2016), where  ≡   /  .They demonstrated that, for high   models, setting  = 1, 10 and ray-tracing images with the  > 1 regions excluded produced a steady decline in the submm total flux and did not produce significant sub-mm size increases.However, setting  = 100 caused an increase in the sub-mm size.The sub-mm total flux also increased in the latter case before fading.
For this study, we chose the electron heating prescription W18 (Werner et al. 2018) because it was the most successful model in matching mm to NIR observations of Sgr A* (Dexter et al. 2020a) and producing NIR centroid motions.Models of particle-in-cell turbulent electron heating and heating during reconnection for an elec-tron/ion plasma is still an active field of study.Future results in this field will provide new prescriptions and more accurate predictions of observables from our models.Here, we also apply the W18 electron heating fraction uniformly regardless of whether reconnection is occurring.Future studies could account for spatial variation in the heating mechanism, instead of solely the fluid magnetization.
In our model, the electrons in the accretion flow are thermal.Relativistic reconnection (Sironi & Spitkovsky 2014;French et al. 2023;Hakobyan et al. 2023) the magnetic Rayleigh-Taylor instability (Zhdankin et al. 2023) are possible sources for non-thermal particles in the accretion flow.Non-thermal particles can significantly increase the NIR total flux while having a less pronounced effect on the submm/mm observables (Dodds-Eden 2009;Scepi et al. 2022).Nonthermal particles may also be the cause of the large mm size on the sky compared to what is calculated in this work.Non-thermal emission may also decrease the true NIR flux if electrons actually emit at higher frequencies (such as the X-ray).
In the images presented in our Results, we also excluded the high  > 1 regions of the accretion flow (-cut = 1).As shown in Figure 2, large regions of the innermost accretion flow ( ≲ 10  ) reach  > 1 values during a magnetic flux eruption.Setting -cut = 10 includes most of the highly magnetized bubble in the image while still excluding the jet.
Increasing -cut to 10 only produces a slight (2 as) increase in the sub-mm size during quiescence compared to its -cut = 1 counterpart.The -cut = 10 images from our fiducial example still show a sub-mm size increase reaching 2x its median value.The gap that forms in sub-mm intensity maps also persist in the high  image, demonstrating this gap is due to the drop in density in the region rather than an artificial omission of emission stemming from our choice of a low  cut.The difference in major axis of the peak sub-mm size between -cut = 10 and -cut = 1 is negligible (< 0.002 as) while the minor axis only increases in the higher -cut case by ∼ 3 as.The -cut = 10 NIR images produce roughly double the total NIR flux compared their -cut = 1 counterparts (see Dexter et al. (2020b) for further discussion on the choice of -cut and its effect on the NIR image.)

Electron Cooling
Electrons radiating at a frequency  lose energy and cool over a timescale   ∝  −1/2 .We do not include synchrotron cooling in our simulations.Recent results suggest that radiative cooling can effect the total luminosity and flux peaks in both the millimeter and NIR wavelengths (Yoon et al. 2020).Non-thermal electrons in the accretion flow will also radiate at high energies and cool.Particle-in-cell studies of electron/positron pair reconnection show that radiative models produce different power laws than their nonradiative counterparts (Werner et al. 2019;Sironi & Beloborodov 2020).
For the NIR photons,  ,    ∼10 minutes.Because  ,    is less than both the flare timescale and the duration of the quasi-circular NIR centroid orbits, cooling may have a significant impact on our results.
The sub-mm synchrotron cooling time, by contrast, is ∼4 hours and much longer than the time scale for the sub-mm size increases (∼ 2 hours) and sub-mm centroid motion.Also, the sub-mm total flux value is manually fixed.However, the sub-mm size increases are highly dependent on the hot, low density bubble from forming in the first place.If the hot NIR electrons cool too quickly in the real accretion flow, the bubble may dissipate more rapidly and not lead to such large sub-mm size increases.

MAD GRMHD Models and Simulation Resolution
Our models are subject to the same numerical computation constraints found in other magnetically dominated GRMHD simulations.Strongly magnetic regions are difficult to evolve but are present in the low density "bubble" that drives size increases and centroid motion.
Extreme resolution (5376×2304×2304) GRMHD simulations of black hole accretion flows produce magnetic flux eruption events that show reconnection features such as plasmoids and X-points (Ripperda et al. 2022).However, the extreme resolution simulations produced lower reconnection rates than its lower resolution counterpart because of the larger numerical resistivity in the latter case (Ripperda et al. 2022).Still, magnetic reconnection in GRMHD simulations is generally slower than in particle-in-cell (kinetic) simulations because of pressure anisotropies that are not well captured in MHD (Bessho & Bhattacharjee 2005;Bransgrove et al. 2021).Therefore, the duration of the magnetic flux eruption events in this work may be overestimated or underestimated.
Numerical diffusivity likely effects the boundary of the highly magnetized bubble as it orbits around and away from the black hole.The temperature increases directly outside the  > 1 region leads to changes in the sub-mm size and NIR total flux.Studying changes in sub-mm and NIR images from higher resolution GRMHD simulations may demonstrate if the sub-mm size increases (and their relationship to the NIR total flux shown in this work) depend strongly on our choice of resolution.At lower resolution (1/2 to 3/4 the number of cells compared to the resolution explored here), Dexter et al. (2020b) found that NIR hot spot formation and rotation around the black hole still occurs.
Despite the success of MAD accretion flows in explaining a variety of Sagittarius A* observations, there still exists some tension between simulated observables and real observations.Event Horizon Telescope Collaboration et al. (2022c) concluded that, although MAD models with inclination angles < 30 • were preferred, no model in their library of static images met all their observational constraints.

Dynamic Imaging with the EHT
EHT imaging algorithms are sensitive to ring morphologies (Event Horizon Telescope Collaboration et al. 2022b).In order to recreate the resolution of the constructed EHT image from interferometric data, we blur the 150 as x 150 as region centered on the black hole in our sub-mm simulated images with a 20 as FWHM Gaussian kernel.We focus on the blurred sub-mm before and during the fiducial NIR flare (Figure 14).Not only do we find that the resulting structure is ring-like in all images, but we can clearly see the sub-mm hotspot rotating around the black hole center.The effective shadow in the middle of the ring is also larger during a sub-mm size increase than during quiescence.Further analysis with dynamic imaging algorithms is required to confirm whether this variability can be reconstructed with current EHT or ngEHT ,  coverage.°.However, the exact value of the inclination angle to Sgr A* is still unclear.We varied the inclination angle for the time interval containing NIR flares from several simulations and studied the effect on the peak sub-mm size increase.The peak sub-mm size increase can remain prominent for very low (10 • ) and very high (73 • ) inclination angles.However, the exact relationship between the sub-mm peak value and the inclination angle can vary dramatically for different NIR flare examples.Figure 15 demonstrates two extremes of this relationship.In one example, the sub-mm peak size is significantly elevated from the median size value but remains flat as a function of inclination angle.In the second example, the sub-mm peak size grows roughly proportionally to the inclination angle.
Since 1.3 mm (230 GHz) is near the synchrotron spectrum peak, the intensity map at 1.3 mm is sensitive to the underlying density distribution of the gas.The actual shape of the the low-density, high temperature region resulting from the magnetic flux eruption event varies considerably.Therefore, the shape of the 1.3 mm emission region will also vary.Although the examples come from two simulations with different spin , it is unclear if  has an strong effect on this relationship between inclination angle and peak sub-mm size.In the pink example, the peak value of the sub-mm maximum size curve grows linearly with angle.By contrast, in the green example, the size increase stays relatively flat for a variety of inclinations.In both examples, the peak size increase is significantly elevated from its median value for all inclination angles.
size can be slightly (∼15%) larger than the observed sub-mm EHT size.By contrast, there is tension between our modeled mm size and the observed mm size.Our simulations consistently produce a maximum quiescent mm size that is about 30% lower than the observed value from the most recent analysis of ALMA data from April 3rd, 2017 (Issaoun et al. 2021).The quiescent eccentricity of the mm (86 GHz) emission is near zero in these simulations whereas the observed value is ∼0.7.Interestingly, the predicted size and eccentricity matches the ALMA data well during a mm size increase following a NIR flare (See Figure 3 from Section 4).Future observations may shed light on the true distribution of the mm size and its variability following NIR flares.
The fiducial NIR flare studied here still produces a NIR centroid track diameter that is a factor of 3 smaller than observed diameter from the GRAVITY 2018 measurement (GRAVITY Collaboration et al. 2018).With an effective diameter of 46as, the NIR centroid motion from the featured example is also the largest in our parameter space.NIR centroid motions with diameters of ∼ 15as − 25as are more common in our simulations but 7x -8x smaller than the observed value.Sources of non-thermal NIR emission missing from our model, such as particle acceleration and synchotron cooling, can account for the discrepancy in the NIR orbit size (Ball et al. 2021).

CONCLUSIONS
GRMHD simulations of magnetically arrested disks undergo magnetic flux eruptions that power NIR flares and may explain the NIR centroid motion observed by GRAVITY.Additionally, the same class of simulations predict significant and frequent sub-mm and mm size increases that occur within a 2 hour window after a NIR total flux peak.The models presented here make a clear, qualitative prediction within reach of the Next Generation Event Horizon Telescope.

Figure 1 .
Figure 1.Dimensionless magnetic flux on the horizon   (purple) and the NIR total flux curve (red) as a function of time for the high-spin ( = 0.9375) prograde simulation.The magnetic flux goes through many cycles of slow buildup followed by rapid dissipation.The dissipation of   , called a magnetic flux eruption event, coincides with a sharp peaks in the NIR total flux curve.

Figure 2 .
Figure2.The density  (in code units), the temperature in units of /   2 weighted by the density , and the plasma  parameter 2 / 2 , where  is the fluid pressure and  2 /2 is the magnetic pressure, all averaged over polar angle .We also show the plasma  parameter  2 /2 averaged over azimuthal angle  (see Equations 10 -13 for how we compute averages).Note that the outermost radius shown in the  parameter plots is less than the outermost radius shown for the other fluid parameter plots.A low-density, high temperature bubble forms in the innermost accretion flow and breaks apart before moving to larger radii during the 2,000 GM/c 3 interval outlined by a black box in Figure1.In the far left column, the accretion flow is in a quiescent state.The NIR total flux is near its median value and the sub-mm size is ∼50 as.The temperature and density is high close to the event horizon at ∼ 2  and then falls off at larger radii.Spiral structure in the density plot is visible.At  = 37,290 GM/c 3 , a large scale (10   ) low density, high temperature, and low plasma  region forms.Also, at this time, the magnetic flux is expelled from the horizon and NIR total flux curve is at a maximum.At  = 37,360 GM/c 3 the low density bubble has increased in size and pushed material out to larger radii.The temperature in this region stays elevated and additional smaller, even higher temperature spots appear throughout the bubble.In the rightmost column, at  = 37,450 GM/c 3 , the sub-mm emission region size is at a maximum.The bubble is dissipating.The temperature in the bubble has started to cool and the density of the bubble has increased.There is a small (10x) increase in temperature in the 15   -25   region.Following the trajectory of the bubble in the last three columns, one can see that it rotates counter-clockwise with the general accretion flow.

Figure 3 .
Figure3.The evolution of key NIR, sub-mm, and mm observables before and during a prominent NIR flare.The 2,000 GM/c 3 interval of time shown here corresponds to the black box outlined in Figure1and is taken from the high-spin, prograde simulation.In this example, the significant NIR flare results in both quasi-circular motion of the NIR centroid and an increase in the      and     of the sub-mm/mm emission region (See Section 2.3 for the definition of      ,       and the centroid position ( Ī , Ī )).The light red shaded regions corresponds to the second, prominent NIR flare; the dotted, vertical blue line represents the sub-mm      peak time.Following a period of quiescence, the NIR light curve first exhibits a small, ∼6 mJy flare near 36,900 GM/c 3 before a second, larger (22.2 mJy) flare that reaches its peak value at  = 37,290 GM/c 3 .The first flare corresponds to a minor eruption of magnetic flux whereas the larger flare occurs during a more dramatic magnetic flux eruption event.Prior to the second flare, both the NIR centroid (middle subplot) and the sub-mm/mm size (lower subplot) undergo only small deviations from their median values.However, during the second flare, the NIR centroid changes its behavior and makes a quasi-circular orbit around the black hole.The sub-mm and mm sizes begin to climb.The sub-mm peak size increase at 37,460 GM/c 3 coincides with a final, lesser NIR total flux peak during the second flare before the total flux quickly drops to its median value.The mm      reaches at maximum only 10 minutes after the sub-mm      peak and outside the flaring period.Then, the sub-mm and mm sizes slowly decrease back to their quiescent values.The entire event, from the beginning of the first flare to the end of the sub-mm/mm size increases, lasts 800 GM/c 3 , or just under 4.5 hours.However, the separation between the peak value of the prominent NIR flare and the sub-mm      peak is only 57 minutes (270 GM/c 3 ).The width of the sub-mm size increase peak is less than 2 hours.We also include the observed sub-mm size from the EHT (Event Horizon TelescopeCollaboration et al. 2022a)  and the observed mm sizes from ALMA(Issaoun et al. 2021) in the bottom panel of the plot.The light blue region around the blue EHT line represents the 68% credible interval; the light grey region around the observed mm      /    lines represents the 95% confidence interval for each measurement.We discuss the agreement between the sub-mm size and our model, as well as the tension between our calculated mm sizes and the observed values, in Section 8.6 of the Discussion.

Figure 4 .
Figure4.Select sub-mm (left) and NIR (right) images from the high-spin, prograde simulation that includes the example NIR flare.All times occur during the 2,000 GM/c 3 interval shown in Figure3.The NIR images are squares of 30  / 2 (150 as on the sky), the sub-mm images are squares of size 60  / 2 (300 as on the sky), and the black hole is centered at the origin in both sets of images.The minor and major axes of the red ellipse in the sub-mm correspond sub-mm       and      , respectively.The green + in the NIR images is the NIR centroid position at the time of the image and past locations of the centroid are represented with a green line.First, a gap in the sub-mm intensity map forms that corresponds to the region of increased NIR emission.The NIR centroid also begins its first, counter-clockwise orbit (from  = 37,200 GM/c 3 to  = 37,290 GM/c 3 ) around the black hole.At time 37,360 GM/c 3 , the gap in the sub-mm intensity map is extremely pronounced and the sub-mm size continues to increase.The NIR intensity in the region 20 as -50 as from the black hole center also increases, pushing the NIR centroid further away from its median value and starting a second, larger orbit around the black hole.By 37,460 GM/c 3 , when the sub-mm emission size is at a maximum, the gap in the sub-mm emission has dissipated but the emission in the 50 as -100 as has increased whereas the emission closer to the black hole center has decreased.As the sub-mm size decreases to its median value during times 37,500 GM/c 3 to 37550 GM/c 3 , the sub-mm intensity near the center of the black hole increases and the extra emission in the outer regions dissipates.The NIR total flux drops significant during this time, but the NIR centroid makes it back to its median value by 37550 GM/c 3 .By 37,900 GM/c 3 , both the sub-mm size and NIR centroid have returned to their median values.

Figure 5 .Figure 6 .
Figure 5. NIR and sub-mm centroid motion from the 2,000 GM/c 3 example which includes the prominent NIR flare, quiescent periods, and a minor flare.Both centroids stay close to their median values during quiescent periods and the first (minor) flare.However, during the second, more prominent flare, both the sub-mm and NIR centroids exhibit similar quasi-circular motion.The NIR centroid lags behind its sub-mm counterpart by ∼ 20 minutes.
Figure 12(right)).The resulting curves show a variety of peak values and values of the characteristic time  * .The peak of the cross correlation coefficients are all near or significantly less than 0.5 and are also

Figure 7 .
Figure 7. Luminosity spectra for the high-spin, prograde simulation at key times before, during, and after the example NIR flare.At  = 36,200 GM/c 3 , the accretion flow emission is in quiescence.The spectrum is flat at low frequencies and falls exponentially at high frequencies.During the peak value of the NIR total flux flare ( = 37290 GM/c 3 ), the spectrum develops an elevated high frequency tail but is only slightly elevated at sub-mm and mm frequencies;   remains near zero.However, during the peak sub-mm size,   = -0.84.A time steps later, at 37,510 GM/c 3 , the NIR flux is no longer elevated but the size is sill 1.6x the median value and   remains negative.The sub-mm size increase is coincident with a shift in the synchrotron spectrum to lower frequencies and a more negative spectral index near 230 GHz.
Figure8.The auto-correlation functions for the NIR Total Flux (red) and the the sub-mm maximum size increase (light blue) over the time 10,000 GM/c 3 to 96,310 GM/c 3 for the fiducial high-spin, prograde simulation.Besides the expected peak at  = 0, there is no significant secondary peak for a  range of from 0 GM/c 3 to 50,000 GM/c 3 , just over half the simulation time.Therefore, there are no significant periods of recurrence for both the NIR Total Flux and the sub-mm size increases.

Figure 9 .
Figure 9.The structure function (, Equation9) for the total NIR flux     , total sub-mm flux   , and the sub-mm max size increase      plotted a log scale for the fiducial high-spin, prograde simulation.Each curve is normalized by its mean value.The structure function for     and   reach a plateau value at ∼5 hours.The structure function for the sub-mm      shows some variation over after reaching a peak value but still stays near the mean, normalized value of 0.18.
Figure10.The Total NIR Flux and sub-mm      for the full (10,000 GM/c 3 to 96,310 GM/c 3 ) high spin, prograde simulation.Even at late times, there are many examples of significant size increases following a period of NIR flares.Some smaller amplitude flares are able to produce large size increases while large amplitude flares can produce less extreme size increases.As explored in the Discussion, the actual size increase observed is dependent on the geometry of the low density, high temperature bubble.The size increases generally have a single peak with similar rising and falling behavior that corresponds to the ejection and then re-accumulation of material near the horizon, respectively.However, the NIR flares demonstrate a variety of morphologies.Some flares exhibit a single, extreme peak followed by lesser peaks (e.g. times ∼57,000 GM/c 3 and ∼92,500 GM/c 3 ) while other flares show several medium amplitude peaks in a row (near 23,000 GM/c 3 and 89,000 GM/c 3 ).

Figure 11 .
Figure11.Examples of a NIR flare followed by a sub-mm and mm emission region size increase for every simulation studied in this work (with the exception of the fiducial simulation already shown in Figure3).Every column is a selected interval of time from a different simulation; the spin  and initial tilt angle  0 of the simulation are written in bold in each top left corner.The light red region highlights the interval of time where the NIR total flux is in a flare state.The light blue, dashed line aligns with the time where the sub-mm      is at a maximum.The black, horizontal solid and dashed lines are the observed median intrinsic (de-scattered) mm (86 GHz) major and minor axis, respectively, fromIssaoun et al. (2021).The grey, shaded region corresponds to uncertainty (95% confidence interval) of the mm observed size.The blue, horizontal line represents the observed sub-mm size and the light blue shaded region is the boundary of the 68 % confidence interval of the observed sub-mm size (Event Horizon TelescopeCollaboration et al. 2022a).In all the examples shown, the sub-mm size increase peak occurs a few minutes before the mm size increase peak.All simulations produce the observed sub-mm size during quiescent periods, but both the mm size and eccentricity is lower than the observed value.

Figure 12
Figure12.a) The cross correlation coefficient as a function of lag time  between the NIR total flux and the sub-mm      for all simulations studied in this paper.The dotted grey line represents a cross-correlation coefficient value of 0.6.Regardless of initial tilt angle  0 or spin , all simulations demonstrate a strong, singular cross-correlation coefficient peak of about 0.6 and a characteristic time delay  * between 20 -60 minutes.Simulations with low spin and low initial tilt have shorter peak lag times than those with high spins or moderate to high initial tilt angles.b) The cross correlation coefficients between the NIR and the sub-mm total flux show both wider peaks and smaller peak values than their counterparts between the NIR total Flux and sub-mm      .There is also greater variability of the cross correlation peak values, the peak lag times, and the peak widths across our parameter space for the NIR and sub-mm total flux versus NIR total flux and sub-mm      .Therefore, a correlation between the NIR and sub-mm total flux exists but is weaker and less consistent than the correlation between the NIR total flux and sub-mm      .

Figure 14 .
Figure 14.Simulated sub-mm (230 GHz) ground truth images using a 20 as FWHM Gaussian kernel for the fiducial NIR flare example shown in Figure 3 and 5.The field of view is 150 as; note that the color bar scale differs for the bottom row of images.The image intensity is plotted on a linear scale.The bright spot in the emission region clearly rotates around the black hole and the emission region also undergoes prominent morphological changes.

8. 6 Figure 15 .
Figure15.Two extremes of inclination angle dependence on the sub-mm (230 GHz) maximum size.Both curves show a size increase that occurred after a flaring event: the green line is taken from the  = 0.9375,  0 = 16 • simulation and the pink line comes from the  = 0.5,  0 = 16 • simulation.In the pink example, the peak value of the sub-mm maximum size curve grows linearly with angle.By contrast, in the green example, the size increase stays relatively flat for a variety of inclinations.In both examples, the peak size increase is significantly elevated from its median value for all inclination angles.

Table 1 .
A list of all simulations in this paper sorted by spin  and initial accretion flow tilt angle  0 .We ignore times prior to   to allow each simulation to evolve from its initial configuration and reach a MAD state.The time   is the duration of each simulation.We often refer to the high spin ( = 0.9375), prograde simulation as the fiducial simulation in this work.
00 − ( 10 / 00 ) 2  11 / 00 − ( 10 / 00 ) ( 01 / 00 )  11 / 00 − ( 10 / 00 )( 01 / 00 )  02 / 00 − ( 01 / 00 ) 2 The probability of detecting at least one point in the fiducial high-spin, prograde simulation where the normalized submm      is above some threshold   given that the NIR total flux is also elevated above the normalized flux     at some point during the same time interval (observation window).We define the sub-mm     =      /(     ) =      /61.5as,where(     ) is the median value of the sub-mm      curve.The quantity     =     /(    ) and the median value of the NIR total flux is (    ) = 2.2 mJy.The solid lines correspond to a long, 100 minute observation window and the dashed lines correspond to the short, 50 minute observation window.The blue lines correspond to  (     > ) and the other lines correspond to  (  >  |    / >  ) for an integer  .Even restricting observations to periods of minor (x2) elevated NIR total flux significantly increases the chances of a positive sub-mm size increase detection.The stronger the flare, the higher the chance of detection.It is possible for the probability of a submm size increase to reach 95% if it occurred within an observation window containing a NIR flare that increased by a factor of x8 from its median value (such flares have already been observed, seeGRAVITY Collaboration et al. Witzel (2018)19).(2019);Witzel(2018)).