New Constraints on Cosmic Particle Populations at the Galactic Center using X-ray Observations of the Molecular Cloud Sagittarius B2

The Sagittarius B2 (Sgr B2) molecular cloud complex is an X-ray reflection nebula whose total emissions have continued to decrease since 2001 as it reprocesses one or more past energetic outbursts from the supermassive black hole Sagittarius A*. The X-ray reflection model explains the observed time variability and provides a window into the luminous evolutionary history of our nearest supermassive black hole. In light of evidence of elevated cosmic particle populations in the Galactic Center, X-rays from Sgr B2 are also of interest as a probe of low-energy (sub-GeV) cosmic rays, which may be responsible for an increasing relative fraction of the nonthermal X-ray emission as the contribution from X-ray reflection decreases. Here, we present the most recent deep NuSTAR and XMM-Newton observations of Sgr B2, from 2018. These reveal small-scale variations within lower-density portions of the complex, including brightening features, yet still enable upper limits on X-rays from low-energy cosmic particle interactions in Sgr B2. We present Fe K$\alpha$ fluxes from cloud regions of different densities, facilitating comparison with models of ambient LECR interactions throughout the cloud.


INTRODUCTION
Centered at ∼100 pc projected distance from the supermassive black hole Sagittarius A* (Sgr A*; Ghez et al. (2008)) at the dynamic center of the Galaxy, and ∼8 kpc from Earth, the molecular cloud (MC) Sagittarius B2 (Sgr B2) is the densest and most massive such object in the Central Molecular Zone (CMZ), a region that extends several 100 pc from Sgr A* and contains ∼10% of the Galaxy's total molecular material (Morris & Serabyn 1996). X-ray observations of Sgr B2 have revealed strong Fe Kα line emission at 6.4 keV (Koyama et al. 1996;Murakami et al. 2001;Koyama et al. 2007;Inui et al. 2009;Terrier et al. 2010;Nobukawa et al. 2011;Terrier et al. 2018) as well as a hard continuum up to ∼100 keV (Sunyaev et al. 1993;Revnivtsev et al. 2004;Terrier et al. 2010;Zhang et al. 2015). These features, which imply energetic, nonthermal interactions capable of ionizing the K-shell electrons of neutral Fe, have made Sgr B2 an object of interest for decades. The X-ray picture is further complicated by the time-varying nature of this emission. Since the peak flux was last observed in 2001, the Fe Kα emission has decreased with every subsequent observation, down to ∼20% of the peak by 2013 (Zhang et al. 2015;Terrier et al. 2018). The hard continuum emission from the complex has correspondingly decreased, by ∼50% from 2003 to 2019 (Kuznetsova et al. 2021).
In the X-ray reflection nebula (XRN) model, the Fe Kα X-rays originate in the reprocessing of external Xrays via K-shell photoionization and subsequent fluorescence of neutral Fe gas while the continuum emission arises from Rayleigh and Compton scattering (Sunyaev et al. 1993;Koyama et al. 1996;Sunyaev & Churazov 1998). Reprocessing of X-rays originating in past flaring activity of Sgr A* is the widely accepted explanation of the time-variable nonthermal emission, where the timevariability emerges as the flares pass in and out of the MC (Sunyaev et al. 1993). Meanwhile, emission due to multiple scattering is expected from the densest regions even after the flare has exited the cloud (Sunyaev & Churazov 1998;Odaka et al. 2011;Molaro et al. 2016).
A short (<10-year) and bright event taking ∼10 − 20 years to traverse the cloud could explain the peak luminosity from the Sgr B2 core as well as the subsequent dimming (see Terrier et al. (2018) and references therein). Though direct observation shows that Sgr A* is presently in a quiescent state (Baganoff et al. 2003;Wang et al. 2013;Corrales et al. 2020), the XRN picture of Sgr B2 and other MCs in the CMZ reveals that Sgr A* has been brighter in the past few hundred years, with at least two short outbursts Clavel et al. 2013;Ponti et al. 2013;Churazov et al. 2017b;Chuard et al. 2018;Terrier et al. 2018).
A portion of the X-ray emission from Sgr B2 could arise from low-energy (<1 GeV, i.e. highly ionizing) cosmic-ray (LECR) electrons or protons, where the Fe Kα line arises from K-shell ionization of neutral Fe and the continuum arises from Bremsstrahlung processes (Valinia et al. 2000;Yusef-Zadeh et al. 2007;Dogiel et al. 2009b;Tatischeff et al. 2012). The observed rates of hydrogen ionization in the Galactic Center (GC) region, which are in excess of local rates by a factor of ∼10, require models with elevated GC LECR populations relative to the local galactic environment (Indriolo & Mc-Call 2012;Le Petit et al. 2016;Oka et al. 2019). Neither LECR electrons nor protons can explain the full timevarying flux. The cooling time for ∼100 MeV protons is longer than the observed timescale of the decrease , and the proton population corresponding to the hydrogen ionization cannot explain the Fe Kα emission in the bright state (Dogiel et al. 2013). Meanwhile, cooling of LECR electrons could explain the time variability (Yusef-Zadeh et al. 2013), but the peak X-ray flux cannot be easily explained by LECR electrons alone (Revnivtsev et al. 2004), requiring a highly-tuned model, e.g. higher metallicity in Sgr B2 than surrounding clouds (Yusef-Zadeh et al. 2013). However, any steady-state LECR population in the cloud contributes a constant nonthermal X-ray flux, in addition to the XRN flux.
Evidence of elevated cosmic-ray particle (CR) populations across a broad energy range in the vicinity of the CMZ (Aharonian et al. 2006;Yusef-Zadeh et al. 2013;Zhang et al. 2014;HESS Collaboration et al. 2016;Heywood et al. 2019;Zhang et al. 2020;Ponti et al. 2021;Huang et al. 2021) motivates a discussion of particle accelerators in the GC and the role of the GC in the dynamics in the Galaxy at large. Despite decades of discussion (Ptuskin & Khazan 1981), these are not presently emphasized in Galactic CR propagation models. The models (i.e. GALPROP (Strong et al. 2007), DRAGON (Evoli et al. 2017)) currently used to predict fluxes at Earth assume CR acceleration in supernova remnants (SNR) distributed throughout the Galactic disk and largely ignore the GC region. A portion of CRs at Earth could be due to extreme processes at the GC (Cheng et al. 2012;Jaupart et al. 2018;Anjos & Catalani 2020), although other recent work suggests that CRs originating in the GC are not able to escape (Huang et al. 2021). LECRs in the Galaxy are excellent tracers for sites of CR acceleration, due to their slow propagation and high rates of energy loss, and are excellent probes of cosmic-ray propagation models (Liu et al. 2021).
Measuring X-ray flux levels enables setting upper limits on ionizing power from LECRs within a given region of Sgr B2 (Valinia et al. 2000;Dogiel et al. 2009a,b;Terrier et al. 2010;Tatischeff et al. 2012;Dogiel et al. 2013;Zhang et al. 2015). The flux levels alone can only produce upper limits due to uncertainties in both the time-varying contribution from primary XRN flares and the contribution from multiple scattering, which can vary with longer timescales (see, e.g., the discussion by Chernyshov et al. (2018)). Limits on X-ray emission from LECR propagation in MCs can probe ambient LECR populations in the CMZ, providing valuable input for CR models, though the ability of LECRs to traverse dense MCs is highly model dependent. The interactions of LECRs within clouds is of particular interest due to the impact of the local ionization environment inside clouds on star formation (Morlino & Gabici 2015).
In this paper, we present recent deep observations of Sgr B2 obtained in 2018 by XMM-Newton and NuSTAR. Section 2 details the observations and data preparation. In Section 3 we show the X-ray morphology of Sgr B2 while in Section 4 we present spectral analysis of the central region. In Section 5 we compare the 2018 flux with earlier data to discuss the continued decrease in Xray reflection since 2001, while in Section 6 we present our main results on upper limits on Fe Kα emission from ambient LECR proton and electron populations in different regions of Sgr B2. Finally, in Section 7 we discuss these results in the GC context. Table 1 lists the observations discussed in this work. We present new observations of Sgr B2, taken jointly by the XMM-Newton and NuSTAR X-ray observatories in 2018. For comparison with the 2018 data, we also use archival XMM-Newton observations of Sgr B2.

XMM-Newton Observations
XMM-Newton consists of three European Photon Imaging Camera (EPIC) instruments: two Metal Oxide Semiconductor (MOS) arrays and a pn array. These cameras detect X-rays from 0.15 − 15 keV with typical energy resolution of ∼2 − 5% and angular resolution of 6 FWHM (Turner et al. 2001;Strüder et al. 2001).
Analysis was performed using the XMM-Newton Extended Source Analysis Software (ESAS; Snowden et al. (2008)) distributed with v.12.0.1 of the XMM-Newton Science Analysis Software. We reduced the data using the standard procedure and filtered the event files to exclude intervals affected by soft proton contamination.
Spectra were extracted with the ESAS mos-spectra and pn-spectra scripts. The MOS1 and MOS2 spectra were combined and all were rebinned with 3σ significance after background subtraction. We analyzed MOS and pn data within the 2 − 10 and 2 − 7.8 keV bands, respectively, where the pn spectra were truncated due to internal lines around 8 keV.

NuSTAR Observations
NuSTAR operates in the 3 − 79 keV band using two focal plane modules (FPMA and FPMB) with angular resolution of 18 (FWHM) and typical energy resolution of 400 eV (FWHM) at 10 keV (Harrison et al. 2013).
The data were reduced and analyzed using the NuS-TAR Data Analysis Software (NuSTARDAS) v.1.3.1 and HEASOFT v.6.24 (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014). They were filtered for periods of high instrumental background due to South Atlantic Anomaly passages and according to a database of bad detector pixels. The data quality was impacted by stray light (unfocused photons arriving directly onto the detector at large off-axis angles) from both bright isolated sources and diffuse X-ray backgrounds. We removed the pixels contaminated by stray light from bright isolated sources using a geometrical model of the telescopes following Krivonos et al. (2014). The FMPB observation was disregarded because the removed pixels covered the Sgr B2 region. In FPMA, pixels were removed as close as ∼50 from the center of Sgr B2.
In contrast to the bright isolated sources, stray light from diffuse backgrounds fills the entire detector area with a non-uniform pattern, the brightness of which limits the signal-to-noise. We therefore use NuSTAR spectra in the range of 10 − 20 keV, where the signal-to-noise ratio is highest.
3. MORPHOLOGY OF X-RAY EMISSION Figure 1 presents the 2018 observations in the 24 ×24 region centered on Sgr B2. The upper panel shows the XMM-Newton pn images in the continuum-subtracted 6.4 keV line and in the 2 − 5 keV and 5 − 10 keV bands. The lower panel shows the NuSTAR FPMA images in the continuum-subtracted 6.4 keV line and in the 10 − 20 keV and 20 − 79 keV bands. The 6.4 keV line images were created by subtracting a continuum band, 5.8 − 6.2 keV, from a 6.2 − 6.6 keV signal band. The 90 (XMM-Newton) and 50 (NuSTAR) source regions used for the primary spectral analysis of the core and envelope are shown in dark blue, where the smaller NuSTAR spectral extraction region is due to stray light contamination (see Section 2.2). The background regions used for spectral analysis with each instrument are in green. The XMM-Newton background region is located outside the Sgr B2 complex. In contrast, the NuSTAR background region is located within the diffuse region of Sgr B2, due to the limited field of view (see Section 2.2). All regions shown in Figure 1 are listed in Table 5 in the Appendix.
The XMM-Newton 6.4 keV line map shows that the core of Sgr B2 is detected at 13σ significance within the envelope region. The core is also detected at >5σ The 2018 X-ray morphology of the 24 × 24 region surrounding Sgr B2 is shown as observed by XMM-Newton pn (top) in the 6.4 keV line (left), 2 − 5 keV (center), and 5 − 10 keV (right) bands; and by NuSTAR FPMA (bottom) in the 6.4 keV line (left), 10 − 20 keV (center), and 20 − 79 keV (right) bands. The 6.4 keV line images are continuum subtracted as in Section 3. Contours (white) of the XMM-Newton 6.4 keV map are overlaid on all images and illuminate the core and envelope of Sgr B2 as well as several substructures, labelled by their Galactic coordinates. The annular stray light pattern observed in all EPIC instruments is most evident in the 2 − 5 keV band (top center, black). The stray light in FPMA from the bright off-axis source is evident in the curved region removed from the top-left of the images, while the diffuse stray light is apparent in the anomalously bright regions of the 10 − 20 keV and 20 − 79 keV bands. Circles (white) indicating the diffuse (9.9 radius) and envelope (2.2 radius) regions of the simplified model are overlaid, while the core (2 − 4 radius) is smaller than the angular resolution of both telescopes. The locations of the brightest (>10 −6 ph cm −2 s −1 in 2 − 7 keV) hard X-ray sources from the Chandra Source Catalog 2.0 (Evans et al. 2018) are shown (green stars), as well as the 90 (XMM-Newton) and 50 (NuSTAR) source regions (lime) and the respective elliptical regions used for background subtraction (yellow). The arrow (lime) indicates the direction to Sgr A*. Color bars indicate intensity in photons per pixel. The images have been smoothed using a 2D Gaussian kernel with standard deviation 3 pixels (XMM-Newton pn; pixel-size 4.3 ) or 5 pixels (NuSTAR FPMA; pixel-size 2.5 ).
significance in the 2 − 10 keV band of XMM-Newton and by NuSTAR from 10 − 20 keV. In the 20 − 79 keV band, the NuSTAR observation is dominated by background, and the core is not significantly detected (<3σ).
In addition to the central core and envelope, Fe Kα emission is detected at >5σ significance from four substructures within the diffuse region of the Sgr B2 in the projected plane. Two of these substructures coincide spatially with the cloud features previously identified by Terrier et al. (2018) as G0.66-0.13 and G0.56-0.11; Zhang et al. (2015) also detected G0.66-0.13 in hard X-rays in 2013. Here we additionally report two substructures, labeled G0.61+0.00 and G0.75-0.01, which were not detected in previous observations. Of these four substructures, only G0.66-0.13 and G0.61+0.00 lie within the NuSTAR field of view. These substructures, which are fainter than the core, are detected by NuS-TAR at 6.4 keV but not resolved above background in the higher energy bands. Figure 2 shows the background-subtracted spectrum of the central region of Sgr B2 as observed in 2018, overlaid with the best fits to three spectral models. All spectral fitting was performed using XSPEC software (Arnaud 1996). In this section, we detail the spectral fitting to models including a phenomenological model in Section 4.1, three XRN models in Section 4.2, and models of LECR-induced X-rays in Section 4.3.   The XMM-Newton data are binned with 3σ significance while the NuSTAR data are binned with 1.5σ significance. The best fit is shown in the solid lines. The contributions of the apec (dotted) and the nonthermal spectral components (dashed; ga, po, LECRe, and LECRp for the respective models) are also shown. All three models show satisfactory agreement with the data. The excess below ∼2.3 keV is compatible with foreground emission as also observed towards Sgr C (e.g. by Ryu et al. (2013)).

SPECTRAL ANALYSIS OF THE SGR B2 CORE
We extracted spectra from XMM-Newton in a 90 source region, which includes the cloud's core and part of the envelope, consistent with Zhang et al. (2015). We used a local background region, which includes any diffuse 6.4 keV emission from the larger GC region (Uchiyama et al. 2013). Constrained by stray light contamination (Section 2.2), we extracted NuSTAR spectra from the inner 50 and from FMPA only. Due to the smaller NuSTAR field of view, background subtraction was performed using an ellipse located within the diffuse region of Sgr B2. The spectral extraction regions are illustrated in Figure 1 and listed in Table 5.
Background subtraction is expected to account for the instrument background as well as diffuse emission from the GC region. Any faint point sources within the selected background region are also subtracted. Due to the different sky regions, flux from the diffuse region of Sgr B2 is subtracted from the NuSTAR but not XMM-Newton spectra, leading to an underestimation of the emission from the core, which is corrected in the next paragraph. No hard point sources were detected in either background region, so above 10 keV, the contribution from both point sources and the diffuse GC emission is expected to be small.
The XMM-Newton and NuSTAR spectra differ in source and background regions used for spectral extraction, in energy band, and in instrument characteristics. Here, the XMM-Newton data are used from 2 to 10 keV to constrain most spectral characteristics. The NuS-TAR spectrum is used from 10 to 20 keV to constrain the spectral index of the continuum, which is observed to be consistent with previous measurements. To facilitate a simultaneous spectral fit, we introduce a multiplicative factor relating the overall normalization of the NuSTAR spectrum to that of the XMM-Newton spectrum as a free parameter in the fit. We estimate a multiplicative factor of ∼0.16. Of this, a factor of 0.49 ± 0.05 in the relative normalization is attributed to the smaller NuS-TAR source region, where the factor, which depends on the distribution of emission in the source area, was calculated based on image analysis of the 6.4 keV band of XMM-Newton. Based on similar image analysis, the location of the NuSTAR background region within the diffuse region of Sgr B2 contributes a factor of 0.66 ± 0.11. Finally, decrease of up to ∼50% is anticipated due to the larger NuSTAR point spread function.
In all models discussed below, we use apec to model thermal emissions remaining after background subtraction, following Zhang et al. (2015). Other works used two apec components (Muno et al. 2004;Walls et al. 2016), where a cooler component at 1 − 2 keV accounts for diffuse GC X-ray emission while a warmer component at 6 − 8 keV accounts for unresolved point sources.
Here, we use a single apec, for direct comparison with results in Zhang et al. (2015). There were no significant differences in nonthermal model parameters between our reported results using the single apec model and fits using two apec components with fixed temperatures.
In Sections 4.2−4.3 we also consider fitted metallicity Z of the cloud relative to solar abundance Z as a metric for the physicality of a fit. In the CMZ, we expect Z/Z in the range of 1 − 2, based on previous measurements (Revnivtsev et al. 2004;Nobukawa et al. 2011;Jones et al. 2011). We nominally assume Z/Z = 2 in the fitting but consider 1 ≤ Z/Z ≤ 2 as reasonable. Table 2. Best-fit spectral parameters are shown for a joint fit of the 2018 XMM-Newton and NuSTAR observations, using the central 90 of Sgr B2 for XMM-Newton and central 50 of Sgr B2 for NuSTAR. We report flux parameters for the 90 region.

Parameter
Unit (257) 1.08 (255) 1.08 (255) 1.11 (257) Note-The goodness of fit is estimated by χ 2 ν and the number of degrees of freedom is given in parentheses. The errors represent 90% confidence. The fluxes and normalizations are for the 90 region. The multiplicative factor relates the flux from the 50 NuSTAR source region to the 90 XMM-Newton region. a The phenomenological model is given by wabs*(apec+wabs*po+ga+ga) in XSPEC. It is characterized by foreground absorption wabs, parametrized by the interstellar hydrogen column density NH (f ), and by internal cloud absorption parametrized by column density NH (i). The thermal apec component is characterized from metallicity Z/Z which could not be constrained and was thus fixed at 2, the plasma temperature kT , and total flux contribution of Fapec. The lines, ga, were fixed at 6.4 keV (neutral Fe Kα) and 7.06 keV (neutral Fe Kβ), with the Fe Kβ flux fixed at 0.15 of Fe Kα. The continuum was modeled by a pure powerlaw po characterized by the photon index Γ pl . We fixed Γ pl = 2 according to the best fit obtained from NuSTAR data alone, to avoid skewing by the higher-statistic XMM-Newton data. b The LECR models are given by wabs*(apec+wabs*LECR) in XSPEC, where LECR is the electron (LECRe) or proton (LECRp) XSPEC model (Tatischeff et al. 2012). The LECR models are characterized by the maximum path length Λ of particles in the cloud, as well as the minimum energy Emin of particles able to traverse the cloud, and the metallicity Z/Z of the cloud. In addition to these cloud parameters, the power law index s of the incident particle population is obtained in the fit. The fit with the LECRp model is reported under two conditions: first with the model parameters allowed to vary, and, second (labelled Z = 1) with parameters constrained to physical values. As described in the text, these LECR spectral fitting results do not clearly correspond to a physical ambient LECR scenario due to the assumptions inherent in the modeling and because the 90 region is not physically motivated as the X-ray production region for ambient LECRs. * Starred parameters were not allowed to vary in the fit.

Phenomenological Model
Throughout this paper, we use a phenomenological model to directly evaluate the 6.4 keV line flux. This model is detailed by Zhang et al. (2015) and given by wabs*(apec+wabs*po+ga+ga). The powerlaw continuum (po) and the neutral Fe Kα (6.4 keV) and Kβ (7.06 keV) lines (ga) expected in both the X-ray reflection and LECR scenarios are included explicitly. The model also accounts for thermal plasma (apec) emission that persists after background subtraction, as well as both intrinsic and foreground absorption (wabs). We use the wabs model, rather than updated models such as tbabs, to facilitate direct comparison with previous works.
We fix the line energies at 6.40 keV and 7.06 keV, noting that best-fit centroid energy of the Fe Kα line is at (6.40 ± 0.01) keV when it is allowed to vary. We also fix the line widths at 10 eV, i.e. much less than the en-ergy resolution of the instruments, and constrain the Fe Kβ normalization at Kβ/Kα = 0.15 (Murakami et al. 2001). We fit the apec plasma temperature but fix the apec metallicity Z/Z = 2 because it is not constrained by the data. The choice of Z/Z impacts the relative flux of the fitted apec and po parameters. The intrinsic and foreground hydrogen column densities, N H (i) and N H (f ), in wabs, are also fitted. After obtaining spectral index Γ ∼ 2 from the NuSTAR data only, consistent with previous measurements Zhang et al. 2015), we fix Γ = 2 in the combined fit. Following Zhang et al. (2015); Hailey et al. (2016), fixing Γ according to NuSTAR data prevents the higher statistics of the XMM-Newton data from skewing Γ based on an energy region where the power law is both degenerate with the thermal emission and more strongly absorbed.
Figure 2 (left) shows the spectral fit for the inner 90 , and the best fit model parameters are in Table 2. We obtained a satisfactory fit with χ 2 ν = 1.07 for 257 d.o.f. The best fit foreground column density, N H (f ) = 0.9 +0.2 −0.1 × 10 23 cm −2 , was higher than the expected value of 0.7×10 23 cm −2 to the GC, and the fitted intrinsic column density N H (i) = 4.6 +0.7 −0.6 × 10 23 cm −2 was comparable to the best fit of 5.0 ± 1.3 × 10 23 cm −2 found with NuSTAR for the same source region (Zhang et al. 2015). We do not expect a physical value for N H (i) in this case because it only represents an average over the region, rather than the complex scattering dynamics in the cloud. We note that all nonthermal parameters are consistent between the values reported here and those obtained in a fit with a two-apec plasma.

X-ray Reflection Nebula Models
We use three self-consistent XSPEC models of X-ray emission in the XRN scenario, including the MyTorus model (Murphy & Yaqoob 2009;Yaqoob 2012), the model developed by Walls et al. (2016), and the uniform Cloud REFLection of 2016 (CREFL16) model (Churazov et al. 2017a), to assess the compatibility of the 2018 emissions with an XRN origin. All three models produce acceptable fits, but most model parameters are poorly constrained. Therefore, we conclude only that the emission spectrum observed in 2018 is not inconsistent with a primarily XRN origin. Details of the fitting with the three spectral models are in Appendix A.

Low-energy Cosmic Ray Models
We use the LECRe and LECRp XSPEC models (Tatischeff et al. 2012) to understand if the spectral characteristics of the Sgr B2 core in 2018 are consistent with a LECR origin. However, we note that the physics of LECR diffusion into dense clouds is highly model dependent (see Section 7.2). The LECR (LECRe and LECRp) models were developed for the Arches cluster, which is smaller and less dense than Sgr B2 and features a stellar cluster that could accelerate CRs. Further, the 90 region is not physically motivated as the X-ray production region corresponding to an ambient LECR population. Therefore, the results of the model fitting in this section do not correspond to the most physical ambient LECR scenario. More robust limits on ambient LECR populations are related to the Fe Kα line fluxes from different cloud regions, reported in Section 6.
The XSPEC model is given by wabs*(apec+wabs*LECR), where the wabs and apec components account for foreground and internal absorption and any plasma emissions, as in Section 4.1. The LECR models assume a MC is bombarded by CRs from an external source whose spectrum follows a powerlaw with index s. The remaining parameters, including the path length Λ of CRs in the X-ray production (nonthermal) region, the minimum energy E min for a LECR to enter the X-ray production region, and the metallicity Z, are properties of the MC. The normalization N LECR describes the injected power dW/dt = 4πD 2 N LECR by LECRs from E min to 1 GeV, given distance D to the MC.
The fitted parameters for the LECRe and LECRp models are in Table 2 and fitted spectra are in Figure 2. Λ could not be constrained by the data and is frozen at 5 × 10 24 H-atoms per cm 2 in accordance with the column density through the core following Tatischeff et al. (2012); Zhang et al. (2015). The metallicity in the LECR model is that of the molecular gas, distinct from that of the apec component.
In the electron (LECRe) case, the fit is satisfactory with χ 2 ν = 1.08 for 255 d.o.f. The best-fit foreground N H (f ) = (0.9 ± 0.1) × 10 23 cm −2 and intrinsic N H (i) = 5.2 +1.2 −1.1 × 10 23 cm −2 column densities are consistent with previous observations. The best-fit plasma temperature is kT = 4.3 +1.1 −0.7 keV, and the cloud metallicity is Z = 1.9 +0.8 −0.4 Z , consistent with the expected range of 1 − 2 Z . The fit favors no lower cutoff on LECR energies in the cloud, with E min = 3.2 +27.7 −2.2 keV, and an elec-  tron spectral index of s = 3.2 +0.8 −0.7 . The fit normalization N LECR = 0.9 +2.0 −0.8 ×10 −6 erg cm −2 s −1 . For D = 8 kpc to Sgr B2, this corresponds to a limit (90% confidence) on the power of LECR electrons, dW/dt < 2.2×10 40 erg s −1 from the central 90 . However, this upper limit is of limited utility given the assumptions in the LECR model and considering that ambient LECRs are not expected to reach the core region of Sgr B2.
In the proton (LECRp) case, the fit statistic for the central 90 is also satisfactory, with χ 2 ν = 1.08 for 255 d.o.f. The best-fit foreground N H (f ) = (0.9 ± 0.1) × 10 23 cm −2 and intrinsic N H (i) = 5.0 +0.4 −1.0 × 10 23 cm −2 column densities are consistent with previous observations. The plasma temperature is kT = 4.3 +1.1 −0.7 keV. E min is completely unconstrained, while the fit favors a similar LECR proton spectrum as in the LECRe case, s = 2.9 +1.6 −1.2 , consistent with the 1.5 < s < 2 expected from diffusive shock acceleration (Tatischeff et al. 2012). The cloud metallicity is fitted as Z < 0.8 Z (90% confidence), which is inconsistent with the expected value of 1 − 2 Z , so this fit does not represent a physical scenario. However by fixing the metallicity at Z/Z = 1 and the proton power law index at s = 1.5, we obtain a similar quality fit with χ 2 ν = 1.11 for 257 d.o.f. and physical model parameters (see Table 2). In this case, N LECR = 1.7 +0.4 −0.2 × 10 −7 erg cm −2 s −1 . The corresponding upper limit on the LECR proton power dW/dt < 1.6×10 39 erg s −1 is subject to the same caveats as in the electron case above.  Figure 1 despite the change in color scale. In these exposure-corrected and continuum-subtracted images, the surface brightness of the core and envelope decreases over time. In contrast to this continued dimming, substructures within the diffuse region, including those identified in Section 3, brighten and dim from one observation to the next.
Here, we discuss this changing Fe Kα brightness and morphology from the cloud overall (Section 5.1) and from the substructures (Section 5.2). Details of the sky regions and spectral fitting are in Appendix B.1 and Table 5. All spectra were extracted using a local background region, such that the reported Fe Kα flux is in excess of the diffuse emission reported by e.g. Uchiyama et al. (2013). Figure 4 presents light curves of the Fe Kα emission from Sgr B2, illustrating the behavior of the diffuse, envelope, and core regions. We use a 6 region (not concentric with the core; defined in Table 5 and illustrated in Figure 3) to illustrate the behavior of the cloud over all. We also analyze the 90 region detailed in Section 4 to probe the behavior of the envelope, and a 10 region to probe the core.  The 6 region includes the Sgr B2 envelope and core and the bulk of the emission from the diffuse portion of the cloud. Unlike the full 9.9 indicated by the simplified model, it is compatible with all four XMM-Newton observations. The resulting light curve illustrates that the total flux from Sgr B2 has continued to decrease with time, by a factor of ∼3.5 since 2001. Relative to the 2012 level, the 2018 emission represents a (29 ± 8)% decrease, indicating that the primary XRN component was contributing to the total nonthermal emission from this region at least as late as 2012. We note that Khabibullin et al. (2021) reported the Fe Kα emission from a comparable sky region using more recent (2019) commissioning data from SRG/eRosita. However, while the 2019 emission level is consistent with the 2018 level reported here, a higher-statistics SGR/eROSITA observation would be required to constrain the light curve.

Time Variability of the Central Region
The 6 region directly corresponds to the sky region detected with the INTEGRAL observatory (which has a 6 resolution), facilitating comparison with the light curve of hard continuum emission reported by Kuznetsova et al. (2021).
The Finally, we show the light curve for the central 10 , which corresponds to the ∼15 half-power diameter of XMM-Newton together with the width of the core. Though the fitting is less significant due to the small source size, we observe a similar pattern of decreasing emission for the core as for the cloud overall. Based on the small relative flux contribution from the core, we conclude that behavior of the 90 region is driven by interactions in the envelope.
The true shape of the XRN light curve depends on the shape of the original flare and the details of the density profile of the cloud (Sunyaev et al. 1993). Here, we fit the intensity I of the emission from the central 90 as a function of the time t since January 1, 2001 as an exponential decrease with a constant offset: .
In context of the light curves, Figure 3 illustrates how the geometry of the emission from the core and envelope has evolved over time. While the 2001 map is brightest on the Sgr A* side of the envelope, the emission is more balanced by 2004, and by 2012 and 2018, the envelope emission is more extended on the opposite side of Sgr A*, as illustrated by the contour lines. This provides further indication that the initial flare from Sgr A* has already passed through some or all of the envelope. Within each substructure, we defined several 40radius circular regions (∼10-year light-crossing time) to illustrate the patterns of light that appear to be moving through the larger substructure. In each substructure, the circles are ordered from least negative declination (A, magenta, i.e. farthest from Sgr A* in the projected plane) to most negative, and the circle A brightens last. For G0.66-0.13, the light curves for circles B and C follow the same pattern as the parent structure, while circle A brightens in 2018, consistent with a flare originating at Sgr A* propagating through the cloud. We note that while these X-ray substructures were identified within the projected area of diffuse region of Sgr B2, we cannot exclude that they may correspond to other structures along the line of sight but outside of the Sgr B2 complex. Efforts to clarify the location of the substructures using line-of-sight velocity maps from the MOPRA 3 mm survey (Barnes et al. 2015) were inconclusive.

LOW-ENERGY COSMIC RAY LIMITS
Theoretical efforts to model the propagation of LECRs into Sgr B2 rely on simplified models of the gas distribution and cloud dynamics. Models are necessitated by the complex gas structure, illustrated in Figure 6. In this section, we have selected several sky regions (illustrated in Figure 6 and detailed in Table 5) that are compatible with the diffuse, envelope, and core components of Sgr B2 in both the simplified model and the observed hydrogen column density while also avoiding the bright substructures identified in Section 5.2. Table 3 presents the Fe Kα flux and surface brightness based on spectral fitting of these sky regions. Details of the data and fitting are in Appendix B.2 and Figure 7. The Fe Kα surface brightness from the representative diffuse, envelope, and core regions of the cloud are (2.9 ± 0.9), (5.7 ± 1.4), and (16 ± 4) × 10 −7 ph cm −2 s −1 arcmin −1 , respectively. Since the 2018 data are the dimmest deep observations of Sgr B2 in this band to-date, these are the current best upper limits on Fe fluorescence due to ionization by LECRs. Note-Data are reported based on the 2018 XMM-Newton observation of Sgr B2. The region boundaries are in Appendix B (Table 5). Regions are circular or annular, with given angular size in radius, unless otherwise specified. Errors and upper limits indicate 90% confidence. The corresponding spectral fits are in Figure 7. * The region selected from the diffuse portion of the cloud is an ellipse, chosen to fall within the diffuse region in both the simplified model and the observed nH and to avoid hard point sources. The reported flux is thus the flux from this region, rather than the total flux from the diffuse region. Due to the more limited field of view of MOS1, we calculated flux using MOS2 and pn only for this region. † Because the actual gas distribution in Sgr B2 is more complicated than the simplified model, two distinct sky regions were used to evaluate the envelope flux from Sgr B2, as shown in Figure 6. The annular region, Env. (0.5 − 2.2 ), represents the bulk of the envelope in the simplified Sgr B2 model. While the total flux measurement from this region may be useful, Sgr B2 has subdominant cores located within this annulus, and portions of this annular region have a column density more similar to the diffuse region, so the surface brightness should be interpreted with caution. On the other hand, the elliptical envelope region, Env. (ellipse), is a region with typical column density for the Sgr B2 envelope. While the flux for this region does not represent the total flux from the Sgr B2 envelope, the surface brightness is typical of the portions of the cloud with the intermediate column density associated with the envelope. 7. DISCUSSION 7.1. Implications for the X-ray Reprocessing Scenario The core and envelope region, previously the brightest part of the cloud, with multiple cores detected (Zhang et al. 2015), is now very faint in Fe Kα fluorescence, as illustrated in Figures 3 and 4. The brightest emission is restricted to the densest core. The pattern of flux decrease from 2013 to 2018 from the core, envelope, and cloud overall indicates that most or all of the major Xray flare previously driving the overall luminosity has passed through the cloud by 2018. The 2018 emission from the central 90 is probably driven by a component that is not decreasing over at least the past several years. Such stationary or quasi-stationary behavior is expected from both multiple scattering and LECRs and could also arise from primary X-ray reflection if the short and bright flare responsible for the time variability was immediately followed by a longer-duration plateau with 0.1 of the peak luminosity. Reflection of the short and bright flare component is not ruled out as a subdominant contributor if it has not fully exited the 90 region. Reflection of secondary flares in diffuse structures along the line of sight could also contribute. The total decrease in the Fe Kα emission from the 6 region from 2004 to 2018 is consistent with the total decrease in the 30 − 80 keV continuum observed by INTEGRAL over the same sky region. The best-fit IN-TEGRAL light curve is a linear decay before 2011 and a constant level thereafter (Kuznetsova et al. 2021). In comparison, Fe Kα light curve reported here was decreasing at least as late as 2012, but due to limitations in the number and significance of the data points, the difference between the two light curves is not significant. A Fe Kα light curve that is consistent with that of the hard continuum indicates that the emission mechanisms for the two energy scales are related. This is expected, since the higher absorption at 6.4 keV is partially compen-  Figure 6. The hydrogen column density as measured by Herschel (Molinari et al. 2011) is shown in log scale and illustrates the complexity of the Sgr B2 structure compared to the simple model (yellow, 9.9 diffuse, and white, 2.2 envelope). The regions used in Table 3 are shown in magenta as the ellipses representing the clean diffuse and envelope regions, and with the 30 circle. The background region is also shown (green ellipse), alongside the brightest hard point sources from the Chandra Source Catalog (green stars). photons above 7.1 keV. If the emission is dominated by multiple scattering, the actual degree of correlation depends on geometry, optical depth and metallicity of the cloud (Sunyaev & Churazov 1998;Odaka et al. 2011). Kuznetsova et al. (2021) suggest multiple scattering of X-rays from the primary external flare as a probable origin of the 30 − 80 keV emission after 2011. Multiple spectral handles can clarify the role of multiple scattering. First, this signal is expected to be more pronounced in the morphology of the hard X-ray continuum than the fluorescent lines, as the absorption cross section is larger than the scattering cross section 10 keV (Odaka et al. 2011;Churazov et al. 2017a). Against a backdrop of fading emission, the 2013 NuSTAR detection (Zhang et al. 2015) of multiple cores above 10 keV suggests that multiple scattering already played an increasingly important role. Continued decrease of the Fe Kα emission once the 30 − 80 keV emission has reached a stationary level, as hinted by the data, would support the multiple scattering origin of the 30 − 80 keV emission. A future NuSTAR observation less contaminated by stray light would constrain the hard-continuum light curves from the core and envelope, allowing the cleanest probe of the multiple scattering scenario. While the 2018 obser-vation was optimized to minimize stray light contamination, cleaner NuSTAR data could be obtained using a contiguous deep observation of a nearby off-source region to measure the stray light background, as in the 2013 observation (Zhang et al. 2015).
In the Fe Kα line, a Compton shoulder feature due to multiple scattering of fluorescence photons is expected on the low-energy side of the line complex (Odaka et al. 2011). While the Compton shoulder is not resolvable with the energy resolution of XMM-Newton, the centroid of the 6.4 keV line is expected to shift toward lower energies as the Compton shoulder becomes a more important contributor to the overall line flux (Khabibullin et al. 2021). The fitted line centroid, 6.40 ± 0.01 keV in 2018 for the central 90 , is consistent with a constant value over the four XMM-Newton observations for both the 6 and 90 regions. Future high-resolution X-ray spectrometers including XRISM/Resolve (Ishisaki et al. 2018) and Athena/X-IFU (Barret et al. 2018) could resolve the Compton shoulder, precisely measuring the relative contribution of multiple scattering (Odaka et al. 2011).
The morphological and brightness variation from several X-ray substructures (reported in Section 5.2) external to the Sgr B2 envelope also reveal implications for the timescale and number of the external X-ray flaring events. The small bright regions, with their relatively short light crossing times, have a faster timescale of emission decrease from X-ray reprocessing and thus better reflect the timescale of the external source than the Sgr B2 envelope . In particular, G0.56-0.11, which was reported as brighting up in 2012 , is even brighter in 2018, with a significant morphologic change. While the 2012 emission is centered on clump B (see Figure 5), the 2018 emission is centered ∼ 13 light years (projected distance) away, in clump A. The light curve of clump B gives an upper limit of ∼14 years for the timescale of the flare illuminating this region. Future observation by XMM-Newton could further constrain the timescales of these flares based on the future behavior of the "A" clumps from each substructure.
Similar light curves to G0.56-0.11 clump B are observed in G0.66-0.13 clumps B and C, suggesting that these two substructures may be illuminated by the same flaring event, if they are a similar distance from Sgr A*. Unfortunately, without knowledge of the line-of-sight positions and detailed n H distributions of these structures, we cannot clearly make this claim. Clavel et al. (2013), Chuard et al. (2018), andTerrier et al. (2018) provided evidence for a minimum of two illuminating events propagating through the CMZ. In the case that these substructures are illuminated by the same event as the core, if we assume Sgr B2 to be at least 50 pc in front of Sgr A* (following Walls et al. (2016); Yan et al. (2017), and references therein), the substructures could be 60 pc behind Sgr B2, farther than the spatial extent of cloud. Therefore, if these substructures are linked to Sgr B2, they are illuminated by a secondary event.

Implications for Low-energy Cosmic Rays
This section discusses the upper limits on the Fe Kα emission from different cloud regions presented in Section 6 and Table 3 in the context of physical models of LECR interactions in Sgr B2.
CR transport is modeled as a diffusive process modified by the effects of elastic and inelastic collisions, energy loss via ionization and excitation, and energy loss via bremsstrahlung and synchrotron radiation in the surrounding medium. While transport into and within MCs is model dependent (Gabici et al. 2007;Tatischeff et al. 2012;Gabici 2013;Dogiel et al. 2015;Morlino & Gabici 2015;Owen et al. 2021), the low relative rates of hydrogen ionization observed within dense MCs indicate that LECRs do not freely traverse these structures (Oka et al. 2005;van der Tak et al. 2006;Indriolo & McCall 2012;Dogiel et al. 2015).
Using a simplified model of Sgr B2, Dogiel et al. (2015) calculate that LECR propagation in the envelope 1 is best described by diffusion on turbulent magnetic fluctuations (Dogiel et al. 1987;Istomin & Kiselev 2013). Meanwhile, diffusion is negligible in the diffuse region, where fluctuations are small. Considering ionization and excitation losses, protons (electrons) with kinetic energy E 20 MeV (1 MeV) traverse the diffuse region to reach the envelope, where they are absorbed within 0.1−0.3 pc (Dogiel et al. 2009a). The ambient LECR proton spectra derived by Dogiel et al. (2015) would deposit the bulk of their energy in the Sgr B2 envelope.
In the LECR proton case, Dogiel et al. (2015) use the hydrogen ionization rates in Sgr B2 and the surrounding environment to estimate the intensity I 6.4 ≈ (3−5) Z Z × 10 −6 ph cm −2 s −1 of the LECR-induced Fe Kα emission from the Sgr B2 complex. The range of 3 − 5 depends on the details of the ambient LECR spectrum and the gas distribution of the cloud.
The measured Fe Kα flux of (10.8 ± 1.2) × 10 −6 ph cm −2 s −1 from the envelope (0.5 − 2.2 ) region in 2018 (Table 3) is comparable to the calculation by Dogiel et al. (2015) if Z/Z = 2. Unless Z/Z > 2, ambient LECR protons in this model cannot explain all of the Fe Kα emission. However, propagation of ambient LECR protons could contribute >50% of the 2018 Fe Kα emission. Kuznetsova et al. (2021) disfavor LECR proton propagation as the sole source of the 30 − 80 keV emission from 2011 to 2019 on the basis of the high overall flux compared to the observed hydrogen ionization rate. As detailed in Section 7.1, the 30 − 80 keV band could have more substantial contributions from multiple scattering compared to the Fe Kα line, and a portion of the flux in the 6 region of INTEGRAL is due to X-ray reflection from substructures. With improved characterization of the flux contribution expected from multiple scattering, enabled by future high-resolution observations of the 6.4 keV line and the hard X-ray flux as discussed in Section 7.1, the LECR contribution could be more tightly constrained even as it coexists with the emission from multiple scattering.
Previous work has demonstrated that an ambient population of LECR electrons cannot explain the full hydrogen ionization rate in Sgr B2. The X-ray continuum from propagation of ambient LECR electrons is expected to have Γ ∼ 1, harder than the Γ ∼ 2 expected from LECR protons or X-ray reflection and observed from Sgr B2. Given the population of ambient LECR electrons that could explain the observed hydrogen ionization, the hard continuum should have been observable as early as 2009 (Dogiel et al. 2015;Zhang et al. 2015;Kuznetsova et al. 2021). While ambient LECR electrons are excluded as the sole origin of hydrogen ionization in Sgr B2, they may still contribute to the nonthermal Xray emission observed in 2018 or in the future. The sensitivity of the 2018 observations to LECR electrons is restricted by the limitations of the NuSTAR observation. With spatial resolution for X-rays up to 79 keV, a future NuSTAR observation less severely contaminated by stray light could resolve this ambiguity.
We additionally considered the annihilation or decay of dark matter as a source of the LECRs responsible for the hydrogen ionization and the nonthermal X-ray emission, concluding that any contribution from dark matter would be small compared to the LECR electron population needed to explain the X-ray emission (details in Appendix C).

CONCLUSION
The X-ray features of Sgr B2 provide a window into past energetic activity of Sgr A*, via X-ray reprocessing in the cloud, and to the GC LECR populations, via ionization and Bremsstrahlung processes.
In this paper, we have presented the 103 ks and 149.2 ks observations of Sgr B2 taken jointly in 2018 by the XMM-Newton and NuSTAR X-ray telescopes, respectively. These data show that the the 2018 Fe Kα emission from the central region is 0.8 ± 0.2 of the 2012 level reported by Zhang et al. (2015), consistent with significant contribution from a stationary component. Meanwhile, they also reveal new brightening substructures, which, assuming they correspond to clumps within the Sgr B2 complex, are illuminated by one or more secondary external flares. Based on both the Fe Kα light curve and the spectral analysis, the 2018 emissions from the central region are consistent with arising primarily from X-ray reprocessing, with possible contributions from the tail of the original flare, multiple scattering albedo, and secondary flares, or with arising primarily from LECR interactions. Thus, the flux levels presented in Table 3 represent best upper limits on fluorescence from LECR interactions within different cloud regions. The Fe Kα emission observed from the Sgr B2 envelope is comparable with expectation from the lowenergy cosmic proton population that would simultaneously explain the observed hydrogen ionization rates in the model of Dogiel et al. (2015).
Future observations of Sgr B2 by XMM-Newton will further constrain the Fe Kα light curves from the envelope, core, and the diffuse substructures, clarifying the origins of the 2018 emission and the duration of the external flare(s) illuminating the bright substructures. Meanwhile, further observations resolving the dense cores can clarify the contribution of multiple scattering and facilitate correspondingly more precise limits on the contributions of LECRs. Two regimes are of particular interest. For hard X-rays above 10 keV, scattering is the dominant photon process, and a NuSTAR observation less contaminated by stray light could detect deviations in the light curve of the densest cores relative to the envelope and diffuse cloud regions. For the Fe Kα line, future high-resolution spectrometers could directly resolve the line features, including the Compton shoulder, expected in multiple scattering. If further decrease of the Fe Kα emission from the envelope is observed, or if a significant portion of the 2018 Fe Kα emission level is definitively attributed to multiple scattering, the LECR proton model of hydrogen ionization (Dogiel et al. 2015) would begin to be constrained.

Facilities: XMM-Newton -The X-ray Multi-Mirror
Mission, NuSTAR -The Nuclear Spectroscopic Telescope Array Mission Software: XSPEC (Arnaud 1996), ESAS (Snowden et al. 2008), heasoft (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014), ciaotools (Fruscione et al. 2006), MyTorus (Yaqoob 2012) SAOImage ds9 (Joye & Mandel 2003), Matplotlib (Hunter 2007), IPython (Perez & Granger 2007), and Astropy (Astropy Collaboration et al. 2013 (257) Note-The goodness of fit is estimated by χ 2 ν and the number of degrees of freedom is given in parentheses. The errors represent 90% confidence. The fluxes and normalizations are for the 90 region. The multiplicative factor relates the flux from the 50 NuSTAR source region to the 90 XMM-Newton region. a Fitting with the phenomenological model (wabs*(apec+wabs*po+ga+ga) in XSPEC) is detailed in Section 4.1 and Table 2. The model parameters are repeated here for comparison with the XRN spectral models. b The MyTorus (Murphy & Yaqoob 2009) model is given by wabs*(apec + MYTS + MYTL) in XSPEC, where MYTS and MYTL are two models representing the scattered continuum and fluorescent line emissions, respectively. MyTorus is parameterized by the photon index Γ of the external X-ray population, the internal column density NH (i), and an overall normalization N , with all parameters coupled between MYTS and MYTL for consistency. c The Walls model is given by wabs*(apec+Walls) in XSPEC (Walls et al. 2016), where Walls is parametrized by the photon index Γ of the external X-ray population, the internal column density NH (i), and the inclination angle θ of the cloud relative to the X-ray source, as well as an overall normalization factor N . d The CREFL16 model is given by wabs*(apec+CREFL) in XSPEC (Churazov et al. 2017a). It is parametrized by metallicity Z/Z , θ, ΓXRN , and τT = 2σT NH (i), where the Thomson cross section σT = 6.65 × 10 −25 cm 2 , as discussed in the text. Here we report NH (i) directly for comparison with the other models. * Starred parameters were not allowed to vary in the fit. spectra are not inconsistent with arising primarily from X-ray reflection. Table 4 summarizes the fitting results.
Previously, Zhang et al. (2015) modeled X-ray reflection from Sgr B2 using the MyTorus X-ray Reprocessing Model (Murphy & Yaqoob 2009;Yaqoob 2012) developed for X-ray reflection from toroid structures surrounding Compton-thick AGN. MyTorus is parametrized by the cloud hydrogen column density N H (i), the spectral index Γ of the external X-ray source, and an overall normalization. The model is additionally parametrized by an inclination angle of the torus relative to the source, which we fix at 0 • to most closely describe the situation of a spherical cloud. MYTorus has additional limita-tions beyond geometry, detailed by Zhang et al. (2015), when applied to model MCs as XRN, and in particular it assumes reflection in a medium with metallicity Z = Z and uniform density. MYTorus is implemented as a series of models in XSPEC. To best approximate the MC scenario, we used two MyTorus components, describing the scattered continuum (MYTS) and the iron fluorescence lines (MYTL), with all parameters coupled between these components so that they were consistent with arising from the same incident X-rays. For consistency with Zhang et al. (2015), we used the models with reflected X-rays up to 500 keV and an energy offset of +40 eV for the fluorescent lines. Thus the XSPEC model is wabs*(apec+MYTS+MYTL), where the MYTorus components are mytorus scatteredH500 v00.fits and mytl V000010pEp040H500 v00.fits, and the thermal plasma apec parameters are as in Section 4.1. The fit quality is marginally acceptable with χ 2 ν = 1.29 for 257 d.o.f. Additionally, both the intrinsic and foreground column densities are fitted much higher than the physical expectation, and the apec component absorbed a flux somewhat larger than expected based on the phenomenological model fit, possibly because it was compensating for low compatibility of the myTorus model with the data. Thus, the 2018 data are not well suited for the MyTorus model with physical parameters. However, due to the limitations of the myTorus model in describing Sgr B2, we do not exclude the XRN scenario based on this result. The Monte Carlo model developed by Walls et al. (2016), hereafter the Walls model, was developed specifically to model MCs as XRN with the particular application of constraining the location of Sgr B2 relative to Sgr A*. Like MyTorus, Walls depends on intrinsic absorption N H (i) of the cloud and spectral index Γ of the external X-ray source. It is additionally parametrized by the angle θ of the spherical cloud relative to the external source, such that θ = 0, which is unphysical given the non-zero projected distance from Sgr B2 to Sgr A*, corresponds to the cloud situated in line between the source and the observer, and the true distance from the X-ray source to the cloud is ∼100 pc/ sin θ. The position θ of the cloud relative to the source strongly impacts the XRN spectrum, as absorption patterns and effective cloud thickness differ dramatically between X-rays reflected to observers at different θ. There are several Walls models with varying cloud metallicities and density profiles. We chose the model with Z = Z for direct comparison with spectral fitting by Walls et al. (2016) but note that with these data the metallicity does not significantly change the spectrum below 10 keV, apart from the overall normalization of incident X-rays. We used a uniform, rather than Gaussian, cloud density profile, which gave a better fit and was found to be more suitable by Walls et al. (2016). Thus the model is wabs*(apec+Walls) where the Walls component is Walls et al 2016 Uniform 1.0.fits.
The fit was satisfactory, with χ 2 ν = 1.09 for 256 d.o.f. We found the best fit θ = 28 +10 −16 , significantly lower than the best fit of θ = 64 +8 −7 found for the bright state (Walls et al. 2016). However, since these data do not represent X-ray reflection in the bright state, we caution against overinterpretation of the result. The foreground column density was consistent with previous results at N H (f ) = 0.9 +0.1 −0.2 × 10 23 cm −2 , and the fitted intrinsic column density was on the high side of previous results, at N H (f ) = 7.9 +3.7 −2.1 × 10 23 cm −2 . While the source XRN spectral index was not constrained, the flux and temperature of the apec component were consistent with expectation, and the data are consistent with XRN origin according to the Walls model.
Finally, the uniform Cloud REFLection of 2016, CREFL16, model (Churazov et al. 2017a) treats the case of XRN in a uniform spherical cloud using the same geometry as Walls et al. (2016) but with additional fluorescence lines (beyond Fe) and treating metallicity Z/Z as a free parameter. In addition to Z/Z , the model parameters include the inclination cos θ to the illuminating source, the spectral index Γ of the illuminating source spectrum, and the optical depth of Thomson scattering, τ T = σ T n H R where σ T ∼ 6.65 × 10 −25 cm 2 is the Thomson scattering cross section, n H is the hydrogen density, and R is the cloud radius, such that N H (i) = 2n H R = 2τ T /σ T . In Table 2 we report θ and N H (i) for direct comparison with other models. Because the CREFL16 fit favored the unphysical result θ = 0, we fixed θ = 64 • as in Walls et al. (2016). Because the best fit gave Z/Z < 1, which is unphysical, we froze the metallicity at solar abundance, the lowest abundance consistent with observation. The fit was acceptable with χ 2 ν = 1.16 for 257 dof., but the apec component absorbed an order of magnitude more flux than expected, suggesting as with Walls that the CREFL16 model is not best suited for these data. Meanwhile, the apec and absorption parameters were consistent with those in the Walls fit.
For all three self-consistent XRN models, χ 2 ν was acceptable, but most model parameters were poorly constrained. In the CREFL16 model, the best fit angle to Sgr B2 was θ = 0, but the fit was still acceptable with the constraint θ = 64 • . While θ = 0 does not reflect the position of Sgr B2, the preference for θ = 0 may reflect a situation in which X-rays observed from Sgr B2 in 2018 are produced after the external X-ray front has partially or fully passed through the cloud. In this case, the external X-rays still being reflected have already traversed a cloud depth larger than the typical reflection depth in the bright state, and the effects of multiple scattering may also become important. The fits in this section illustrate that the 2018 flux from the core of Sgr B2 can be described by the XRN spectral models. However, even if X-ray reflection dominates the 2018 flux, the XRN models do not account for time-dependent effects as the external front passes through the cloud, so we cannot use these results to constrain the external X-ray spectrum.
B. SPECTRAL FITTING OF SGR B2 REGIONS Table 5 lists the sky regions used in the analysis of 1) time variability of the Sgr B2 structures in Section 5 (fitting details in Section B.1) as well as 2) LECRs in different portions of the cloud in Section 6 (fitting details in Section B.2). Circular and annular regions centered on the Sgr B2 core are not listed in the table, but are centered at RA=17 h 47 m 19.992 s , Dec=−28 • 23 07.08 .

B.1. Analysis of Sgr B2 Time Variability
The light curves in Section 5 (Figures 4 and 5) were produced using spectral fits to data extracted from four different XMM-Newton observations, extracted using a local background region selected for each observation. The source regions, and the background region for each observation, are given in Table 5. To produce each light curve, we first selected the observation when that substructure was brightest. If the 2001 observation was brightest, as for the regions in Figure 4, we selected the 2004 observation due to the low statistics of the 2001 observation. We performed a spectral fit to the phenomenological model as in Section 4.1, except that we froze N H (f ) = 0.9×10 23 cm −2 for consistency between regions and there was no multiplicative factor in the absence of NuSTAR data. Because the structures are too dim in some epochs to constrain a good fit, we then fixed all model parameters other than the normalizations before fitting the spectra from the other observations. Though constraints based on the large size of Sgr B2, hard point sources, and different focused fields of view necessitated use of different sky regions for background subtraction in each observation (see Table 5), cross checks using different background regions in the same observation yielded consistent Fe Kα flux measurements. Figure 7 shows the spectral fits with the phenomenological model used to obtain the Fe Kα fluxes reported in Table 3. Because of the restricted NuSTAR field of view, these data were based on fits to the 2018 XMM-Newton observations only, with spectra extracted using the source and background regions in Table 5. Spectral fitting with the phenomenological model is as in Section 4.1, except that we froze N H (f ) = 0.9×10 23 cm −2 for a fair comparison between regions and there was no multiplicative factor in the absence of NuSTAR data. The resulting fits were satisfactory as evaluated by χ 2 ν . The Fe Kα flux was evaluated using the cflux command in XSPEC. While the choice of N H (f ) introduces a systematic uncertainty into the total Fe Kα flux calculation, it is expected to be small compared to the sta-tistical uncertainty, as a 10% change in N H (f ) produces a ∼1.5% change in the fitted line flux at 6.4 keV. The surface brightness in Table 3 was calculated using the region area obtained by image analysis (i.e. accounting for missing pixels).

C. LIMITS ON DARK MATTER ANNIHILATION
The annihilation of dark matter has been of interest as a possible source of GC CRs (Murgia (2020) and references therein) including GC LECRs (e.g. Linden et al. (2011)). Here we consider the case of LECRs produced within Sgr B2 in the annihilation of dark matter as the nonthermal source responsible for the observed Fe Kα emission.
Any model of DM annihilating to low-energy electrons or protons would produce a LECR population within GC MCs and thus induce Fe Kα fluorescence. In contrast to the ambient LECR scenario of Section 7.2, LECR illumination would be expected throughout the Sgr B2 core and envelope. Considering the production of LECR from dark matter, the generic dark matter annihilation rate per unit volume is: where the thermally averaged dark matter annihilation cross section σv and dark matter mass M DM parametrize the model, and the local dark matter energy density ρ is estimated as ρ ∼ 100 GeV/cm 3 in Sgr B2, given ∼100 pc distance from the GC (Linden et al. 2011).
To estimate an overly optimistic scenario for the dark matter origin of LECR electrons with GC MCs, we assumed that all annihilation energy is transferred to electron kinetic energy. We also assume that all electrons produced in the core and envelope are absorbed within those regions, which is reasonable considering the diffusive dynamics of LECRs within the envelope of Sgr B2 in Section 7.2 and is explicitly treated by Gabici (2013). We considered the 90 sky region surrounding the core of Sgr B2, which represents a spherical volume of ∼165 pc 3 given a distance to Sgr B2 of ∼7.9 kpc. Comparing our fitted maximum power of LECR electrons in the central 90 of Sgr B2, dW/dt < 2.2 × 10 40 , to the annihilation rate in Eq. (C1), we calculated that a dark matter model with σv larger than 3 × 10 −22 cm 3 s −1 (M DM ∼ 500 keV) to 6×10 −19 cm 3 s −1 (M DM ∼ 1 GeV) is necessary to produce a LECR electron population sufficient to explain the bulk of the Fe Kα emission observed from Sgr B2.
This minimum possible σv is several orders of magnitude larger than both the existing best limits on σv for Table 5. Sky regions used for spectral extraction in Section 5 (upper), in Section 6 (middle), and for local background subtraction throughout this work (lower). Note-All coordinates are in terms of right ascension (R.A.) and declination (Dec.) using the J2000 system. Circular and annular regions centered on the Sgr B2 core are not listed in the table but assume the Sgr B2 complex is centered on RA=17 h 47 m 19.992 s , Dec=−28 • 23 07.08 . * No single suitable region located outside of the spatial extent of Sgr B2 and unaffected by hard point sources was also compatible with the field of view of all four XMM-Newton observations. Thus, multiple background regions were used. XMM-Newton Background A was used for all analyses with the 2018 data (0802410101). The alternate XMM-Newton Background B was used for the 2012 and 2004 observations (0694640601 and 0203930101) while XMM-Newton Background C was used for the 2001 observation (0112971501).
weakly interacting massive particle (WIMP) dark matter in this mass range annihilating to all visible products and the thermal relic annihilation cross section predicted by theory (Leane et al. 2018). Because this calculation assumes that all dark matter annihilation energy is transferred to the cloud via ionization and excitation by CRs, it is perfectly general to the case of hidden sector dark matter models as well. As with the WIMP case, existing limits on σv for hidden sector dark matter (Elor et al. 2016) are several orders of magnitude smaller than compared to a dark matter model necessary to explain the fluorescence in Sgr B2. Thus, annihilation of dark matter cannot produce a LECR electron population ca-pable of producing the observed Fe Kα fluorescence in Sgr B2.  Figure 7. The 2018 spectral fits that yield the data presented in Table 3, are shown. Background-subtracted spectra are shown for the Diffuse (upper left), Envelope (ellipse, upper right), Envelope (0.5 − 22 , lower left), and Core (0.5 , lower right) regions (given in Figure 6), fitted with the phenomenological model. The models are fitted simultaneously to data from the XMM-Newton MOS1 (black, 2 − 10 keV) MOS2 (red, 2 − 10 keV) and pn (green, 2 − 7.8 keV). The Diffuse region is not compatible with the MOS1 field of view. The data are binned with 3σ confidence. The best fit is shown in the horizontal lines. The contributions of the apec (dotted) and the nonthermal spectral components (dashed; ga and po) are also shown. Deviations between the data and the model are most pronounced at low energies.