© The Authors 2023

1 Introduction

There is growing evidence that the Galactic center (GC) plays an active role in the acceleration of cosmic rays (CRs). Very-high-energy gamma-ray observations (0.1–100 TeV) of the central molecular zone (CMZ) by the High Energy Stereoscopic System (HESS; Aharonian et al. 2004, 2009; H.E.S.S. Collaboration 2016, 2018a), the Major Atmospheric Gamma-Ray Imaging Cherenkov (MAGIC; MAGIC Collaboration 2020), and the Very Energetic Radiation Imaging Telescope Array System (VERITAS; Adams et al. 2021) all support a particle acceleration scenario. In these observations, extended emission correlates with the gas morphology in the CMZ, suggesting that a central and continuous source of protons with energies up to 1 PeV diffuses out and collides with protons in the ambient gas, producing neutral pi-mesons (π⁰) that decay into observable gamma rays (H.E.S.S. Collaboration 2016). The Fermi bubbles (FBs), two ~10 kpc-scale structures observed to be emanating from the GC by the Fermi Large Area Telescope (Fermi LAT) at 0.1–100 GeV energies (Su et al. 2010; Ackermann et al. 2014), are additional evidence of the particle acceleration. To explain them, either inverse Compton emission produced by high energy electrons or the proton-proton interaction just mentioned can be invoked, both corresponding to particle acceleration from the GC (Ackermann et al. 2014).

Currently, the Galactic sources of particles with ≈1 PeV energies, so-called PeVatrons, are still under debate (Aharonian et al. 2019). Supernovae have been proposed as the most important Galactic PeVatron, where CRs are impulsively accelerated on the shock front created by its supersonic outflow (Berezhko 1996; Hinton & Hofmann 2009). Additionally, stellar winds of Wolf-Rayet (WR) stars in compact and massive clusters of young stars could persistently accelerate CRs due to the collective effect of the shock fronts created by the entire WR population (Klepach et al. 2000). Finally, accretion onto a compact object could produce high energy particles (Hinton & Hofmann 2009). For the CMZ gamma-ray emission, CRs could be accelerating in the vicinity of the supermassive black hole Sgr A* (H.E.S.S. Collaboration 2016), within the Arches (AC), Quintuplet (QC), and nuclear (NSC) clusters of young massive stars (Aharonian et al. 2019), or due to an unresolved population of millisecond pulsars (Guépin et al. 2018). Despite many possibilities, the CR origin in the GC is still unknown. Our preliminary results show that the source is inside the NSC if the CMZ is a continuous structure, or from both the NSC and the supernova Sgr A East if the CMZ has an inner cavity (Scherer et al. 2022, hereafter Paper I).

Additionally, CR transport within the CMZ is not completely understood either. CR transport in the Galactic disk is known to be dominated by diffusion, where the diffusion coefficient is constrained by CR abundance and their fragmentation timescale (Longair 2011; Gaisser et al. 2016). However, in highly dense zones, such as the CMZ, it is theoretically expected (Ormes et al. 1988) and constrained from gamma-ray observations (Abeysekara et al. 2017; Aharonian et al. 2019) that diffusion should be two orders of magnitude lower than in the Galactic disk. The main reason is the extreme high-turbulent magnetic fields associated to the high density regions, which create a slow diffusion. Also, the outflow velocities of the FBs (Bordoloi et al. 2017) suggest that advective transport toward the Galactic poles may be important for the propagation of CRs within the CMZ.

The main features of the ~ 1 PeV CRs accelerated within the CMZ have been constrained indirectly via gamma-ray flux measurements and considering the gas column density (H.E.S.S. Collaboration 2016, 2018a). The gamma-ray luminosity emitted by proton-proton interaction is directly proportional to the ambient gas density (Aharonian 2004; Longair 2011), but the CMZ 3D morphology has not been completely determined yet (Sofue 2022; Henshaw et al. 2023). From molecular rotational line emission, for example CO lines (Oka et al. 2012) and CS lines (Tsuboi et al. 1999), different maps of the CMZ column density have been obtained, which have been used as input for hydro-dynamical simulations, to quantitatively compare the emission and absorption lines, or to search for the best orbits to explain the observed radial velocities and 2D morphology. From those analyses, a variety of CMZ shapes have been inferred by different studies, including two spiral arms (Sofue 1995; Ridley et al. 2017), a twisted elliptical ring (Molinari et al. 2011), an open elliptical stream (Kruijssen et al. 2015), either an elliptical ring if the star formation is minimal or a fragmented ring if the star formation is intense (Armillotta et al. 2019, 2020), or a bar-like structure (Sawada et al. 2004; Yan et al. 2017).

In Paper I, we explored the impact on the CR analysis of the CMZ having or not a large inner cavity with respect to the line of sight distribution. Using a simple CR model consistent with the H.E.S.S. Collaboration (2016) results, we found that a disk-like gas distribution reproduces the gamma-ray observations considering only stellar winds from WR stars in the NSC. However, a structure with a significant inner cavity requires a composite CR population, coming from the same stellar winds and the supernova Sgr A East. In this paper, we aim to reproduce the HESS observation of the CMZ considering more realistic CR dynamics and CR sources with injection rates derived from their observed kinetic energies. We compute the CR transport considering a differentiated diffusion (inside and outside the CMZ) and polar advection, for both gas distributions proposed in Paper I, and we explore the CR contribution from the NSC, AC, QC and Sgr A East. Additionally, we show that upcoming observations with the Cherenkov Telescope Array (CTA; CTA Consortium 2019) will verify the proposed transport model and constrain its ingredients.

2 Methodology

We computed a CR transport model using a Monte Carlo method. We considered a mix of continuous CR sources, namely the NSC, AC and QC, plus impulsive CR injection from Sgr A East. We modeled transport across two diffusion zones – low diffusion within the CMZ, and regular Galactic-disk diffusion outside the CMZ. Additionally, we added transport by polar advection. We computed the CR dynamics and calculated the gamma-ray luminosity in a 3D grid, using the same two gas distributions considered in Paper I – a CMZ with and without an inner cavity. Finally, we estimated the average CR energy indirectly to find the model that most closely matches the HESS measurements, and proposed observations with CTA to confirm and fine-tune the models.

2.1 Cosmic-ray model

We modeled the CR dynamics using the standard diffusion equation for protons (Aharonian 2004; Strong et al. 2007), neglecting reacceleration as its energy gain is much slower than regular acceleration in shock fronts (Longair 2011; Blasi 2013), (1)

where ψ(r, p, t) is the CR distribution, r is the position vector, p is the total CR momentum at position r, t is the time, D(E_p) is the diffusion coefficient, E_p is the energy of ultrarelativistic particles, V is the advection velocity, and Q(r, p, t) is the CR source term. In our model, we considered only particles of E_p = 1 PeV, so we transferred all energy of CRs above 10 TeV to these mono-energetic particles. We simplified the CR intrinsic spectrum because the solution of Eq. (1) considered by HESS was obtained with mono-energetic particles subject to standard diffusion. In this way, we can easily compare our new results to those by HESS and also to our own from Paper I. When assessing future CTA observations, in Sec. 3.2, we distributed the energy across the relevant domain.

2.1.1 Sources

We considered three continuous and one impulsive source, with CR injection rates derived from their observed kinetic energies. Additionally, we included another impulsive source whose properties are not well-constrained observationally.

For the continuous sources, which correspond to stellar clusters with colliding WR winds, we defined the injection rate as Q_con = A₀(E_p/E₀)^−Γg(t), where A₀ is the CR injection rate, E₀ is the energy normalization, Γ is the spectral index of the CR intrinsic spectrum, and g(t) is a Heaviside function in the time interval t_age < t, where t_age is the source age. A₀ was chosen by assuming that a canonical 1% of the stellar winds’ kinetic power goes into acceleration of CRs above 10 TeV (Aharonian et al. 2019) and (E_p/E₀)^−Γ = 1 for mono-energetic particles with all CR energy above 10 TeV. The proposed continuous CR sources are: NSC, modeled as a injection rate of A₀ = 3.4 × 10³⁶ erg s⁻¹, considering that the total mass-loss rate of its WR population is ≈10⁻³ M_⊙ yr⁻¹ and their typical stellar wind velocity is ≈1000 km s⁻¹ (see Table 1 of Cuadra et al. 2008, based on Martins et al. 2007); AC, modeled as an injection rate of A₀ = 2.0 × 10³⁶ erg s⁻¹, considering that the total mass-loss rate of its WR population is ≈2.0 × 10⁻⁴ M_⊙ yr⁻¹ and their typical stellar wind velocity is ≈1850 km s⁻¹ (see Table 2 of Martins et al. 2008); and QC, modeled as an injection rate of A₀ = 1.1 × 10³⁶ erg s⁻¹, considering that the total mass-loss rate of its WR population is ≈8.1 × 10⁻⁴ M_⊙ yr⁻¹ and their typical stellar wind velocity is ≈670 km s⁻¹ (see Table 2 of Rockefeller et al. 2005b). For each star cluster, we considered a t_age = 10⁶ yr, since the WR winds remain approximately constant over that timescale (Calderón et al. 2020), the actual age of the young stellar populations in each cluster is of that order of magnitude (Genzel et al. 2010), and an accurate assumption of t_age makes no difference in this timescale within the CMZ. The source locations are shown in Fig. 1, where the line of sight positions of AC and QC are taken from Hosek et al. (2022).

For impulsive sources, which correspond to supernova explosions, we defined the injection rate as Q_imp = B₀(E_p/E₀)^−Γδ(t − t_age), where B₀ is the injection energy and δ(t) is a Dirac delta function (with units of s⁻¹). Again, we assumed that a canonical 1% of the supernova explosion energy goes into acceleration of CRs above 10 TeV (Aharonian et al. 2019) and (E_p/E₀)^−Γ = 1 for mono-energetic particles with all CR energy above 10 TeV. The proposed impulsive CR source with known kinetic data is Sgr A East, which was modeled as: injection energy of B₀ = 1.5 × 10⁴⁹ erg, considering an explosion energy of ~1.5 × 10⁵¹ erg (Rockefeller et al. 2005a; Fryer et al. 2006). For Sgr A East, t_age = 1700 yr (Rockefeller et al. 2005a; Fryer et al. 2006), and its location is shown in Fig. 1, where the line of sight coordinate is from Rockefeller et al. (2005a).

Finally, motivated by the mismatch between our models and the observations at large radii that we obtained in Paper I (see also Sect. 3.5 below), we considered an additional not well-constrained impulsive source associated with Sgr D. Radio observations have identified three structures close to Sgr D (G1.1−0.1, G1.0−0.2 and G0.9+0.1), which are consistent with supernovae (Heywood et al. 2022). Some studies indicate that these structures are independent of the CMZ, appearing near it only in projection (Kauffmann et al. 2017). However, other studies indicate that according to their distances along the line of sight, these structures could be GC objects (Sofue 1995). Specifically, Sofue (1995) located G1.1−0.1, G1.0−0.2 and G0.9+0.1 at 8.9 kpc, 8.0 kpc and 9.6 kpc, respectively, with an error of ≈30%. Given that uncertainty, all of them could be within or outside the CMZ. We therefore arbitrarily assumed that only the middle structure, G1.1−0.1, is inside the CMZ and located on its far side, as suggested by Sofue (2022). We modeled G1.1−0.1 (SN Sgr D) as: assumed injection energy of B₀ = 5 × 10⁴⁸ erg, considering an explosion energy of 5 × 10⁵⁰ erg, which is consistent with the less energetic Galactic supernovae (Acero et al. 2016). Moreover, for SN Sgr D we assigned t_age = 10⁴ yr, according to its diameter (Sofue 2022) and indicate its location in Fig. 1.

Fig. 1 CMZ gas distribution in our model, schematic representation of CR dynamics, and CR source location. Top panel: CMZ particle column density spatial distribution computed from CO (J = 3–2). The CMZ was covered between −1.07° < l < 1.73° and −0.42° < b < 0.40°, within a velocity range |υ| < 220 km s−1 (Oka et al. 2012). The color scale represents the column number density of H2, and a low diffusion of CRs is considered within this region. The cyan outer area corresponds to regions where no gas column density was measured, for which a high diffusion of CRs is computed. The green area denotes the zone with polar advection. Bottom panel: CMZ line of sight distribution boundaries for both the disk and ring configurations, according to Paper I. The red ellipse is the external boundary for both scenarios, and the blue ellipse is the internal boundary for the ring scenario. A low diffusion of CRs is defined inside the red ellipse for the disk and between the red and blue ellipses for the ring. High diffusion of CRs is computed over the cyan outer area and inside the red ellipse for the ring scenario. RE1 and RE2 are the semi-major and semi-minor axis for the external ellipse, and RI1 and RI2 are the semi-major and semi-minor axis for the internal ellipse. Gray zones demarcate the area projected along the line of sight analyzed in the HESS work, where RP1 and RP2 are the limits for the projected radii from Sgr A*.

2.1.2 Differential diffusion

As the CMZ is a very high density region, a low diffusion coefficient is appropriate. On the other hand, outside of the CMZ conditions should be similar to those in the broader Galactic disk. Therefore, we considered a low isotropic CR diffusion inside the CMZ defined as D_in(E_p) ~ 10²⁶(E_p/10¹⁰ eV)^δ cm² s⁻¹ (Ormes et al. 1988; Abeysekara et al. 2017; Aharonian et al. 2019), and a high isotropic CR diffusion outside the CMZ described by D_out(E_p) ~ 10²⁸(E_p/10¹⁰ eV)^δ cm² s⁻¹ (Aharonian 2004; Strong et al. 2007). For both regimes, we considered δ = 0.5 for a Kraichnan spectrum for magnetic turbulence (Aharonian 2004).

We defined the CMZ borders in Galactic longitude and latitude by the limits of the CO (J = 3−2) observation detailed in Oka et al. (2012), and the line of sight edges by the distributions proposed in Paper I, which are shown in Fig. 1. The upper panel shows the area without gas observation in cyan, which is considered to have the diffusion coefficient of the Galactic disk, while where the observed column density is shown we used the diffusion coefficient for high densities. The lower panel shows the boundaries of the CMZ as a ring or disk according to Paper I, where the red ellipse is the outer border of the ring and disk, and the blue ellipse is the inner border of the ring. Within these borders we considered a low diffusion, and outside a high diffusion.

2.1.3 Polar advection

The observation of the FBs suggest polar advective transport. According to Bordoloi et al. (2017), the outflow velocity of the FBs is ≈ 1000–1300 km s⁻¹, their age is ~6−9 Myr, and, at small Galactic latitudes (b < 15°), the FBs have a beam angle of ≈ 120° from Sgr A*. To quantify the relative importance of advection and diffusion transport, it is useful to define a dimensionless number q = VR_trans/D, where R_trans is the transport radius. The CR dynamics are dominated by advection if q ≫ 1 and dominated by diffusion if q ≪ 1 (Owens & Jokipii 1977). In our case q ≈ 0.8, considering a CMZ thickness of ≈75 pc [≈5°, Oka et al. 2012, for a GC distance of 8.5 kpc] and the inner diffusion coefficient for 1 PeV particles. As q is somewhat smaller than unity, advection can contribute considerably but not completely dominate the CR dynamics in the area of the FBs. In our model, we defined the polar advection as a ballistic motion with a velocity of 1000 km s⁻¹ toward the north Galactic pole for all particles at an angle of 120° above Sgr A*, and to the south Galactic pole for all particles at an angle of 120° below Sgr A*, as is shown in green on the top panel of Fig. 1.

2.1.4 Energy loss

π⁰ decay is by many orders of magnitude the dominant process of energy loss for CRs protons (Aharonian et al. 2009; Longair 2011). This process occurs when high energy protons, or heavier nuclei, collide with protons in the ambient gas. These inelastic collisions produce secondary particles, such as π⁰ (e.g., p + p → p + p + π⁰), which then decay into gamma rays (π⁰ → 2γ; Aharonian 2004; Longair 2011; Gaisser et al. 2016). However, the particle’s average time within the high density volume, that is the CMZ, is much shorter than the characteristic cooling time. As such, the proton cooling is neglected in the transport calculation (third term on the right of Eq. (1)), and only used to generate gamma-ray synthetic maps (Sect. 2.3).

2.1.5 Numerical simulation

We simulated mono-energetic particles of 1 PeV using the Monte Carlo method in a 3D domain to obtain the modeled CR energy density (w_CR), where w_CR ∝ ψ. We executed 10 CR simulations of 10⁶ test particles, for each source and CMZ shape (disk or ring). We discretized the results on a 3D grid centered on Sgr A*, where w_CR was computed in each bin. The grid covers Galactic longitudes between −1.07° < l < 1.73°, Galactic latitudes within −0.42° < b < 0.40°, and ±132 pc along the line of sight. We assumed a GC distance of 8.5 kpc, with cubic bins of size 4 pc × 4 pc × 4 pc, resulting in an angular resolution of ≈0.03°. Finally, for each source and CMZ model, all 10 simulations were averaged for every bin.

2.2 3D CMZ distribution

We considered the same 3D gas distributions as we did in Paper I. For more details, see Sect. 2.2 of that paper, which we summarize here.

We developed two models to distribute the CMZ gas along the line of sight based on Kruijssen et al. (2015) and Launhardt et al. (2002). We considered that the CMZ could be an elliptical ring or an elliptical disk to simulate a CMZ with (Sofue 1995; Kruijssen et al. 2015) and without (Sawada et al. 2004; Yan et al. 2017) an inner cavity. Notice that the latter option includes also models that do consider a cavity (Ridley et al. 2017; Armillotta et al. 2019), but which is small or shallow enough so it can be neglected for our purposes. Within the disk or ring, we distributed uniformly the column density of molecular gas along the line of sight domain, to obtain an approximation of the 3D CMZ distribution. The total CMZ column density was measured from CO (J = 3−2) lines observed by the Atacama Submillimeter Telescope Experiment (ASTE; Kohno et al. 2004) and published by Oka et al. (2012)¹. We adopted a mass conversion factor CO-to-H₂ of X_CO = 1.8 × 10²⁰ cm⁻² K⁻¹ Km⁻¹ s (Bolatto et al. 2013) and a ratio of CO(J = 3−2)/CO(J = 1−0) = 0.7 (Oka et al. 2012) at the GC.

The top panel of Fig. 1 shows the particle column density of the CMZ, while the bottom panel shows its adopted distribution along the line of sight, where the red ellipse is the common external border, and the blue ellipse is the internal border for the ring model. Finally, we computed the gas particle density (n_H) using the same 3D grid described in Sect. 2.1.5, considering the gas to be static.

2.3 Gamma-ray synthetic maps

As described in Sect. 2.1.4, gamma rays are produced from proton-proton interactions. We computed the gamma-ray luminosity per unit volume (L_γ/V_bin) in the same grid of Sect. 2.1.5, using the values of w_CR and n_H obtained in each grid bin as (Fatuzzo et al. 2006; H.E.S.S. Collaboration 2016) (2)

where κ_π is the fraction of kinetic energy of high-energy protons transferred to π⁰ production, σp–p is the cross section for proton-proton interaction and η_N is the gamma-ray contribution from heavier nuclei in CRs and ambient gas. For protons with energies in the GeV-TeV range, κ_π ≈ 0.18 (Fatuzzo et al. 2006), for mono-energetic particles of 1 PeV, σp–p ≈ 53 mb (Aharonian 2004), and η_N ≈ 1.5 (H.E.S.S. Collaboration 2016). Next, we integrated L_γ/V_bin along each line of sight to obtain gamma-ray synthetic maps, where the gamma-ray luminosity per bin cross section was converted a gamma-ray surface flux computed at the Earth reference frame.

Our model neglects the gamma-ray foreground and background of the CMZ (Abramowski et al. 2014) emitted by the isotropic Galactic CRs sea (Blasi 2013). However, this large-scale emission is considered to be relatively weak and should have a negligible contribution in the central region (H.E.S.S. Collaboration 2018a), such that it will not change our main results.

Fig. 2 Gamma-ray synthetic maps computed as observed from the Earth. Left panels: synthetic maps for the CMZ as a disk and different CR sources (indicated on each panel). Top four right panels: synthetic maps for the CMZ as a ring and different CR sources (indicated on each panel). All maps have been smoothed with a 0.03° Gaussian function to adopt the best CTA beamwidth. Bottom right panel: regions where we contrasted our models with the H.E.S.S. Collaboration (2016) results.

3 Results

First, we analyzed the results considering only the sources whose identity and energy is well constrained (i.e., NSC, Sgr A East, AC, and QC). For those sources, we contrasted the results with both HESS observations and our Paper I models. Also, we studied the future constraints to be obtained by CTA observations. Finally, we added the proposed injection from SN Sgr D to improve the model at larger Galactic longitude.

3.1 Model results and current observation

Figure 2 shows the results of our models considering the four sources with CR injection derived from observations. Left panels show the results of the CMZ-as-a-disk model, considering NSC (top panel), NSC and Sgr A East (second panel), NSC, AC and QC (third panel), and NSC, Sgr A East, AC and QC (fourth panel) as the CR sources. The top four panels on the right show the result with the same sources but considering the CMZ as a ring. All maps have been adjusted to CTA resolution, that is have been smoothed with a 0.03° Gaussian function to adopt the best CTA beamwidth.

Comparing all maps, they all show a central overdensity of gamma-ray flux and lower emission at high longitudes. However, models where the CMZ is considered as a disk and Sgr A East is included as a CR source show more emission at larger longitudes. In the case of continuous sources (NS, AC and QC), the particles injected in the past can escape the CMZ due to the external high diffusion, and the gamma-ray emission close to the sources is high because “young” particles have only traveled in the low-diffusion region. In the case of the impulsive source, two different effects are experimented by the particles on either gas configuration. Given the low t_age for Sgr A East, the particles are concentrated around the source due to the low diffusion in the disk shape. On the other hand, when the CMZ is a ring, the particles move away quickly, and then get confined when entering the CMZ, as there is not enough time for them to leave again the low-diffusion ring, creating an extended over density.

In order to contrast all models with the HESS observations, we integrated the gamma-ray surface flux over the areas shown in the bottom right panel of Fig. 2. Thus, we obtained gamma-ray luminosities (L_γ) for the same regions selected in H.E.S.S. Collaboration (2016). These areas correspond to three annular sectors centered in Sgr A* (l = −0.056°, b = −0.04588°), with inner and outer radii of 0.1°/0.15°, 0.15°/0.2°, 0.2°/0.3°, and excluding the region between the angles +10^° and −56^° from the positive Galactic longitude axis, plus seven circular regions with 0.1° of radius centered in b = −0.04588° and l = −0.656°, −0.456°, 0.344°, 0.544°, 359.344°, 0.744°, 0.944°, and 1.144°. The comparison of Lγ for all models with the HESS data points is show in Fig. 3, where luminosities are radially projected from Sgr A*, therefore two data points are plotted at ≈60 pc and ≈90 pc due to the observation areas selected by HESS. Only the model with the CMZ as a ring and with injection from all sources (NSC, Sgr A East, AC and QC) reproduces satisfactorily the data for projected radii R ≤ 150 pc. This match is relevant because the CMZ morphology, CR diffusion coefficients, CR polar advection and CR sources have been constrained from observations and previous literature. In conclusion, more realistic CR dynamics suggest that the CMZ presents an inner cavity and that the GC PeVatron can be reproduced by the four modeled sources.

In Paper I, we obtained two scenarios consistent with the H.E.S.S. Collaboration (2016) results but considering a simplified CR transport (i.e., one isotropic diffusion coefficient, no advection, and mono-energetic particles of 1 PeV), which is generally accepted but does not represent properly the GC environment. Here, in the left panel of Fig. 4, we contrasted the gamma-ray emission from the best-fitting model of this work with the previous results. Additionally, we computed the average energy density along the line of sight (w_CR) within the same areas used to calculate L_γ. For that, we considered only the gas column density in the CMZ shown in the top panel of Fig. 1 (i.e., ignoring the CMZ line of sight distribution, as done typically in the literature), n_H = N_H/l_los, and V = A l_los, where A is the observed gamma-ray area and l_los is the length of the domain on the line of sight, so Eq. (2) is rewritten as: (3)

In the right panel of Fig. 4 we contrasted the w_CR for the CMZ ring and disk models without differential diffusion and polar advection developed in Paper I, and the model with the CMZ as a ring, considering differential diffusion, polar advection, and CR injection from all sources (NSC, Sgr A East, AC and QC). In both panels, we observed that in general all models are similar and fit the observations, with only some differences at R ≈ 60 pc, where the new model deviates a bit more. At a projected radius of R ≈ 175 pc, all models predict a luminosity of a factor of several below the data (see Sect. 3.5 below). All in all, considering that we already have these three satisfactory models of Fig. 4, we maintain that fine tuning them is not currently needed. Rather, we investigated whether future CTA observations will allow us to choose which of the models is the best match.

Fig. 3 CMZ gamma-ray luminosity profiles, computed using the model with realistic CR dynamics introduced in this paper. The luminosities are calculated by integrating the flux in the regions shown in the bottom right panel of Fig. 2. Left panel: gamma-ray profiles considering the CMZ as a disk, for all models shown in the left panels of Fig. 2. Right panel: gamma-ray profiles considering the CMZ as a ring, for the models shown in the top four right panels of Fig. 2. Black crosses denote the observed gamma-ray luminosity computed by H.E.S.S. Collaboration (2016) in the CMZ. Yellow squares, blue circles, gray diamonds and red triangles, represent the gamma-ray luminosity considering the NSC, NSC+Sgr A East, NSC+AC+QC, and NSC+Sgr A East+AC+QC, as CR sources, respectively.

Fig. 4 Gamma-ray luminosity and average CR energy density for our satisfactory models. Left panel: gamma-ray luminosity profiles from Sgr A* along the Galactic disk, taken from Fig. 3 and from Paper I. Right panel: corresponding average CR energy densities. Black crosses denote the observed gamma-ray luminosity (left panel) and the indirect observations of CRs (right panel) computed by H.E.S.S. Collaboration (2016). Yellow squares and cyan diamonds correspond to the CMZ ring and disk models without differential diffusion and polar advection developed in Paper I. On the disk model, the only CR source is the NSC, while on the ring model the CR sources are the NSC and Sgr A East. Red triangles correspond to the new model with the CMZ as a ring, considering differential diffusion, polar advection, and CR injection from all sources (NSC, Sgr A East, AC and QC).

3.2 Future observation by CTA

CTA will soon be the most powerful telescope for very-high-energy gamma-ray astronomy, covering the energy range from 20 GeV to 300 TeV. CTA will improve the sensitivity by one order of magnitude at 1 TeV compared to current observatories, with an angular resolution of ≈0.03° and an energy resolution of ≈8%. Very relevant for our study, one of its key science projects will spend 525 h of exposure on the GC (CTA Consortium 2019). In this context, we compared the features of the three models that match HESS observations (see Fig. 4) using the expected CTA sensitivity and angular resolution.

3.2.1 Energy distribution of synthetic maps

As explained in Sect. 2, to calculate the CR dynamics we concentrated all the energy of CRs above 10 TeV into mono-energetic particles of 1 PeV, so our computed gamma-ray luminosity represents the total emission from high energy particles over 10 TeV. To analyze the synthetic maps as real gamma-ray observations, we needed to distribute this luminosity in the energy domain, specifically between 1 TeV and 100 TeV, because E_γ ≈ 0.1E_p where E_γ is the gamma-ray energy (Aharonian et al. 2009; Longair 2011). H.E.S.S. Collaboration (2016) observed that the CMZ energy spectrum follows a power-law without any cutoff, which was defined as dN_CMZ/dEdAdt = 1.92 × 10⁻¹²(E_γ/TeV)^−2.32 TeV⁻¹ cm⁻² s⁻¹, where dN_CMZ/dEdAdt is the number of gamma rays emitted from the CMZ per unit of energy, area and time, as a function of E_γ. Then, we defined the total gamma-ray flux (F_γ) as (H.E.S.S. Collaboration 2018b) (4)

where Φ₀ is the flux normalization. Finally, we discretized all F_γ of each pixel over five bins per energy decade and adjusted Φ₀ for every F_γ.

3.2.2 Gamma-ray absorption

Gamma-ray emissions can be absorbed during their propagation in the Galaxy by pair production. Specifically, high-energy photons above ≈10 TeV and ≈200 TeV interact with the Galactic interstellar radiation field (ISRF) and the cosmic microwave background (CMB), respectively, producing an electron-positron pair (γ + γ → e⁻ + e⁺; Moskalenko et al. 2006; Popescu et al. 2017). Therefore, we attenuated F_γ according to the gamma-ray transmittance (i.e., the fraction of photons reaching the Earth) as a function of E_γ for gamma-ray sources located at the GC computed by Moskalenko et al. (2006), where transmittance decreases to ≈0.8 for gamma rays between 50 and 100 TeV.

3.2.3 Data analysis of synthetic maps

According to the methodology described in Mohanty et al. (1998), we converted F_γ to photon counts following (5)

where N_ON is the number of gamma rays observed by CTA, ∆t is the observation time, A₀ is the telescope effective area, and ϵ_γ is the data taking efficiency. For the CTA telescope, A₀ = 5 × 10¹⁰ cm², ϵ_γ = 0.7, and we considered the typical observation time to be 50 h (CTA Consortium 2019). We computed N_on for each modeled pixel of angular resolution 0.027°. To estimate the background, we used the GC gamma-ray observation of the H.E.S.S. Galactic plane survey (H.E.S.S. Collaboration 2018b) and the field-of-view background method, which is appropriate for extended sources (Berge et al. 2007). We selected an area centered on Sgr A* of 4.5°, which is equal to the minimum CTA field of view (CTA Consortium 2019), and excluded known gamma-ray sources under the assumption that they will be easily isolated by CTA. From the entire resulting area, we obtained a background count of N_OFF = 259 006 in a total solid angle of Ω_OFF = 0.00273 sr for the same 50 h of observation. Finally, we estimated the number of excess events over the background (N_γ) for each pixel as (H.E.S.S. Collaboration 2018b) (6)

where Ω_ON is the pixel solid angle and α = Ω_ON/Ω_OFF = 8 × 10⁻⁵. Then we computed the gamma-ray statistical significance for each event (N_ON, N_OFF) according to Li & Ma (1983), where the event standard deviation (σ) is defined by (7)

To estimate the influence of the systematic uncertainty and statistical error, we computed the signal-to-background ratio (SBGR) and the signal-to-noise ratio (S/N) as (8) (9)

where N_bg = αN_OFF, N_CR is the uncertainty in counts associated to the statistical error of the CR model, and N_CTA is the uncertainty in counts associated to the telescope systematic error. Finally, we considered that gamma-ray emission is actually detected when N_γ ≥ 10, σ ≥ 5, SBGR ≥ 0.05 and S/N ≥ 5 (Actis et al. 2011), considering an uncertainty in the background of 1% (Funk et al. 2013).

Figure 5 shows the synthetic maps of gamma-ray surface flux, considering 50 h of observation with CTA and E_γ ≥ 1 TeV, where contours demarcate the detection criteria. The clearest differences between the maps are with regard to the extension and the CR sources. In order to compare in detail our models, Fig. 6 shows selected longitudinal profiles with prominent features, where F_γ was obtained between −0.797° < l < 1.337° and for b = 0.094°, 0.013°, −0.067° and −0.229°. A Al and Ab of 0.03°, comparable to the CTA beamwidth, was adopted for all cases. In addition, we did not consider the central 0.2° around Sgr A* because HESS observes here a point-like source with a different energy spectrum than the CMZ (H.E.S.S. Collaboration 2016). Moreover, the amount of gas in this region, including the circumnuclear disk (Genzel 1988), is still under debate (Chen et al. 2016). In all panels, the differential diffusion effect is observed at l ≈ 1.2° and l ≈ −0.6°, where F_γ decrease due to the external high diffusion. The effect is similar in both directions but at a different relative radius with respect to Sgr A* (l ≈ −0.05°), due to the asymmetric distribution of the CMZ gas. In the bottom-right panel, the polar advection effect becomes noticeable, where F_γ decreases at Sgr A* for the model that includes advection and a roughly flat profile is observed for the other models. This is due to advection allowing CRs to escape faster than in other areas of the CMZ. When the CR dynamics is dominated by diffusion or advection it will decrease ∝ 1/r or ∝ 1/r² respectively (H.E.S.S. Collaboration 2016), such that the effect is more prominent farther away from Sgr A* and NSC. In the top-right and bottom-left panels, the CR injection from AC (l ≈ 0.12° and b ≈ 0.02°) and QC (l ≈ 0.16° and b ≈ −0.06°) are resolved as F_γ peaks. On the other hand, the models without differential diffusion and polar advection show a central small bump at different longitudes (l ≈ −0.05° for the disk model and l ≈ −0.08° for the ring), which traces the central w_CR peak. This discrepancy is due to their respective gas distributions along the line of sight.

3.3 CMZ gamma rays over 50 TeV

The CR spectrum at 1 PeV energies is still not well-constrained observationally. Although the HESS analysis is consistent with protons of ~ PeV, observation of gamma rays close to 100 TeV would be unequivocal evidence for a GC PeVatron. Since such particle acceleration is assumed in our model, we analyzed if CTA would detect gamma rays over 50 TeV. At gamma-ray energies over ≈10 TeV, the photon detection method is dominated by N_γ ≥ 10, where the telescope array needs to cover ≥km² (Funk et al. 2013) or a very long observation time is needed (Eq. (5)). In this context, H.E.S.S. Collaboration (2016) did not observe any gamma rays over ≈50 TeV with ≈200 h of exposure time, while our models with 50 h of observation with CTA are dominated by photons between 1 TeV and 10 TeV. To address this issue, we repeated the data analysis (see Sect. 3.2.3), but only considering gamma rays above 50 TeV, 525 h of observation with CTA, and a background that follows a power-law with a spectral index of −2.57 (Bernlöhr et al. 2013). Figure 7 shows the gamma-ray flux as observed by CTA for photons over 50 TeV. In all models, the CMZ is partially filled by gamma rays and extended emission is observed in the central area, where gamma rays between 50 and 100 TeV are detected even considering the gamma-ray absorption by pair production. CTA will be able to trace CRs of ≈ PeV energy, and their detection will cover the area of all proposed sources, with AC and QC being individually resolved.

Therefore, CTA will directly confirm whether clusters of young massive stars constitute PeVatrons. Such a detection would provide relevant information about the Galactic origin and acceleration of CRs up to 1 PeV without all the complications of the NSC. In the latter a point-like gamma-ray source is observed (H.E.S.S. Collaboration 2016), but its origin remains ambiguous due to other potential options: Sgr A* could inject CRs (Aharonian & Neronov 2005); gamma-ray emission from a dark matter peak is theoretically expected (Belikov et al. 2012); an unresolved population of millisecond pulsars could accelerate particles (Guépin et al. 2018; Bower et al. 2018); there is a pulsar wind nebula on the line of sight (Wang et al. 2006); stellar winds could be interacting further with an outflow coming from the supermassive black hole vicinity (Cuadra et al. 2015; Russell et al. 2017); and we suggested an additional CR injection from Sgr A East (Paper I).

Fig. 5 Synthetic maps of gamma-ray flux, considering 50 h of observation with CTA and Eγ ≥ 1 TeV. Top and middle panels: gamma-ray flux of CMZ models without differential diffusion and polar advection, reported in Paper I. Bottom panel: gamma-ray flux of CMZ model with differential diffusion and polar advection. CR sources and the CMZ morphology are indicated on each panel. White, yellow, cyan and green contours are lower limits to Nγ, σ, SBGR and S/N, respectively. All maps have been smoothed with a 0.03° Gaussian function to adopt the CTA beamwidth.

3.4 Gamma rays from Sgr A*

As mentioned previously, H.E.S.S. Collaboration (2016) reported a point-like source consistent with emission from the location of Sgr A*. The main difference between Sgr A* observation and the gamma rays from the CMZ is the energy cutoff at 10 TeV observed from this point-like source. We integrated the gamma-ray surface flux over the same area that was observed by HESS (i.e., a circular region of radius 0.1° centered in Sgr A*), and our model predicts a luminosity four times lower than the luminosity obtained from integrating the observed spectrum between 1 and 100 TeV. Also, if we consider that our prediction follows the CMZ observed spectrum, the flux magnitude of the cutoff is above our modeled luminosity.

We explored including a new CR source at the location of Sgr A* so that the model better matches the point-like source observation, however, any additional CR population modified strongly the profiles shown in Fig. 4 when interacting with the CMZ gas. Therefore, our NSC+Sgr A East+AC+QC model (with the CMZ as a ring) requires an additional gamma-ray source with an energy cutoff at 10 TeV, an extension smaller than 0.2°, and located at the Sgr A* coordinates.

3.5 CR injection from SN Sgr D

In Paper I, we attributed the discrepancy at R ≈ 175 pc to more realistic CR dynamics and/or an additional CR source. As is evident from Fig. 4, more realistic dynamics did not solve this problem. Therefore we analyzed the possible injection of CRs from SN Sgr D to fit this detection. In Fig. 8, we added to the NSC+Sgr A East+AC+QC model (CMZ as a ring) the SN Sgr D simulation. The sum reproduces the observed gamma-ray emission satisfactorily and improves our model, although assuming arbitrary yet reasonable values for the energy of the explosion and the line-of-sight distance for that source.

To estimate the impact of adding SN Sgr D to our best model with differential diffusion and polar advection, we recomputed the synthetic maps and profiles observed by CTA, finding that SN Sgr D does not modify our previous analysis. The maps and profiles considering SN Sgr D are shown in Figs. A.1 and A.2.

Fig. 6 Surface gamma-ray flux profiles for successful models. Profiles are obtained between 0.797° < l < 1.337° and b = 0.094° (top-left panel), 0.013° (top-right panel), −0.067° (bottom-left panel) and −0.229° (bottom-right panel). The profile is binned adopting the CTA angular resolution, 0.03°. Lines show gamma-ray flux considering 50h of observation with CTA and Eγ ≥ 1 TeV. Yellow and cyan profiles represent the models proposed in Paper I with the CMZ as a disk and a ring. The red profiles are from the model considering differential diffusion, polar advection, and CR injection from NSC, Sgr A East, AC, and QC (all sources). The central 0.2° around Sgr A* was excluded.

4 Discussion

In Paper I, we found two models for the CR dynamics in the GC that were consistent with HESS gamma-ray observations. Those models however were too simplistic, and neglected several expected physical effects, such as differential diffusion and polar advection, and CR sources, such as the Arches and Quintuplet stellar clusters. These ingredients were included in the models presented here, and we obtained a new model that also recreates the HESS data. This match is remarkable, as it is based on a priori reasonable assumptions, without fine-tuning any parameters. For all physical processes and sources we used fiducial values based on observations and/or the current literature. However, other combinations of sources and processes could also fit the data, and of course there are uncertainties in the fiducial values we used (e.g., the rate of kinetic energy going into acceleration of CRs above 10 TeV). Therefore, we do not claim that our model corresponds to reality, and instead simply argue that upcoming CTA observations will be able to constrain the models in unprecedented detail.

According to our observational analysis, the different HESS-matching models will show different features when imaged by CTA. The CR realistic transport model proposed in this work will present the most different signature, where CTA will be able to constrain the differential diffusion effect close to the CMZ edge and verify the polar advection to high or low latitudes at the longitude of Sgr A*. Moreover, thanks to its fine angular resolution of ≈0.03° (CTA Consortium 2019), it will be possible to detect the gamma-ray peaks associated to AC, QC, and a composite source of both NSC and Sgr A East, independently, as well as any other potential CR acceleration sites so far not considered. Additionally, using these future observations, we could fine tune our model to determine observationally the CMZ diffusion coefficient, polar advection magnitude, the source injection rates, and the gas distribution along the line of sight.

Regarding gamma-ray emission over 50 TeV, we verified the main results for more strict transmittances (e.g., ≈0.75 for gamma rays between 50 and 100 TeV; Popescu et al. 2017;Porter et al. 2018). The conclusions are the same, CTA will directly confirm PeVatrons in the GC, and the observed morphology will only turn out to be less extended due to more intense ISRF. Additionally, the spectrum that CTA will observe should be corrected due to gamma-ray absorption to obtain the intrinsic CR spectrum to ~1 PeV energies.

As shown above, our model remains consistent with the current observations of Sgr A*, and requires an additional pointlike gamma-ray source at that location. This emission has been argued to arise from the vicinity of the supermassive black hole (Aharonian & Neronov 2005) or the ambient gas in the inner few pc (Chernyakova et al. 2011), within our modeled cavity.

While almost all our sources are already well-constrained from observations, SN Sgr D remains poorly determined. Its currently assumed properties are consistent with its diameter (Sofue 2022) and with the energetics of the weaker supernovae in our galaxy (Acero et al. 2016). According to our simulations, a verification of this source by CTA may be very difficult due to its low CR injection rate. Additionally, a more accurate determination of its line-of-sight location is necessary.

The model presents some limitations and uncertainties, which could affect our result. Most importantly, we followed H.E.S.S. Collaboration (2016) and used the standard Galactic CO-to-H₂ mass conversion factor for the CMZ. However, other studies estimate it to be lower by up to an order of magnitude (Sodroski et al. 1995; Dahmen et al. 1998; Sofue 2022). In Fig. 9, we recomputed the NSC+Sgr A East+AC+QC model (CMZ as a ring) setting the CO-to-H₂ mass conversion factor to X_CO = 1.8 × 10¹⁹ cm⁻² K⁻¹ Km⁻¹ s, that is one order of magnitude lower than the Galactic value, considering the same CR dynamics and injection as our satisfactory model from Sect. 2 (brown squares), a 10 times higher CR injection rate (orange circles), a ten times lower diffusion outside the CMZ (light blue x-marks), a ten times lower diffusion inside the CMZ (green triangles), and a ten times lower diffusion both inside and outside the CMZ (purple diamonds). Analyzing these models, it is evident that the CR energy density inferred by HESS is in strong tension with the lower value of the gas conversion factor and therefore strong additional sources would be required to match the gamma-ray observations. For instance, the model with the same sources but ten times higher injection rates, produces the same luminosity as in our fiducial model. We neglected additional potential CR sources, such as the Sgr B1 and Sgr B2 star formation complexes (Ginsburg & Kruijssen 2018; Henshaw et al. 2023) or other supernova remnants that are visible in radio and X-ray maps (e.g., Heywood et al. 2022), because our aim in this study was to compute a fiducial model supported by known CR sources with well-constrained kinetic energies and low uncertainty on whether they are within the CMZ (Crocker et al. 2007; Protheroe et al. 2008; Acero et al. 2016; Kauffmann et al. 2017). In future works, we will explore the impact of a more exhaustive characterization of these potential sources. On the other hand, considering only the known sources, the CMZ mass conversion factor would need to be close to the standard Galactic one.

Other uncertainties include that the positions of AC and QC along the line of sight are not well constrained. However, repeating the simulations with the locations proposed by Kruijssen et al. (2015) on the near side of the CMZ, we find that the results do not vary substantially, as the clusters are still embedded in high-density gas. Additionally, CR diffusion and advection areas were delimited from molecular emission lines and gamma-ray observations, respectively, which can only approximate the actual CR dynamics. Moreover, the transition between zones could be smoother.

Finally, our results are consistent with recent theoretical models of the CMZ and Westerlund 1 (a cluster of young massive stars with a gamma-ray emission similar to the CMZ; Aharonian et al. 2019). For the CMZ, a composite population of CRs within the GC was constrained to reproduce current observations (Becker Tjus et al. 2022), and for Westerlund 1, CR acceleration at the wind termination shock within massive stars clusters and a low diffusion in high density zones were derived to be consistent with gamma-ray observations (Bhadra et al. 2022).

Fig. 7 Same as Fig. 5, but considering 525 h of observation with CTA and Eγ ≥ 50 TeV. The additional pink dotted contours correspond to the maps perimeter without considering gamma-ray absorption, while red dashed contours mark the perimeter of the synthetic maps of Fig. 5.

Fig. 8 Same as Fig. 4, but now plotting the gamma-ray luminosity and the indirect observations of CRs for the CMZ ring, considering differential diffusion, polar advection, and CR injection from SN Sgr D (green squares). Magenta circles are the sum of both components (all sources+SN Sgr D).

Fig. 9 Same as the model of all sources and CMZ as a ring in Fig. 4, but setting the CO-to-H2 mass conversion factor to one order of magnitude lower than the Galactic value. We computed models considering the same previous dynamics and CR injection (brown squares), ten times more CR injection (orange circles), ten times less diffusion coefficient outside the CMZ (light blue x-marks), ten times less diffusion coefficient inside the CMZ (green triangles), and ten times less diffusion coefficient inside and outside the CMZ (purple diamonds).

5 Conclusion

We developed a Monte Carlo CR diffusion model for the GC CRs considering differential diffusion, polar advection and several sources. We studied the gamma rays emitted via protonproton collision with the CMZ gas, creating synthetic gamma-ray maps. We considered two different CMZ 3D models, with either a ring or a disk shape, that simulate a significant inner cavity or not. Our simulation results are consistent with the former, implying that the CMZ has an inner cavity and that the CR injection is from the combination of the NSC, Sgr A East, AC and QC. The resulting gamma-ray luminosity was found using fiducial parameters and normalizations derived from observations and the current literature, including the standard Galactic CO-to-H₂ conversion factor. If a CO-to-H₂ conversion factor one order of magnitude lower is considered, the same sources would need injection rates that are ten times higher to match the observed gamma-ray luminosity. Additionally, an impulsive emission from SN Sgr D can explain the relatively high gamma-ray emission detected by HESS at a large projected radius from Sgr A*. Our conclusion on the shape of the CMZ is consistent with several kinematic and dynamical gas models that require an inner cavity, but stands in contrast to models based on observations of molecular lines and some dynamical gas models, which suggest a continuous structure. Finally, our predictions do not overproduce the current gamma-ray spectrum observed from Sgr A*, for which additional physical mechanisms have been proposed.

To verify whether this model, or the models without differential diffusion and polar advection propounded in Paper I, explain better the gamma-ray emission, we studied the signatures that CTA could detect in the future. Considering 50 h of observation, CTA will distinguish those models, constraining the effects of differential diffusion, polar advection, CR sources, and the CMZ morphology. Moreover, with the planned total 525 h of observation, CTA should detect gamma rays close to 100 TeV, which will be unequivocal evidence of a GC PeVatron, and even confirm whether clusters of young massive star without a super-massive black hole can accelerate particles up to 1 PeV.

Acknowledgements

We thank the anonymous referee, and AS’s thesis reviewers, Mario Riquelme and Rolando Dünner, for constructive comments that helped us improve the paper. This project was partially funded by the Max Planck Society through a “Partner Group” grant. A.S. acknowledges the help and useful comments by Brian Reville, the hospitality of the Max Planck Institute for Nuclear Physics, where part of the work was carried out, and funding from the Deutscher Akademischer Austauschdienst (DAAD). We thank Roberto Lineros for useful discussions at an early stage of this project. A.S. and J.C. acknowledge financial support from FONDECYT Regular 1211429. F.E.B. acknowledges support from ANID-Chile BASAL CATA FB210003, FONDECYT Regular 1200495 and 1190818, and Millennium Science Initiative Program – ICN12_009. The Geryon cluster at the Centro de Astro-Ingenieria UC was extensively used for the calculations performed in this paper. BASAL CATA PFB-06, the Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL 130008 provided funding for several improvements to the Geryon cluster.

Appendix A Maps and profiles considering SN Sgr D

Fig. A.1 Same as bottom panel of Figs. 5 (left panel) and 7 (right panel), but now adding SN Sgr D.

Fig. A.2 Same as Fig. 6, but now plotting all sources (NSC, Sgr A East, AC, and QC) plus SN Sgr D.

References

All Figures

Fig. 1 CMZ gas distribution in our model, schematic representation of CR dynamics, and CR source location. Top panel: CMZ particle column density spatial distribution computed from CO (J = 3–2). The CMZ was covered between −1.07° < l < 1.73° and −0.42° < b < 0.40°, within a velocity range |υ| < 220 km s−1 (Oka et al. 2012). The color scale represents the column number density of H2, and a low diffusion of CRs is considered within this region. The cyan outer area corresponds to regions where no gas column density was measured, for which a high diffusion of CRs is computed. The green area denotes the zone with polar advection. Bottom panel: CMZ line of sight distribution boundaries for both the disk and ring configurations, according to Paper I. The red ellipse is the external boundary for both scenarios, and the blue ellipse is the internal boundary for the ring scenario. A low diffusion of CRs is defined inside the red ellipse for the disk and between the red and blue ellipses for the ring. High diffusion of CRs is computed over the cyan outer area and inside the red ellipse for the ring scenario. RE1 and RE2 are the semi-major and semi-minor axis for the external ellipse, and RI1 and RI2 are the semi-major and semi-minor axis for the internal ellipse. Gray zones demarcate the area projected along the line of sight analyzed in the HESS work, where RP1 and RP2 are the limits for the projected radii from Sgr A*.

In the text

Fig. 2 Gamma-ray synthetic maps computed as observed from the Earth. Left panels: synthetic maps for the CMZ as a disk and different CR sources (indicated on each panel). Top four right panels: synthetic maps for the CMZ as a ring and different CR sources (indicated on each panel). All maps have been smoothed with a 0.03° Gaussian function to adopt the best CTA beamwidth. Bottom right panel: regions where we contrasted our models with the H.E.S.S. Collaboration (2016) results.

In the text

Fig. 3 CMZ gamma-ray luminosity profiles, computed using the model with realistic CR dynamics introduced in this paper. The luminosities are calculated by integrating the flux in the regions shown in the bottom right panel of Fig. 2. Left panel: gamma-ray profiles considering the CMZ as a disk, for all models shown in the left panels of Fig. 2. Right panel: gamma-ray profiles considering the CMZ as a ring, for the models shown in the top four right panels of Fig. 2. Black crosses denote the observed gamma-ray luminosity computed by H.E.S.S. Collaboration (2016) in the CMZ. Yellow squares, blue circles, gray diamonds and red triangles, represent the gamma-ray luminosity considering the NSC, NSC+Sgr A East, NSC+AC+QC, and NSC+Sgr A East+AC+QC, as CR sources, respectively.

In the text

Fig. 4 Gamma-ray luminosity and average CR energy density for our satisfactory models. Left panel: gamma-ray luminosity profiles from Sgr A* along the Galactic disk, taken from Fig. 3 and from Paper I. Right panel: corresponding average CR energy densities. Black crosses denote the observed gamma-ray luminosity (left panel) and the indirect observations of CRs (right panel) computed by H.E.S.S. Collaboration (2016). Yellow squares and cyan diamonds correspond to the CMZ ring and disk models without differential diffusion and polar advection developed in Paper I. On the disk model, the only CR source is the NSC, while on the ring model the CR sources are the NSC and Sgr A East. Red triangles correspond to the new model with the CMZ as a ring, considering differential diffusion, polar advection, and CR injection from all sources (NSC, Sgr A East, AC and QC).

In the text

Fig. 5 Synthetic maps of gamma-ray flux, considering 50 h of observation with CTA and Eγ ≥ 1 TeV. Top and middle panels: gamma-ray flux of CMZ models without differential diffusion and polar advection, reported in Paper I. Bottom panel: gamma-ray flux of CMZ model with differential diffusion and polar advection. CR sources and the CMZ morphology are indicated on each panel. White, yellow, cyan and green contours are lower limits to Nγ, σ, SBGR and S/N, respectively. All maps have been smoothed with a 0.03° Gaussian function to adopt the CTA beamwidth.

In the text

Fig. 6 Surface gamma-ray flux profiles for successful models. Profiles are obtained between 0.797° < l < 1.337° and b = 0.094° (top-left panel), 0.013° (top-right panel), −0.067° (bottom-left panel) and −0.229° (bottom-right panel). The profile is binned adopting the CTA angular resolution, 0.03°. Lines show gamma-ray flux considering 50h of observation with CTA and Eγ ≥ 1 TeV. Yellow and cyan profiles represent the models proposed in Paper I with the CMZ as a disk and a ring. The red profiles are from the model considering differential diffusion, polar advection, and CR injection from NSC, Sgr A East, AC, and QC (all sources). The central 0.2° around Sgr A* was excluded.

In the text

	Fig. 7 Same as Fig. 5, but considering 525 h of observation with CTA and Eγ ≥ 50 TeV. The additional pink dotted contours correspond to the maps perimeter without considering gamma-ray absorption, while red dashed contours mark the perimeter of the synthetic maps of Fig. 5.
In the text

	Fig. 8 Same as Fig. 4, but now plotting the gamma-ray luminosity and the indirect observations of CRs for the CMZ ring, considering differential diffusion, polar advection, and CR injection from SN Sgr D (green squares). Magenta circles are the sum of both components (all sources+SN Sgr D).
In the text

Fig. 9 Same as the model of all sources and CMZ as a ring in Fig. 4, but setting the CO-to-H2 mass conversion factor to one order of magnitude lower than the Galactic value. We computed models considering the same previous dynamics and CR injection (brown squares), ten times more CR injection (orange circles), ten times less diffusion coefficient outside the CMZ (light blue x-marks), ten times less diffusion coefficient inside the CMZ (green triangles), and ten times less diffusion coefficient inside and outside the CMZ (purple diamonds).

In the text

	Fig. A.1 Same as bottom panel of Figs. 5 (left panel) and 7 (right panel), but now adding SN Sgr D.
In the text

	Fig. A.2 Same as Fig. 6, but now plotting all sources (NSC, Sgr A East, AC, and QC) plus SN Sgr D.
In the text