Flares in the Galactic centre – II. Polarization signatures of flares at mm-wavelengths

ABSTRACT

Recent polarimetric mm-observations of the galactic centre showed sinusoidal loops in the |$\mathcal {Q{\small --}U}$| plane with a duration of one hour. The loops coincide with a quasi-simultaneous X-ray flare. A promising mechanism to explain the flaring events are magnetic flux eruptions in magnetically arrested accretion flows (MAD). In our previous work, we studied the accretion flow dynamics during flux eruptions. Here, we extend our previous study by investigating whether polarization loops can be a signature produced by magnetic flux eruptions. We find that loops in the |$\mathcal {Q{\small --}U}$| plane are robustly produced in MAD models as they lead to enhanced emissivity of compressed disc material due to orbiting flux bundles. A timing analysis of the synthetic polarized light curves demonstrate a polarized excess variability at time-scales of |$\simeq 1~\rm h$|⁠. The polarization loops are also clearly imprinted on the cross-correlation of the Stokes parameters which allows us to extract a typical periodicity of |$30~\rm min$| to |$1~\rm h$| with some evidence for a spin dependence. These results are intrinsic to the MAD state and should thus hold for a wide range of astrophysical objects. A subset of general relativistic magnetohydrodynamic simulations without saturated magnetic flux, namely, single temperature standard and normal evolution models, also produces |$\mathcal {Q{\small --}U}$| loops. However, in disagreement with the observations, loops in these simulations are quasi-continuous with a low polarization excess.

1 INTRODUCTION

The centre of our galaxy harbours one of the brightest radio sources in our sky, which is thought to coincide with the location of a putative supermassive black hole called Sagittarius A* (Sgr A*). The source has been extensively studied since its discovery in the 1970s (Balick & Brown 1974). More recently, the Event Horizon Telescope Collaboration reconstructed the first direct image of Sgr A* (Event Horizon Telescope Collaboration 2022a, c). The image shows a ring-like structure with a darkening in its centre, interpreted as gravitational lensing of the black hole’s event horizon – an effect often referred to as the black hole’s ‘shadow’ (Bardeen 1973; Luminet 1979; Falcke, Melia & Agol 2000). The size of the ring is |$50 ~\mu$|as, which agrees within 10 per cent with the prediction for a 4.1 × 10⁶ M_⊙ Kerr black hole whereby the mass is well constrained through the monitoring of stellar orbits (Gravity Collaboration ; Do et al. 2019a).

Sgr A* exhibits variability on all frequencies at time-scales from minutes to hours (e.g. Marrone et al. 2008; Witzel et al. 2018; Haggard et al. 2019; Iwata et al. 2020; Event Horizon Telescope Collaboration 2022b; Wielgus et al. 2022b). The strongest variable signatures are manifested as X-ray flares whereby the flux has been seen to increase several hundred-fold over a time-scale of ∼20 min (Haggard et al. 2019). The GRAVITY experiment observed Sgr A* during near infrared flares and found a bright component orbiting clockwise around the compact object at approximately ten gravitational radii (GRAVITY Collaboration 2018, 2023). Besides centroid motion, they also measured the polarization during the flare, which shows periodic variation in Stokes |$\mathcal {Q}$| and |$\mathcal {U}$|⁠, indicative of a hotspot that rotates through a poloidal magnetic field.

Similar |$\mathcal {Q{\small --}U}$| loops have been found at mm-wavebands (Marrone et al. 2006; Wielgus et al. 2022a). The most recent work by Wielgus et al. (2022a) analysed the Atacama Large Millimeter Array (ALMA) light curve obtained during the EHT campaign in 2017 April (see also Wielgus et al. 2022b). The data of April 11th commences just after an X-ray flare and indicates a drop of mm flux and a sinusoidal pattern in both Stokes |$\mathcal {Q}$| and |$\mathcal {U}$|⁠, leading to clockwise motion in a |$\mathcal {Q{\small --}U}$| diagram (hereafter: ‘the ALMA loop’). In Vos, Mościbrodzka & Wielgus (2022), the authors argue that the observation can be explained by a hotspot orbiting the black hole through a poloidal magnetic field. Although this simplified model successfully matches many aspects of the observation, it is uncertain if accretion flows will give rise to similar |$\mathcal {Q{\small --}U}$| loops if the flow geometry is more complex and dynamically varying. We must rely on general relativistic magnetohydrodynamic (GRMHD) simulations to answer this question.

Although Sgr A* is one of the brightest radio sources in the sky, its total bolometric luminosity is orders of magnitude lower than its Eddington luminosity such that Sgr A* is in the regime of Radiatively Inefficient and Advection Dominated Accretion Flows (RIAF/ADAF; Ichimaru 1977; Narayan & Yi 1994). Recently, the EHT measurements (Event Horizon Telescope Collaboration 2022a) allowed to put constraints on the state of the accretion flow which favour a magnetically arrested disc (MAD) model with an accretion rate |$(5.2{\small --}9.5)\times 10^{-9}\, {\rm M_\odot \, yr^{-1}}$|⁠. Similarly, the EHT measurements for M87* also favour the MAD state (Event Horizon Telescope Collaboration 2019, 2021).

In a MAD, the magnetic flux on the horizon has reached a maximum for a given accretion rate at which the accretion disc can be magnetically choked (Igumenshchev, Narayan & Abramowicz 2003; Tchekhovskoy, Narayan & McKinney 2011), in contrast to weakly magnetized accretion flows (SANE; e.g. Narayan et al. 2012; Porth et al. 2019). When a MAD model reaches its saturation state, the jet’s magnetic pressure and the disc’s ram pressure are in quasi-equilibrium. Accretion proceeds through spiral modes and the dynamics are highly intermittent: as more magnetic flux is accreted on to the black hole, eventually the jet will expel parts of its magnetic flux back into the disc in a bursty flux eruption. Whether the flux eruption is caused by an interchange instability between plasma on horizon penetrating field lines and the disc (e.g. McKinney, Tchekhovskoy & Blandford 2012) or via reconnection in an equatorial current sheet (Ripperda et al. 2022) or both has not been settled conclusively. In any case, the flux eruption generates an orbiting flux tube which is subsequently dissolved in the accretion flow, most likely by Rayleigh–Taylor-type instabilities (Ripperda et al. 2022; Zhdankin, Ripperda & Philippov 2023). After the flux eruption, accretion of matter and flux proceeds and the magnetic flux builds up until the next flux eruption event, leading to a repeating cycle of dissipation and accumulation of magnetic flux. The dynamics of orbiting flux tubes were studied in some detail in our first paper of this series (Porth et al. 2021) and by others (e.g. Begelman, Scepi & Dexter 2022; Chatterjee & Narayan 2022).

In contemporary MAD GRMHD simulations, the recurrence time of flux eruptions is ≈1000 GM/c³ corresponding to |$\simeq 5\, 1/2~\rm h$| for Sgr A*. Since near-infrared flares occur on a similar time-scale, i.e. 2–4 per day (e.g. Genzel et al. 2003; Dodds-Eden et al. 2011) – with roughly one in four also having an X-ray counterpart (e.g. Baganoff et al. 2001; Hornstein et al. 2007; Boyce et al. 2019) – flux eruptions seem a very good candidate for the physical process triggering such flares. Hence there is great interest in magnetic flux eruptions to explain AGN flares and several models have been applied to M87* as well as Sgr A* (Chatterjee et al. 2021; Dexter et al. 2020; Ripperda, Bacchini & Philippov 2020; Porth et al. 2021; Scepi, Dexter & C. 2022; Jia et al. 2023; Zhdankin et al. 2023). For the case of Sgr A*, the dissipation of magnetic energy can accelerate particles to non-thermal energies giving rise to efficient synchrotron emission up to X-rays. Additionally, the ejected flux tube orbits in the disc, which might account for the centroid motion observed by GRAVITY. While several studies have addressed the behaviour of total intensity at various wavelengths and various prescriptions for the underlying particle distribution, the polarization signatures of flux eruptions in MAD simulations are not yet well understood.

In this paper, we use GRMHD simulations in the SANE and MAD regime to characterize their linear polarization time-domain signatures with a particular focus on |$\mathcal {Q{\small --}U}$| loops. In Section 2, we will describe our GRMHD models and our polarized radiative transfer code. In Section 3, we will describe our results. In Section 4, we will discuss our findings and summarize our conclusions.

2 METHODS

2.1 General relativistic magnetohydrodynamic models

We analyse a number of ideal GRMHD simulations both in the SANE and in MAD regimes carried out with the BHAC code¹ (Porth et al. 2017; Olivares et al. 2019). We initialize our simulations from a Fishbone & Moncrief (1976) torus solution for different black hole spins a_* ∈ { − 0.9375, −0.5, 0, 0.5, 0.9375}, where a_* ≔ J/(Mc), with J being the black hole angular momentum, M the black hole mass, and c the speed of light. Negative spin indicates counter-alignment between the angular momenta of the accreting matter and black hole as customary. The simulations are part of the EHT simulation library and have previously yielded model comparison for the EHT 2017 Sgr A* campaign (Event Horizon Telescope Collaboration 2022d). The simulations are performed in spherical horizon penetrating Kerr–Schild coordinates r ∈ [1.19, 2500]r_g, where r_g ≔ GM/c², with G being the gravitational constant. We perform our analysis over the time-span |$t\in [15,30]\times 1000\, r_{\rm g}/c$| corresponding to roughly |$3.5\, \rm d$| when scaled to Sgr A*. Through three static mesh refinement levels, we obtain an effective grid resolution of N_r × N_θ × N_ϕ = 384 × 192 × 192 cells which is sufficient to resolve the magnetorotational instability throughout the evolution (Porth et al. 2019). For details of the employed numerical methods and setup, we refer to Porth et al. (2017), Olivares et al. (2019), Mizuno et al. (2021), and Event Horizon Telescope Collaboration (2022d).

2.2 General relativistic ray tracing

The GRMHD simulations provide a dynamical 3D model for the structure of the accretion flow. To compare our models to observations, we post-process the GRMHD simulations with our general relativistic ray-tracing code RAPTOR² to compute synthetic polarized images and light curves. RAPTOR computes null geodesics in curved space–times from a virtual camera outside our simulation domain. The camera is assigned pixels, and each is given an initial wave vector, after which the geodesic equation is solved to find their corresponding trajectories. Along these paths, we then solve the polarized radiative transport equation to obtain synthetic images for all four Stokes parameters, |$\mathcal {I}$|⁠, |$\mathcal {Q}$|⁠, |$\mathcal {U}$|⁠, and |$\mathcal {V}$|⁠. Observables are produced for a fiducial EHT frequency of |$229.1~\rm GHz$| and we consider thermal synchrotron emission (thus e.g. no bremsstrahlung), absorption, Faraday rotation, and conversion. The emission and absorption coefficients are obtained from the fit functions given by Dexter (2016) and rotativities from Shcherbakov (2008). In this study, we ignore effects caused by the finite propagation speed of the rays and assume the MHD state remains fixed during raytracing (known as the ‘fast light’ approximation). We also apply a ‘sigma-cut’ where we treat regions with magnetization σ ≔ B²/ρ > 1 as vacuum, where B is the magnetic field strength as measured in the fluid-frame, and ρ denotes the rest mass density. This is to exclude emission and absorption from the magnetized funnel region which is artificially mass-loaded by the GRMHD floor treatment (see also the discussions in Chael, Narayan & Johnson 2019; Porth et al. 2019). Our space–time is described by a Kerr–Schild metric for a rotating black hole. We set our virtual camera at a distance of d_camera = 10⁴r_g and appropriately re-scale the flux density ∼1/r² to the distance of Sgr A* – here assumed to be 8.127 kpc (Reid et al. 2019). The field of view of the camera is set to 44r_g, and we sample the image with 400² pixels.

2.3 Thermodynamics

Since the GRMHD simulations only provide information on the dynamically dominant ions, we must introduce parametrizations for the electron properties. First, assuming a pure hydrogen gas, charge-neutrality of the plasma dictates that n_p = n_e, where n_p is the proton number density and n_e is the electron number density. Secondly, we compute the electron temperature as

R then sets the ratio between the electron to proton temperature.³ For this ratio, we follow the prescription of Mościbrodzka, Falcke & Shiokawa (2016) by setting

where β = P_gas/P_mag is the ratio of thermal pressure P_gas and magnetic pressure P_mag. Via this prescription, we can vary the uncertain emission location using the factor R_high, where low values of R_high lead to relatively hot disc electrons promoting the emission from the disc region. Conversely, large R_high values lead to emission from the more magnetized regions near the jet wall. The following values are adopted in our study: R_high ∈ {1, 10, 40, 160}. As the simulations are all scaled to the same average 230 GHz flux density by changing the density normalization of the GRMHD data, increasing R_high leads to larger plasma densities. This is because larger values of R_high correspond to colder disc electrons which lower mm wavelength emission, consequently a larger density normalization is required to arrive at a given 230 GHz flux density. It is important to keep density normalization in mind when interpreting the results as it has consequences e.g. on the Faraday depth and polarization degree of the solutions.

During a magnetic flux eruption, particle acceleration due to magnetic reconnection can alter the shape of the electron distribution function (eDF; Ripperda et al. 2022; Hakobyan, Ripperda & Philippov 2023; Zhdankin et al. 2023). Non-thermal particles gain importance in radio- and IR frequencies (Özel, Psaltis & Narayan 2000; Davelaar et al. 2018) and are usually invoked to explain the inverted IR spectrum during strong flares (Do et al. 2019b) as well as the highest X-ray fluxes (Haggard et al. 2019). However, how and where to assign non-thermal eDFs in our GRMHD domain is very uncertain; parametrizations of the eDF can be done in various ways, which would enlarge our parameter space substantially. In this work, as a first step, we will hence only discuss mm-emission which is caused by the thermal core of the eDF and postpone a full eDF study to future work.

Lastly, the GRMHD simulations which are run in dimensionless geometric units (G = M = c = 1) under a ‘test-fluid’ assumption (neglecting the contribution of the plasma to the gravitational field) need to be scaled to concrete physical units. Thus we introduce scaling factors for plasma mass, length, and time. Units of length are scaled with |${\mathcal {L}}=r_{\rm g}$|⁠. Units of time are scaled with |${\mathcal {T}} = r_{\rm g}/c$|⁠, and the mass unit is set by the scaling factor |${\mathcal {M}}$|⁠. The mass unit scales plasma density via |${\mathcal {M}}/{\mathcal {L}^3}$| and is related to the mass accretion rate via |$\dot{M} = \dot{M}_{\rm sim}{\mathcal {M}}/{\mathcal {T}}$|⁠, where |$\dot{M}_{\rm sim}$| is the mass accretion rate in geometric units. The black hole mass is tightly constrained by observations of stellar orbits in the galactic centre; see e.g. GRAVITY Collaboration (2019) and Do et al. (2019a). Consequently, we set the black hole mass to M = 4.14 × 10⁶ M_⊙ which fixes |$\mathcal {L}$| and |$\mathcal {T}$|⁠. Since the plasma density is poorly constrained, for each parameter combination (R_high, i, a_*), we use the density normalization |${\mathcal {M}}$| as a free parameter to match the average flux density of our models to the observed flux density of |$2.4~\rm Jy$| (Event Horizon Telescope Collaboration 2022a).

3 RESULTS

3.1 Polarization loops during flux eruptions

To study the effect of flux eruptions on the polarization signatures, we first analyse the counter-rotating a_* = −0.9375 MAD simulation since it features the strongest and most frequent flux eruptions in our sample. The radiation transfer parameters are set to i = 10° and R_high = 1 such that we look nearly down the black hole spin axis and the emission is disc-dominated. We monitor the evolution of the polarized flux in the |$\mathcal {Q{\small --}U}$| plane over a sliding time window of 180 GM/c² ≃ 1 h, similar to the period observed by Wielgus et al. (2022a). The cadence is set by the output frequency of the GRMHD simulations to 10 GM/c², which results in sufficient data points to see coherent structures in the polarization time domain. As an example, we show three loops associated with flux eruptions in the top panel of Fig. 1. All selected cases show periods of ≃1 h and move in a counter-clockwise direction. We recover a range of polarized flux from 0.2 to 0.3 Jy resulting in varying loop radii. While the ALMA loop was significantly off-centre at |$(-0.19,0.01)~\rm Jy$| (as indicated by x on the figure), all our loops are centred close to the origin of the |$\mathcal {Q{\small --}U}$| plane. Our t₃ loop, similar to the ALMA loop, recovers a secondary ‘tiny loop’ close to the origin; see Vos et al. (2022) on an explanation of how such features arise in hotspot models. Apart from an overall flip (which can be obtained by instead observing at i = 170°), and offset, this ‘Pretzel’-shaped loop qualitatively shows similar behaviour to the ALMA loop.

We find a correlation between the formation of the loops and the slope of the horizon penetrating magnetic flux |$\Phi _{BH}:=1/2\int _0^{2\pi }\int _0^\pi |B^r|\sqrt{\gamma }d\theta d\phi$| where γ is the determinant of the three-metric and B^r is the radial magnetic field component as seen from the Eulerian observer. As magnetic flux is being accumulated at the event horizon (a positive slope, |$\frac{d}{dt} \Phi _{BH}\gt 0$|⁠), prominent large polarization loops are rare. However, flux eruptions (a negative slope, |$\frac{d}{dt}\Phi _{BH}\lt 0$|⁠) are almost always accompanied by loops in the |$\mathcal {Q{\small --}U}$| plane; see also Appendix A where we show the full sample of loops. It is worth noting that the ALMA loop occurred on the rising slope of flux density (corresponding to the re-accretion phase) where loops are less pronounced in our simulations.

In the second panel of Fig. 1, we simultaneously monitor the evolution of 230 GHz flux with the horizon penetrating magnetic flux Φ_BH. Overall, the data shows 13–15 large flux eruptions characterized by a sudden drop in Φ_BH which results in the typical recurrence time of |$\approx 1000 \rm M$|⁠. Generally, the post-eruption radiation flux is lower than the flux in the build-up phase. This behaviour was also noted by Jia et al. (2023). However, we also find several instances where the 230 GHz flux density sharply peaks during flux eruption. In our sample, the ‘peaking’ flux eruptions occur especially in earlier times where the intensity and optical depth are also the highest.

The intensity variations can be attributed to the evolution of the magnetic flux eruption, as illustrated in Fig. 2: during the expulsion of magnetic flux (hence for |$\frac{d}{dt}\Phi _{BH}\lt 0$|⁠), a cavity of low density plasma appears in the image at the photon ring (top-left panel). As more flux escapes from the black hole, the cavity grows in size and starts interacting with the accretion flow (top row, second, and third panel). The initial emission spike (see panel 2 of Fig. 1) is attributable to plasma heating which can be due to reconnection in the MAD equatorial current sheet and due adiabatic compression of the ambient plasma. Subsequently, as Fig. 2 illustrates, emission is most pronounced at the back of the flux bundle. This can be explained by emitting disc material that is compressed against the slower rotating flux bundle. The rotation of this localized emission region also leads to a change in the EVPA compared to the time average values (cf. black versus red ticks in Fig. 2). As the flux bundle expands in the accretion flow, the region of enhanced emissivity adiabatically cools and becomes sub-dominant (bottom row). Disc material is pushed away, leaving behind a large cavity of underluminous (but hot) jet plasma (bottom right panel). In our simulation, this emission ‘hole’ is responsible for the lower post-flare radiation flux. The EVPA has now performed a full loop while the flux bundle has progressed by roughly half of a circle. As the hole is dissolved, the intensity recovers to pre-flare levels on a range of time-scales (⁠|$1{\small --}5 ~\rm h$|⁠) accompanied by re-accretion of the magnetic flux. It is interesting to note that a similar recovery of the mm-flux on a time-scale of |$\sim 2~\rm h$| was also observed after the X-ray flare on 2017 April 11 (Wielgus et al. 2022b).

We further monitor the linearly polarized flux |$\mathcal {P}$| and spectral index α where F_ν ∝ ν^α in panel 3 of Fig. 1. Both quantities, but α in particular, correlate with the magnetic flux. The evolution of spectral index and polarization degree indicates a changing optical depth, consistent with a transition from optically thick to optically thin emission during the flux eruption event. In the aftermath, as the inner accretion disc is filled back with emitting plasma, intensity and spectral index climb back to pre-flare values. Similar trends, that is a spectral index decrease by 0.5 ± 0.2 across the mm peak as well as an excess of linear polarization during the flare (increasing from 9 to 17  per cent) were also observed by Marrone et al. (2008) in their light curve of the 2006 simultaneous X-ray and mm wavelength flare. Since the cooling time for electrons emitting in the mm band is much longer than the mm-flare duration, this evolution must be driven by adiabatic expansion of the flaring region, just like in our simulations which neglect radiative cooling.

Regarding the ALMA loop with data taken just |$\simeq 1~\rm h$| after the X-ray peak, it is noteworthy that no sub-mm flare was recorded on April 11, rather the |$230~\rm GHz$| lightcurve starts at a low value of |$\simeq 1.8\rm Jy$| followed by a recovery to usual quiescent values |$\sim 2.5~\rm Jy$|⁠. At the same time, the spectral index increases from ≃ −0.2 back to its quiescent value which is consistent with zero. In the context of our simulations, one interpretation of this spectral timing behaviour is that the ALMA loop was observed during the re-accretion phase characterized by increasing flux and spectral index on a time-scale of hours (cf. the large flux eruption event in Fig. 1 before |$20\, 000~\rm M$|⁠). An alternative interpretation to the recovery of the mm-flux density during the ALMA loop is that intense heating during the X-ray flare (coincident with the flux eruption) has shifted the characteristic frequency of the emitting electrons out of the mm-band. The observed recovery is then governed by the radiative cooling of the X-ray and IR emitting electrons. If the cooling time is fast enough this could lead to an increase of mm-flux during the ‘loopy phase’ just after the flux eruption (see also corresponding discussion in Wielgus et al. 2022a). This scenario should be tested with future simulations that take into account radiative cooling.

The behaviour of the |$\mathcal {Q{\small --}U}$| parameters during flaring events can also be seen in the time series of these fluxes. Focusing on the loop formed at t₃, in Fig. 3, we show the time evolution of the linear polarization components |$\mathcal {Q}$| and |$\mathcal {U}$| in comparison with the magnetic flux. A telltale sinusoidal pattern emerges immediately after the local peak of the magnetic flux, which corresponds to the loopy feature in the |$\mathcal {Q{\small --}U}$| plane. Identification of such features in polarized light curves should be fairly straightforward.

3.2 Inclination dependence

To test whether the observation of the rotating EVPA can place constraints on the inclination, we investigate the dependence of the polarization loops on the viewing angle. To this end, in Fig. 4, we show the loop formed in t₃ for three different inclinations with respect to the black hole spin axis: 10°, 30°, and 50°. Due to the loss of symmetry, the time-averaged polarization degree typically increases with larger viewing angles, however, this does not uniformly lead to an additional offset of the polarization loops from the origin. For example, for the t₃ event, only the 50° case is markedly offset and not contain the origin in the |$\mathcal {Q}{\small --}\mathcal {U}$| plane. In fact, the loop-averaged polarization in the t₂ event is even slightly closer to the origin when inclination is varied from 10° to 50°. The lack of systematic offset for low inclinations can be interpreted in the following way: during the flux eruption, the polarized emission from the loop boundary (see also Fig. 2) is dominant over any background ‘shadow’ contribution which is thus unable to significantly offset the average polarization. The size of the loops shrinks as we incline the line of sight; for inclinations larger than 50°, the linear polarization components show more stochastic behaviour, and loops are less clear. The loss of coherent loops for larger viewing angles can also be quantified by the corresponding decrease of correlation between |$\mathcal {Q}$| and |$\mathcal {U}$| discussed further in Section 3.5.

3.3 R_high dependence

We now turn to the dependence of the loops on the R_high parameter. While the fiducial simulation with R_high = 1 results in disc-dominated emission due to the large plasma temperature in the disc, increasing R_high progressively moves the emission site towards the jet-sheath. Furthermore, since the lower disc temperature is compensated by an increase in the number density of cold electrons (i.e. a larger |$\mathcal {M}$| is needed to arrive at the 2.4 Jy average flux density), there is also an increase of the Faraday depth with increasing R_high. The effect this has on the t₁ − t₃ loops is shown in Fig. 5. As can be gathered from the Figure, the changing emission site has a large effect on the appearance in the polarization plane. While for the t₁ event, the R_high = 1, 10 cases still have similar appearance and period, the t₂ and t₃ events show a more complex behaviour with the R_high = 10 case being smaller with added stochastic component. Upon increasing R_high further, the stochastic component takes over and the association between flux eruptions and prominent loops in the Q–U plane breaks down.

3.4 Polarization loops in non-flaring SANE simulations

In the previous section, we show that the loops can form during magnetic flux eruptions, in particular when the emission is disc-dominated; this implies that the loops in the |$\mathcal {Q}$| and |$\mathcal {U}$| plane could provide evidence of Sgr A* being in a MAD state. To test this hypothesis, we analyse our SANE simulations for the same spin, inclination and R_high. SANE models do not exhibit magnetic flux eruptions, as the magnetic flux settles down at a much lower value of around Φ_BH ∼ 0.25. In contrast to our hypothesis, the SANE R_high = 1 model does form loops. This behaviour can be seen in Fig. 6; on the top-left panel, several consecutive loops with a period of ∼1 h are visible. Another aspect of this behaviour can be traced in the bottom of Fig. 6 where the time series of the |$\mathcal {Q,U}$| (blue and red lines) is overplotted with the magnetic flux (green line). Multiple sinusoidal patterns are visible in |$\mathcal {Q,U}$|⁠; however, the magnetic flux does not exhibit substantial dissipation events, in contrast to the MAD models. This indicates that the loops must have a different physical origin.

To elucidate the origin of these features, in the top right panel of Fig. 6 we show the instantaneous map of the polarized flux |${\mathcal {P}}$| at 25 400 M, overlayed with the instantaneous electric vector (white) and time-average electric vector (red). In this snapshot, there is only a very small deviation between the average and instantaneous angle of the electric vector. However, the lengths of the white ticks change since the linear polarization is not constant as a function of time. The polarized intensity shows an enhancement of emission due to azimuthal m = 1 and m = 2 modes in the image plane (such that the intensity is enhanced once or twice on a circle around the origin), resulting in an asymmetry in the emission morphology. As these modes orbit around the black hole they change the orientation of the locally emitted EVPA in a periodic fashion, similar to the orbiting flux tube in the MAD case. However, in contrast to the flux tubes in the MAD case, the spirals are longer-lived, therefore generating multiple consecutive loops in the |$\mathcal {Q{\small --}U}$| plane. Additionally, due to the lower contrast of the emitting features (and since the overall EVPA distribution is nearly rotationally symmetric), the resulting net polarization is smaller than in the MAD case.

When moving the emission away from the disc when R_high ≥ 10, the net polarization decreases further. Additionally, the trajectories in the |$\mathcal {Q{\small --}U}$| plane becomes significantly slower and stochastic: in this case, we do not obtain any polarization loops consistent with the observations, see Appendix  B for more details. This behaviour is not seen in MAD simulations which show loops for any of the tested R_high parameters since the emission in MADs is generally less sensitive to R_high.

3.5 Timing properties of Q–U loops

To further analyse the systematic signal in the |$\mathcal {Q}$| and |$\mathcal {U}$| time series, we compute power spectrum densities of |$\mathcal {I}, \mathcal {Q}$|⁠, and |$\mathcal {U}$|⁠. These are shown in Fig. 7 for both MAD and SANE discs with different spins at i = 10°. For the case of a hot disc (R_high = 1), at around |$f = 1 \ \rm [1/h]$| a broad peak emerges for |$\mathcal {Q}$| and |$\mathcal {U}$| PSDs while the power spectrum for |$\mathcal {I}$| follows the typical red-noise accretion variability spectrum (e.g. Hogg & Reynolds 2016; Porth et al. 2019; Bollimpalli et al. 2020).

This qualitative behaviour is generally present for both MAD and SANE discs, however, in the cold-disc SANE case with R_high = 160 the |$\mathcal {Q},\, \mathcal {U}$| PSDs are broader and peak at much lower frequencies |$\sim 0.3/\rm h$|⁠. This can be explained by the emission site: for SANE, R_high = 160 the emission is located at larger distances in the jet. Since the relativistic jet rotation is expected to decrease with distance ∝ 1/r (Komissarov et al. 2009; Komissarov, Porth & Lyutikov 2015), this corresponds to larger periods compared to the disc mid-plane emission site. Accordingly, for the zero spin case with R_high = 160, we do not observe any significant periodicity due to the absence of jet rotation. With increasing R_high, coherent polarization patterns are additionally scrambled due to the increasing Faraday depth which adds a stochastic component thus hiding periodic signals.

By contrast, since R_high has less impact on the emission site in MAD discs which are highly magnetized also in the mid-plane (c.f. fig. 4 of Event Horizon Telescope Collaboration 2019), the PSDs for MAD discs are insensitive to the values of R_high explored in our work.

To further quantify the structural variability, we compute the normalized cross-correlation between |$\mathcal {Q}$| and |$\mathcal {U}$| for models with various spins. We use the Scipy Signal package to correlate the quantities

where |$\overline{\mathcal {Q}}$| and |$\overline{\mathcal {U}}$| are the time averages, N is the number of samples and |$\sigma _{\mathcal {Q}}$|⁠, |$\sigma _{\mathcal {U}}$| are the standard deviations of the fluxes. The considered time interval is the entire duration of the light curves |$t\in [15\, 000,30\, 000]\rm M$|⁠. Figs 8(a) and (b) depict the cross-correlations of MAD and SANE discs. For both disc configurations, the cross-correlation is antisymmetric about the origin, and the maximum and minimum occur at similar time lags for each spin – regardless of the disc being MAD or SANE. The sense of rotation is indicated by the slope of the cross-correlation through the origin, all our simulations feature counter-clockwise rotation and are thus sensitive to the rotation of the disc, not the black hole. For the co-rotating case, a_* = +0.9375 the maximum cross-correlation is at |$\tau = -20~\rm M$| and for the negative spin a_* = −0.9375 it is at |$\tau = -40 ~\rm M$|⁠. This implies that the linear polarization flux components vary with a time delay of 7 and 13 min, respectively. Interpretation of this result is straightforward if we model Q ∼ sin (2πt/P) and U ∼ cos (2πt/P) (cf. Fig. 3): at zero time-shift τ = 0, we obtain a loop with minimal cross-correlation; once U is shifted by a quarter period ±P/4, both signals vary in-phase/antiphase. While the peaks are fairly broad, this implies a typical periodicity of 27 to 54 minutes. The universal behaviour of the cross-correlations shows that |$\mathcal {Q}$|-|$\mathcal {U}$|loops are a robust prediction for low inclination MAD and R_high = 1 SANE accretion torus simulations. Upon increasing the inclination, the correlation of |$\mathcal {Q}$| and |$\mathcal {U}$| decreases: for example, while the maximum normalized cross-correlation for the MAD a* = 0.5 run with i = 10° is ≃ 0.6, it decreases to ≃ 0.5 for i = 30° and to ≃ 0.2 for i = 50°. Our analysis shows that the SANE R_high = 1 case shows a similar dependence on inclination. This is consistent with the loss of polarization loops for increasing inclination discussed in Section 3.2.

4 DISCUSSION AND CONCLUSIONS

One of our most striking results is that for viewing angles ≤50°, rotating EVPAs manifested as loops in the |$\mathcal {Q}-\mathcal {U}$| plane are copiously produced by GRMHD simulations in the MAD regime, in particular when the emission is disc-dominated. Loops are largely associated with violent flux eruptions leading to enhanced emissivity of disc material that is compressed against the orbiting flux bundles.

The loops in linear polarization also have a clear signature in the cross-correlation of |$\mathcal {Q}$| and |$\mathcal {U}$| as the correlation is maximized by phase shifting one of the observables by a quarter of a period. This allows us to deduce a typical loop period between 27 and 53 min. There is a tendency for smaller phase shifts (i.e. faster periodicity) in co-rotating cases compared to counter-rotating cases which imply that the loop periods are sensitive to black hole spin. While the cross-correlation peaks in our data are fairly broad which makes it difficult to determine exact phase shifts, the clear trends we obtain suggest that polarization time-domain data can provide an interesting avenue to constrain orbital motion and ultimately black hole spin.

Another striking feature shows up in the light curves of different Stokes parameters: while total intensity displays typical red-noise accretion variability, we observe a broad but significant variability excess in |$\mathcal {Q}$|and |$\mathcal {U}$| at periods of |$\approx 1\rm \, h$|⁠. This polarized excess variability can be explained as follows: while all Stokes vectors vary according to the fluctuations imprinted by the mass accretion rate (corresponding to radial motions), |$\mathcal {Q}$| and |$\mathcal {U}$| are additionally sensitive to the azimuthal motions of emitting features such that they sample the orbital dynamics of the innermost accretion flow. It will be very interesting to search for similar time-domain signatures in observed linear polarization light curves.

In our sample, persistent polarization loops are also visible for SANE cases with disc-dominated emission (R_high = 1). In this case they are produced by m = 1 and m = 2 modes in the polarized intensity most likely corresponding to spiral overdensities in the accretion flow. Since the non-axisymmetric modes have relatively low contrast in polarized intensity and the EVPA distribution is close to rotationally symmetric, a low net polarization results in the SANE cases. Together with the long lifetime of the spiral features, we obtain near continuous loops with small radii in the |$\mathcal {Q}{\small --}\mathcal {U}$| plane.

That being said, it is important to stress that both from theoretical grounds which demand T_e < T_ions and from observational modelling of Sgr A* (Event Horizon Telescope Collaboration 2022d), these models (SANE, R_high = 1) are disfavoured. Already at R_high = 10, we do not recover any polarization loops with hourly periods any more. This suggests that the observed loops by Wielgus et al. (2022a) are generated by an accretion flow that is in the MAD regime.

For the quiescent state of SgrA*, the observations of the EHT typically favour models where the electron temperature in the disc is low compared to the proton temperature, so R_high > 1. In our analysis, we found that |$\mathcal {Q}{\small --}\mathcal {U}$| loops are suppressed for large R_high values, which at first glance suggests a tension between our models and the observations. However, in our work, we assumed the electron distribution function to be a relativistic thermal Maxwellian. This is an oversimplification of the complex electron thermodynamics during a flare since we expect the electrons to be accelerated via magnetic dissipation in the equatorial current sheets that form during a MAD flux eruption (Ripperda et al. 2022). Additionally, as the resulting flux tube orbits through the higher density accretion flow, Rayleigh–Taylor instabilities could provide additional heating of the electrons (Zhdankin et al. 2023). Both of these effects would push the electron energies around the flux tube to larger values compared to quiescence, a behaviour that is mimicked by the R_high = 1 models, which has a more disc-dominated emission morphology. The addition of power laws in the electron distribution function will be studied in future works.

MAD flux eruptions strongly perturb the magnetic field in the emitting region, resulting in local polarization vectors which differ significantly from an average ‘background’ magnetic field/EVPA configuration. Emitting spiral features in SANE simulations on the other hand simply enhance the local polarized emissivity, while keeping the direction of the EVPA largely intact. This means that inferring magnetic field topology from ‘hotspot’ models should be done with some care since (MAD-) flares are capable of significantly perturbing the magnetic field in the emitting region. For the more docile R_high = 1 SANE simulations; however, our results lend some support to treating the magnetic field as a fixed background as commonly done in hotspot models (e.g. Vos et al. 2022). This in principle allows us to relate the orbital frequency of emitting features with the period in the |$\mathcal {Q}{\small --}\mathcal {U}$| plane: in the SANE case which is dominated by the toroidal magnetic field component, as the background EVPA rotates twice along one orbit (cf. Fig. 6), the typical orbital period corresponds to twice the polarization period (see also the discussion in GRAVITY Collaboration 2018; Vos et al. 2022).

Regardless of the inclination, in our study, the polarization loops almost always enclose the |$\mathcal {P}=0$| origin, whereas the ALMA data reported in Wielgus et al. (2022a) features a clear offset with |$\overline{\mathcal {P}} =0.19 \rm ~ Jy$|⁠. In the latter study, the authors argue that the total polarization flux is a sum of the two components: |$\mathcal {P} = \mathcal {P}_{\rm hsp} + \mathcal {P}_{\rm sh}$|⁠, where hsp stands for the varying hotspot and sh stands for the background black hole shadow emission with |$\mathcal {P}_{\rm sh} \approx 0.2 \rm ~ Jy$|⁠. Since our model includes the entire emitting region within |$0.2\rm ~mas$|⁠, any background shadow emission is consistently taken into account. In addition, for an inclination <50°, the mean linear polarization degree in all our MAD models is relatively low, in particular below |$5{\small --}10~{{\ \rm per\ cent}}$| fractional linear polarization commonly observed from Sgr A* (see appendix H in Collaboration 2024) for a quantitative overview of the BHAC models. While the obtained polarization fractions are quite sensitive to the uncertain electron thermodynamics, the lack of an overall offset in our simulations could indicate that another polarized emission component outside our field of view might be required.

Despite the different underlying assumptions between the GRMHD simulations and the hotspot models of (Vos et al. 2022; Wielgus et al. 2022a) similar features can be produced by both models. Specifically, our t₃ loop shown in Fig. 1 is very similar to the pretzel-looking loop shown in fig. 6 of Vos et al. (2022). The period of both loops is about |$180\rm ~M\, \sim 60~\rm min$| and similar periods are reported by GRAVITY for a number of flares (GRAVITY Collaboration 2023).

The features observed so far (GRAVITY and ALMA loops) move in a clockwise direction whereas loops and emitting region in our simulations rotates counter-clockwise. We have verified that clockwise rotation is obtained when we choose an inclination >90°; this data is shown in Appendix  C. Apart from the directionality, our conclusions remain unaffected when flipping the viewing angle.

There are several recent studies that investigate GRMHD simulations in the MAD regime to model black hole flares (e.g. Chatterjee et al. 2021; Ripperda et al. 2022; Scepi et al. 2022; Jia et al. 2023). Although we focus on the polarization properties, it is interesting to also compare the behaviour of total intensity during the flares. In particular, Jia et al. (2023) found that except for very cold disc electrons (R_low, R_high) = (100, 100), the 230 GHz flux dims during eruptions. We instead find that sharp flux increases during the eruption are recovered, in particular when the flux is above |$\approx 2~\rm Jy$| (cf. Fig. 1). The difference is likely explained by the larger optical depth obtained in our simulations which have been scaled to Sgr A*, not M87* as in the comparison paper. Indeed, as Jia et al. (2023) discussed, the dimming can be understood from the temperature dependence of the emissivity alone. This explanation, however, only holds in the optically thin regime. In Sgr A* models, the optical depth is close to one in quiescence and can reach τ = 10 at the eruption site. Since we model thermal emission and the Planck function increases monotonously with temperature for any frequency, flux increase is instead expected during optically thick flaring events.

On a related note, the |$230~\rm GHz$| spectral indices obtained in our fiducial model (α ∈ [−1, −0.2]) are too small to match the observations which find α ≥ 0 in quiescence and only slightly negative ≈−0.2 post flare (Goddi et al. 2021; Wielgus et al. 2022b). This implies that our models might still be too optically thin. In fact, as also shown by Ricarte et al. (2023), for the adopted electron temperature prescriptions following equation (2) with R_high ∈ {1, 10, 40, 160} the spectral indices are negative for essentially all the SANE and MAD models considered here. We suggest that models with colder electrons should be tried for Sgr A* which will increase the plasma density and optical depth of the models and thereby increase the spectral index.

Since EHT and GRAVITY should be able to resolve the centroid of the flaring emission on a scale of |$\sim 10 \rm ~M$|⁠, it will be interesting to assess the imprint of the flux eruptions on interferometric observations. In particular, it is important to predict the radius and motion of the emission centroid. We leave such an investigation to future work. Furthermore, including radiative cooling in the GRMHD simulations might be worthwhile since fast-cooling heated electrons can cool back into the mm-emitting range on a short time-scale, thereby changing also the mm-light curves.

ACKNOWLEDGEMENTS

We would like to acknowledge Hector Olivares and Christian Fromm for stimulating discussions during this work. We thank Maciek Wielgus for useful comments on a draft. JD is supported by NASA grant NNX17AL82G and a Joint Columbia/Flatiron Postdoctoral Fellowship. Research at the Flatiron Institute is supported by the Simons Foundation. YM is supported by the National Key R&D Program of China (grant no. 2023YFE0101200), the National Natural Science Foundation of China (grant no. 12273022), and the Shanghai municipality orientation program of basic research for international scientists (grant no. 22JC1410600). This research has made use of NASA’s Astrophysics Data System.

Software:python (Oliphant 2007; Millman & Aivazis 2011), scipy (Jones et al. 2001), numpy (van der Walt, Colbert & Varoquaux 2011), and matplotlib (Hunter 2007).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: SAMPLE OF LOOPS FOR THE ERUPTION AND RE-ACCRETION PHASES

The fiducial inclination for our analysis is i = 10°, which reveals loops that rotate counter-clockwise. In Fig. C1, we show that flipping the inclination to i = 170° still generates loops that rotate clockwise and are slightly different in shape compared to their counterparts in Fig. 1. The majority of the loops in the |$\mathcal {Q{\small --}U}$| plane form during the eruptive period of the magnetic flux. In order to be able to inspect this correlation in more detail, we show the |$\mathcal {Q{\small --}U}$| plane during a selection of the eruptive and accumulating periods of the magnetic flux. For an overview of the sample, see Fig. A1,

The maxima of the magnetic flux are marked in red triangles and the minima in blue. These time markers are the initial time of the inspection when we begin to monitor the |$\mathcal {Q{\small --}U}$| as the magnetic flux evolves for about 1 h. In Fig. A2, we see the corresponding |$\mathcal {Q{\small --}U}$| for the selection of eruptive periods; in almost all panels, we can recognize a loopy feature. On the other hand, for the re-accretion periods in Fig. A3, most panels depict a more stochastic pattern than the loop-shaped curves.

APPENDIX B: R_high DEPENDENCE OF SANE LOOPS

The response of the polarization signal upon changing the electron temperature parameter R_high is shown in Fig. B1 for a fiducial SANE simulation. We choose an interval of an hour starting at |$25260~\rm M$|⁠, thus within the time-span also shown in Fig. 6. Hourly polarization loops are only recovered with our most disc-dominated model R_high = 1. For larger values, motion in the |$\mathcal {Q}$|-|$\mathcal {U}$| plane is slower with a stronger stochastic component. There are two main underlying effects behind this: (i) the change of the emission region which moves to larger distances hence slower time-scales, (ii) the increase of the Faraday depth when cooling of the electrons in the disc contributing to more stochastic polarization signals.

APPENDIX C: FLIPPED VIEWING ANGLE

Here, we show the selected flux eruption events from Fig. 1 when the viewing angle is flipped with respect to the equatorial plane. Fig. C1 thus adopts a viewing angle of 170°, showing that the loops now move clockwise in the |$\mathcal {Q{\small --}U}$| plane. Since the track is sensitive to small-scale asymmetries in the accretion flow, the loops are not exactly mirrored versions of Fig. 1; however, their overall extent and period remains similar.