Implications from 3D modeling of gamma-ray signatures in the Galactic Center Region

The total density profile of the GC region is modeled to be composed of three individual components:

1. 1.

   local Molecular Clouds (MCs) - the properties of the well-known clouds in the GC region are summarized in Table 1 in Guenduez et al. (2020).

2. 2.

   the gas structure of the innermost 10 pc of the GC region is taken into account as discussed in Ferrière (2012).

3. 3.

   an analytical description of the intercloud medium is given in Ferrière et al. (2007). As compared to that original publication, the normalization factor is scaled down in Guenduez et al. (2020), as the subtraction of the sub-structures that are considered separately (see item (1) and (2)) is necessary. This procedure results in the diffuse densities for H₂ and H as $n_{0,H_{2}}=91.1$ cm^-3 and $n_{0,H}=5.3$ cm^-3. We use these numbers for the diffuse gas component in this work as well.

The total magnetic field is given by the superposition of three B-field components derived from

1. 1.

   the diffuse Inter Cloud Medium (ICM), $\@vec{B}_{\text{ICM}}$

2. 2.

   eight known non-thermal filament (NTF) structures, $\@vec{B}_{\text{NTF,i}}$ with $i=1,\ldots 8$;

3. 3.

   twelve local MC regions, $\@vec{B}_{\text{MC,j}}$ with $j=1,\ldots 12$.

The orientation of the large-scale magnetic field components in the ICM and the NTF regions is predominantly poloidal. The field is therefore modeled with an analytical, divergence-free and isotropic poloidal field model that is presented in Ferrière & Terral (2014) as model C. Normalization, opening of the field lines and lengths scales are adapted individually for the eight NTFs and for the ICM as described in Guenduez et al. (2020). An average value of $\overline{B}_{\rm ICM}\approx$10 $\mu$G determines the field strength in the ICM (Ferrière, 2009) and the typical parameters for the NTF regions can be taken from Table 1 in Guenduez et al. (2020).

The MC regions, on the other hand, are typically expected to have a horizontal field orientation. The ratio $\eta=B_{r}/B_{\phi}$ between the radial and azimuthal field is set to $\eta=0.5$ as suggested in Guenduez et al. (2020). The field structure is then derived by using Euler’s $\alpha$ and $\beta$ potentials, which deliver a divergence-free expression of the magnetic field using $\@vec{B}=\nabla\alpha\times\nabla\beta=(B_{r},\,B_{\phi},\,0)$. The magnetic field strength in the MC region can be taken from Table 1 in Guenduez et al. (2020).

The total field in the CMZ is then obtained by a superposition of the magnetic fields:

Figures 1 and 2 display the 3D configuration and strength of the total mass distribution and magnetic field in the CMZ, respectively. These configurations are used in our simulations, using the gas distribution for modeling hadronic interactions and the large-scale magnetic field model for the three-dimensional propagation of the charged particles.

The CMZ is home to more than 200 O-type stars (Kauffmann, 2017). These types of stars are rare and bright, and their temperature exceeds 3$\cdot 10^{4}$ K. For instance, 4 of the 90 brightest stars that are visible to the eye from Earth are O-type stars. They are prominent candidates for future violent supernova explosions due to their high mass. Thus, these types of stars represent energy reservoirs for CR acceleration in the CMZ. In total, three well-known star clusters in the GC exist: first the Central cluster, second the Quintuplet cluster and third the Arches cluster. Accordingly, the photon field is enhanced from FIR to UV wavelengths. Although in the CMZ the star and photon density is the highest in the Milky Way, star formation is suppressed, indicated by the lack of H₂O and methanol masers in the GC (Kauffmann, 2017). Rathborne et al. (2014) show this suppression applies to the entire CMZ.

The background photon density has been discussed to play a crucial role for gamma-gamma interactions in Porter et al. (2018). In addition, Inverse Compton scattering of high-energy electrons of the photon field can lead to the production of high-energy photons. And finally, high-energy photons undergo pair-production via their interaction with the ambient low-energy photons. In order to include these processes, this work makes use of a 3D photon distribution model as described below. Proton-gamma interactions are negligible considering the spatial extent and the energy range in the CMZ and are not included in this work.

Sophisticated three-dimensional models of the interstellar radiation field (ISRF) in the Galaxy are presented by Freudenreich (1998) and Robitaille et al. (2012), both investigated in the context of cosmic-ray propagation in Porter et al. (2017, 2018). While the model of Robitaille et al. (2012) generally describes the global Galactic emission best, it assumes local symmetry in the GC region. Freudenreich (1998) on the other hand considers a spatially asymmetric distribution of the photon field in the GC. This is why we apply the ISRF model developed in Freudenreich (1998) in our work. In the description of Freudenreich (1998), the spatial description of the energy density is given by the function

with $R_{\mathrm{end}}=3.128$ pc, $H_{\mathrm{end}}=0.461$ pc and the normalization factor $n_{0}$ determined by the energy density and $R_{s}$ the bar radial coordinate

The photon energy density as a function of the wavelength is displayed in Figure 3 and the resulting integrated energy spectrum as a function of spatial coordinates in Figure 4.

Using this model, Porter et al. (2018) investigate the influence of gamma-ray attenuation due to electron-positron pair generation in the GC. They assert an attenuation of around $10\,\%$, which would lead to an increase in the true cut-off energy by a factor of $2.1$ as compared to what is measured with H.E.S.S. The authors assume an injection of photons from a single centralized source, which neglects the fact that the photons are more likely emitted diffusely through the entire CMZ. As the propagation length through the photon field in the GC region becomes shorter if they are produced with a certain distance from the center, it is likely that the attenuation is rather lower than the expectation presented in Porter et al. (2018). We will test the influence of the photon distribution in this paper.

The simulation in this work is constrained by the following features:

1. 1.

   The minimum Lorentz factor is fixed to $\gamma_{\mathrm{min}}=1000$ so that the injected protons approximately move at the speed of light.

2. 2.

   The simulation volume is restricted to a box of $200\times 400\times 120$ pc. All particles leaving the simulation volume are lost.

3. 3.

   The maximum trajectory length $d_{\mathrm{max}}$ of injected CRs is related to the propagation time as $T_{\mathrm{max}}=d_{\mathrm{max}}/c$. A stationary solution is necessary in order to reproduce the observed diffuse gamma-ray emission (Abdalla et al., 2018). It is therefore useful to consider the CR propagation time required to leave the region of interest. Figure 6 shows the different timescales for the propagation. The diffusion timescale (green, dashed line) at relativistic energies is given as $\tau_{\mathrm{diff}}\simeq 3.3\cdot 10^{12}\cdot(E/4\ \mathrm{GeV})^{-1/3}$ s for traversing a region with radius $R=200$ pc as done in (Guenduez et al., 2018, 2019) and dominates the other timescales of advection (orange, dotted line), adiabatic losses (red, dot-dashed line) and hadronic interactions (blue, solid line). The maximum trajectory length for CR protons considering $\gamma_{\mathrm{min}}=1000$, i.e. $E\simeq 1$ TeV, yields

   $$d_{\mathrm{max}}=\tau_{\mathrm{diff}}\cdot c=5.1\ \mathrm{kpc}.$$

   (22)

   In fact, with these settings, almost all (more than 99$\,\%$) of the CRs injected from a centralized source leave the GC region. The amount of detected CRs after leaving the region of interest as a function of the maximum trajectory length is presented in Fig. 5.

4. 4.

   The source distribution $S$ is specified for the GC region as discussed in Section 3.3.

In order to optimize the choice of the minimum step length, test simulations are performed with $10^{4}$ protons of $1$ TeV energy, using $\epsilon=0$. Figure 7 shows test particle trajectories for minimum step lengths of $l_{\min}$/pc=$0.1,\,0.5,\,1.0,\,2.0$, respectively, going from top to bottom. We use this interval to be sensitive to the smallest structures in the simulation, which are the smallest MCs, in the order of $0.1-1.7$ pc. The left panel shows the $Y-Z$ distribution, i.e. the edge-on view of the GC, whereas the right panel shows the $X-Y$ distribution, i.e. the face-on view of the GC. The particles are injected directly into the GC and are propagated in the magnetic field configuration presented in Guenduez et al. (2020), with a detection taking place every $1$ pc. As the smallest MC is on the order of $0.1$ pc, this is the smallest minimum step size that we consider. As the distributions do not show any significant deviations, we use the large value $l_{\rm step,min}=1$ pc in order to save computational time.

In order to determine a useful maximum step length, the poloidal magnetic field in the ICM region is considered as it shows the least curvature and thus allows for the largest steps. While in these regions the magnetic field is quite uniform at a length of $\sim 200$ pc, the requirement of at least 25 steps before a particle leaves the simulation volume provides better reliability. This reduces the step length to $8$ pc and yields $l_{\rm step,max}=5.0$ pc. With the implementation of the SDE Module in CRPropa3.1 (Merten et al., 2017) the step size is chosen adaptive to ensure small interaction probability in each step and keep the position error due to the magnetic field integration below $10^{-4}$ per step.

Sources in the GC that have the potential to accelerate particles up to knee energies or beyond are:

1. 1.

   The central supermassive black hole Sgr A^∗  located at $l=359.94^{\circ}$ and $b=-0.046^{\circ}$ (Ghez et al., 1998; Petrov et al., 2011);

2. 2.

   The SNR Sgr A East with its center at a distance of 2 pc from Sgr A^∗ (Haardt et al., 2016);

3. 3.

   The SNR G0.9+0.1 at a distance of $\sim$130 pc from Sgr A^∗ (Abdalla et al., 2018);

4. 4.

   other SNRs such as G359.0 -0.9 (LaRosa et al., 2000), G359.10 -0.5 (LaRosa et al., 2000), G0.30 +0.04 (Kassim & Frail, 1996; LaRosa et al., 2000), G0.9 +0.1 Abdalla et al. (2018) and Sgr D (Sidoli et al., 2001);

5. 5.

   The population of unidentified pulsars that could contribute to the gamma-ray excess in the GC observed in the Fermi-LAT energy range¹¹1The pulsar (PSR) J1746-285 has already been detected in the H.E.S.S. energy range (Abdalla et al., 2018) and is consistent with emission from the pular-wind nebular (PWN) G0.13-0.11..

In this simulation, for simplicity, we use protons as the primary cosmic-ray component. We neglect the contribution from the elements of helium and higher mass numbers at this point, as the error is expected to be rather small, as in the energy range of interest, the flux is dominated by protons, see e.g. Becker Tjus & Merten (2020) for a review.

The minimum energy is set to the lower threshold energy that applies to CRPropa, which is written for ultra-relativistic particles in the limit $E\sim p\cdot c$ only. For protons, this results in a minimum energy of $E_{\min}=1$ TeV. This means that the secondary photon flux that can be considered must have a lower limit in energy $>100$ GeV. We therefore restrict the interpretation of our results to the H.E.S.S. data and refrain from comparing to Fermi-LAT data at this point.

The maximum energy is assumed to correspond to the maximum expected Galactic CR energy of $1$ PeV, i.e. the energy at the knee. The spectral index is taken from the results of the theory of stochastic acceleration in shock vicinities and therefore set to $\alpha=2.0$ for the simulations. As the SDE method allows for the reweighting of the results, this spectral index can be varied in the analysis in order to find the best fit.

The GC is enhanced in gas density and bright in background photons over a wide range of wavelengths. This abundance of background gas and photons induces different processes, i.e. hadronic interactions, electromagnetic pair production and inverse Compton scattering. These processes are taken into account in the numerical simulation approach. The related modules in CRPropa are called HI (hadronic interactions), EMPairProduction, EMPP (electromagnetic pair production), and EMInverseComptonScattering, EMIC (electromagnetic inverse Compton scattering). The module HI requires an input function concerning the gas density distribution. EMPairProduction and EMInverseComptonScattering have equivalent input parameters that require the photon density as a function of the wavelength and a scaling-grid obtained from the spatial distribution of the photon field. In addition to the GC background photon field presented in Section 2.3 (from now on called GCB), CMB and Cosmic Radio Background (CRB) (Protheroe & Biermann, 1996) radiation are considered. The latter two are homogeneously distributed and already included in the public version of CRPropa 3.1.5. All interaction modules are labeled so that secondary particles can be traced back to the generation process or the photon field. The corresponding average optical depths indicating the relevance of a particular process are presented in Figure 10. The optical depth for $HI$ is calculated for a target medium with a gas density of $10^{2}$ cm^-3 ($HI-1$) and $10^{4}$ cm^-3 ($HI-2$).

The CRPropa module Observer allows to customize the output file. Here, the user is able to construct an observer surface that detects particles only when they cross the surface. The surface may have a paraxial or a spherical shape. However, because this work is focused on secondaries, the output is not restricted to an observer sphere. All secondary particles are recorded immediately after their generation and remain an active part of the simulation.

This work takes the following information of the particles, subsequently referred to as Candidates, into account: current position, current energy, current direction of motion, trajectory length, serial number, particle type, and particle type of the parent particle, position of their source, and the generation process. Because the generation process is tagged, the gamma-ray attenuation can be considered by subtracting all gamma rays, which scatter into electron-positron pairs via the module EMPP.

In Abdalla et al. (2018), the gamma-ray count profiles are given as the number of photons detected in a certain longitudinal or latitudinal interval. In this section, we produce count maps from our simulation data: The longitudinal profiles are integrated over the latitudes from $-0.3^{\circ}$ to $-0.3^{\circ}$ and the latitudinal profiles over the longitudes from $-0.5^{\circ}$ to $0.5^{\circ}$. Because the ratio of perpendicular to parallel diffusion $\epsilon$ has a significant impact on the latitudinal and longitudinal profile, this parameter can be fixed by comparing simulations of different $\epsilon$ values with the count maps as measured by Abdalla et al. (2018). In order to reproduce the measurements, the simulation results must further be smeared to the H.E.S.S. angular resolution of 0.077^∘222The angular resolution is achieved by considering a point spread function corresponding to a $68\,\%$ containment radius.. For this purpose, the results from the simulations are treated with 2D-Gaussian smoothing with a sigma corresponding to a radius of 0.077^∘. Furthermore, the simulated data are normalized to the amplitude value at the center of the longitudinal profile. The significance of the detection is the largest at the center, so that the background noise can be neglected. We use the same binning as H.E.S.S. for our simulation data. The simulated latitudinal and longitudinal count maps for five different source distributions and four different $\epsilon$ values in each plot are presented in Figure 11.

The uppermost panels show the latitudinal (left) and longitudinal (right) profiles that are expected for the [Sgr A^∗] model, the second row shows the model [3sr], the third one [3sr+uPSR], the fourth [uPSR] and the lowest panel [hom]. In each plot, the simulation results are shown for four different $\epsilon$ values, i.e. $\epsilon=0.001,\,0.01,\,0.1,\,0.3$. In each plot, the H.E.S.S. data are shown as filled squares.

For all source distributions, $\epsilon=0.01$ and $\epsilon=0.001$ fit slightly better to the observations as compared to the results obtained from $\epsilon=0.1$ and $\epsilon=0.3$ and thus, a low influence of perpendicular diffusion is expected. More importantly, it can be seen that the source distribution has a significant impact on the latitudinal and longitudinal profiles. It appears that only the source distributions [Sgr A^∗] and [3sr] deliver a good description of the data for both the latitudinal and longitudinal profiles. Our results show that the [3sr] scenario provides the best fit. Adding the pulsar distribution does not change the result for the longitudinal profile significantly³³3the longitudinal part actually fits the data slightly better for the pulsar distribution, but the latitudinal profile is enhanced such that the prediction overshoots the data significantly in particular for positive $Z$ values. The reduced source distribution with only the central source, [Sgr A^∗], narrows down the longitudinal profile instead. For the homogeneous distribution of sources, the latitudinal profile is largely overestimated and the longitudinal profile does not fit well, in particular for positive values of $Y$. It is interesting to note that the three-dimensional propagation model comes to the same conclusion as the one-dimensional approximation used in Abramowski et al. (2016), i.e. that a three source scenario is the best fit for the data.

Thus, for the energy spectrum and the count maps that we present, we restrict ourselves to show the results for the [3sr] scenario.

For the interpretation of the longitudinal emission profile, it can be summarized that the complex and asymmetric structure, with a number of smaller peaks around the central one, is difficult and none of the models fit all peaks simultaneously. The first peak is related to the high gas density at the center and all models can reproduce it relatively well. The second observed peak location, however, is shifted with respect to the peak predicted in this work. While the peak in our simulations is caused by the dense gas of the dust ridge clouds, the observed peak rather coincides with the location of Sgr B2. These findings indicate that the gas density taken from Goldsmith et al. (1990) and used in our simulations is not compact enough. If the mass that is connected to Sgr B2 is distributed onto a smaller volume with a higher density, the second peak at the location of Sgr B2 could be enhanced. For clarification, further observations by instruments such as CTA will provide significantly better spatial resolution and cover a broader energy range. Improvements on the gas density measurements will also help to better quantify the fluxes.

The peak at $Y\sim 60$ pc in all source setups is due to the compact MCs dust ridges A - F. Their compactness seems to be sufficient and reasonable for the source setups that provide the best fit, i.e. [Sgr A^∗] and [3s]. The enhanced emission at $Y\sim-100$ pc to $Y\sim-30$ pc is not reproduced by this model.

For the future, for a reliable identification of a full source setup, the detailed knowledge on the MCs is crucial. The observed third peak is not seen in the simulation results because no relevant sources or MCs are located between $l\sim-40$ pc and $l\sim-80$. The MC Sgr C at $\sim 82$ pc has a radius of approximately 1.7 pc and alone would not be sufficient to explain the third peak. Thus, in the given vicinity, additional MCs or unresolved active sources could be present. Alternatively, the size of Sgr C and thus the mass could be underestimated in the observational results. For a convenient comparison, Tables 2 and 3 list the resulting chi-square tests adapting the null hypothesis that the observed data are described by the simulation.

To summarize, a more realistic gas distribution is necessary in order to explain all of the features seen in the gamma-ray data. The H.E.S.S. collaboration has developed an empirical model for the gas distribution, see Abdalla et al. (2018) and references therein, based on the gamma-ray distribution. This, however implies the knowledge of the source distribution. As both cannot easily be disentangled from each other, increasing the knowledge on the gas distribution itself would be of high importance.

The energy spectrum of the CMZ as measured by H.E.S.S. is of particular interest, as it is the first to extend to almost $100$ TeV photon energy. For a hadronic interpretation, that implies the existence of cosmic rays with $\sim$ PeV energies or higher, thus providing the first measurement of a possible PeVatron in the Galaxy.

The simulation data that we present in the following are integrated over the so-called Pacman region, which corresponds to a ring-like, symmetric region of a radius of 0.45^∘ around Sgr A^∗, where the inner 0.15^∘ and a section of 66^∘ are excluded (Abramowski et al., 2016). A schematic view of the Pacman region together with seven more regions of detection that form the basis of the spatially resolved signal presented in Abramowski et al. (2016) and Abdalla et al. (2018) are shown in Fig. 12 .

In the simulation, we inject a cosmic-ray energy spectrum $Q\propto E^{-2}$ and we reweight this spectrum to different energy behaviors $Q\propto E^{-\alpha}$ with $\alpha\in$[2.0, 2.05, …, 2.4]. The energy interval binning is adjusted to match H.E.S.S. data presented in Abramowski et al. (2016). The normalization factor is determined by assuming that the sum of the data bins equals the sum of simulated data bins, i.e. the surface below the data points is the same. The photon energy spectra resulting from the simulations performed here are shown in the four panels of Fig. 13 - the uppermost panels are for $\epsilon=0.001$ (left) and $\epsilon=0.01$ (right) and the lower ones for $\epsilon=0.1$ (left) and $\epsilon=0.3$ (right). In each panel, the result is shown for five different spectral indices, i.e. $\alpha=2.0$ (circle), $\alpha=2.1$ (pentagon), $\alpha=2.2$ (triangle up), $\alpha=2.3$ (triangle down) and $\alpha=2.4$ (diamond). The H.E.S.S. data are shown as well (red squares). The best fit to the data is achieved for the combination $\epsilon=0.001$ and $\alpha=2.3$, where we base the goodness of the fit on chi-square tests under the null hypothesis that the observed data are described by the simulation. The results of this test are presented in Table 4. An increase of $\epsilon$ seems to decrease the best-fit spectral index. The simplistic assumption by Abramowski et al. (2016), i.e. $\alpha=\alpha_{\gamma}-0.1$, requires a proton spectral index of $\sim$ 2.4, which is slightly steeper than what we find.

Different leptonic and hadronic processes can contribute to the gamma-ray energy spectrum as discussed above. Recently, Porter et al. (2018) predict a significant attenuation of high-energy photons due to their interaction with the ambient photon field. In their calculation, they assume gamma-ray emission from one centralized source. By contrast, this work considers the entire CMZ, and the results of this chapter can be used to test the previous prediction. Moreover, this work takes into account gamma-ray production via inverse Compton scattering by secondary electrons.

Figure 14 displays the total differential photon flux from the best-fit simulation ($\epsilon=0.001$ with source index $\alpha=2.3$) as well as the contribution from proton-proton interactions ($HI$, circle) and from electromagnetic inverse Compton ($EMIC$, pentagon).

The figure shows that the dominant contribution above $100$ GeV photon energy indeed comes from the hadronic interaction channel. Inverse Compton (IC) scattering has a relatively high contribution at $\sim 100$ GeV yet due to the steepness of the spectrum, the relative contribution of IC decreases rapidly with energy. It is possible that primary electrons still have an impact on the photon spectrum until TeV energies near the center of the Galaxy, but this question is not investigated in this paper, as we do not propagate the electrons here.

Figure 14 even shows the fraction of the originally produced photon flux that is attenuated through gamma-gamma interaction (EMPP, triangle up).

The attenuation reduces the flux by $1\%$ to $10\%$ and occurs mainly at energies below $\sim 7$ TeV. This attenuation is significantly lower compared to the results presented in Porter et al. (2018). A reason could be that the photons in our simulations are produced in a relatively broad region in the CMZ, while Porter et al. (2018) assume a single central source. Thus, the photons in our simulations traverse a smaller column than in Porter et al. (2018), which leads to a smaller attenuation effect.

We also perform a fit to both photon, neutrino and electron data from the simulation. We fit the following function for all three spectra:

with $i=\gamma,\nu,e$. Using a source emission $\mathrm{d}N/\mathrm{d}E\propto E^{-\alpha}$ with $\alpha=2.3$ and strong parallel diffusion ($\epsilon=0.001$) the best-fit adjustment for the gamma-ray flux in Fig.  14 ($\epsilon=0.001$) delivers $N_{0}^{\gamma}=(1.37+\pm 0.03)\cdot 10^{-9}$ eV^-1 cm^-2s^-1, $\alpha_{\gamma}=2.20\pm 0.02$ and $E^{\gamma}_{\mathrm{cut}}=42.4\pm 17.6$ TeV. It should be noted that we do not include a cutoff in our simulations. The cutoff in this spectrum is consistent with being above the last significant simulation data point. A discussion of the cutoff in the data can be found in Abramowski et al. (2016). We refrain from discussing this point in detail as we consider the data not to deliver significant information about a cutoff at this point. With CTA data, this will change in the future.

The results for the neutrino–spectra are discussed in Section 4.4.

Figure 16 shows the count maps of the photons that are produced in the GC via the interaction of the $10^{6}$ induced cosmic-ray protons with the ambient medium for $\epsilon=0.01$ (upper panel) and $\epsilon=0.001$ (lower panel). We find that the center emission is very concentrated and declines strongly toward high Galactic latitudes, basically disappearing at values of $|b|>0.2$. The decline toward larger Galactic longitudes is slower; the signal extends up to $l\sim 1.0^{\circ}$ and disappears at $l<-0.5^{\circ}$ at negative latitudes.

In order to properly compare the simulated spatial structure of the gamma-ray map to the emission observed by H.E.S.S., Fig. 17 shows our results using Gaussian smoothing that corresponds to the H.E.S.S. angular resolution of 0.077^∘. The figure also shows the $>4.5\,\sigma$ and $>7.5\,\sigma$ significance levels for the H.E.S.S. detection as orange and red lines, respectively. For comparison, the predicted lines of equal gamma-ray counts are added in black. While we do not make a prediction of the normalization itself, it can be seen that the shape of the observed distribution is matched well by the simulation data.

The region with the highest significance of detection, enclosed by the red line, in general shows quite good agreement with our simulation data. Even in this representation, the simulation with $\epsilon=0.001$ fits slightly better than $\epsilon=0.01$ or larger epsilon values. As this is consistent with the findings in the previous subsections, we show the following plots for $\epsilon=0.001$ only. The discrepancy is more significant for positive longitudes, consistent with what we found in the results shown in Fig. 11. One reason for this could be existing Molecular Clouds that are not considered in this work due to the lack of data. The existence of such additional MCs has been discussed by Oka et al. (2001), but it remains difficult to detect a full sample at this point. In particular, the cataloged masses show inconsistencies with other observations, see e.g. Kauffmann et al. (2017); Ferrière et al. (2007); Goldsmith et al. (1990), which might be due to an overestimation obtained from the usage of virial masses that rather represent an upper limit than an exact measurement. Neglecting the Galactic background diffuse large-scale contribution, as described by Abramowski et al. (2014), can also partially explain the discrepancy.

Figure 18 shows the gamma-ray map with a smearing of $0.03^{\circ}$ which corresponds to the approximate resolution of CTA south for $>$ TeV energies (Acharya et al., 2013). At such a high resolution, the structure of the gas can be identified well. Beyond the central source, there are two specific localized regions that show up at the CTA resolution which are less pronounced for the H.E.S.S. data: at a longitude of $b\sim 0.5$, the six dust ridge clouds (A-F) are illuminated. At a longitude of $b\sim-0.5$, the cloud Sgr C shows up. It can be expected that CTA will discover even more molecular clouds, as the analysis of the comparison of our simulation data with the H.E.S.S. results in previous subsections indicates that there is still some significant amount of MC gas that could exist, but is not identified at this point.

The count map of co-produced neutrinos via proton-proton interactions is presented in Fig. 19. The map shows all neutrino flavors, but considering neutrino oscillations from the ratio at a production of $(\nu_{e}:\nu_{\mu}:\nu_{\tau})=(1:2:0)$. With neutrino oscillations, the ratio is $(\nu_{e};\nu_{\mu}:\nu_{\tau})=(1:1:1)$ so that without considering the normalization, the picture shown does represent the flux of each neutrino flavor individually. This is important, since the best pointing for neutrinos is received for charged current interactions of muon neutrinos in the Antartic Ice instrumented by the IceCube Neutrino Observatory.

While the count map is quite similar to the gamma-ray one, it is not exactly the same. There are two reasons for this: first of all, Inverse Compton processes contribute to the gamma-ray map. For the hadronic part, neutrino- and gamma-ray maps should basically be close-to identical at production. Afterwards, gamma rays are absorbed by gamma-gamma interactions, which changes the count map with respect to the one at production. These two effects alter the gamma-ray map with respect to the neutrino one.

A smeared version with an angular resolution of $0.1^{\circ}$ is shown in Fig. 20. This resolution corresponds to the one of IceCube for through-going tracks at $\sim 2$ TeV (Aartsen et al., 2017). Considering only kinematic effects, this resolution is already obtained at $0.1$ TeV, and the resolution is in general increasing with increasing energy, as a result of the enhanced forward scattering at higher energies (Learned & Mannheim, 2000; Becker, 2008).

As can be seen from Fig.  19, it will be difficult to see sub-structures, even if the number of detected neutrinos would be enough to have a spatial resolution. So, IceCube right now will have difficulties in resolving these structure, with future observatories and detections at $\sim 100$ TeV$--$ PeV, the resolution becomes better and the substructures might emerge, if the event rates are large enough (see further discussion below).

The neutrino energy spectrum is calculated in the same way as the gamma-ray energy spectrum (see Section 4.2). In Fig. 21, the simulated energy spectrum of muon neutrinos and anti-muon neutrinos is shown, derived from the flux of all neutrinos by considering neutrino oscillations from a flavor ratio at production of $\left.(\nu_{e}:\nu_{\mu}:\nu_{\tau})\right|_{\rm source}=(1:2:0)$ to the one at Earth of $\left.(\nu_{e}:\nu_{\mu}:\nu_{\tau})\right|_{\oplus}=(1:1:1)$. We show the simulation results for primary spectra of $\alpha=2.1$ (squares) and $\alpha=2.3$ (circles).

The power-law fit (Equ. 42) for the neutrino flux induced by a source population with a primary energy spectrum of $E^{-2.1}$ and $E^{-2.3}$ yields the parameters given in table 5.

As for the gamma-ray spectrum, the cutoff-energy is not an effect of a cutoff in the simulation, but it is located at energies that are higher than the simulated, significant data points. It is basically a lower-limit for a cutoff in the spectrum.

An IceCube $90\,\%$ C.L. upper flux limit for the position of Sgr A^∗ ($\sin(\delta_{SgrA*}=-29.01^{\circ})=-48.5$) is given in Aartsen et al. (2020b) for a spectral index of $\alpha=$2.0 (dashed, grey line in Fig. 21). The lower energy threshold of the analysis is relatively high, as the GC region is located at a declination quite far above the horizon at the South Pole, i.e. $\delta=-30^{\circ}$. This implies that there is a large background of atmospheric muons, which is reduced by applying the high-energy cut. A different method is to only consider events that start inside of the detector. This approach was first developed for high energies (”high energy starting events”, HESE) and lead to the first detection of the diffuse neutrino flux (Aartsen et al., 2013). The extension of this method to lower energies followed shortly after (Aartsen et al., 2015) and provides another channel of detection of the Galactic Center region.

The expected neutrino flux from the Galactic Center emission can be converted into a number of neutrinos to be detected in different analyses with IceCube based on the effective areas. Here we consider only the muon-neutrinos and their effective area, as the direction uncertainty for the other flavors is to large. For the HESE sample it is published for the declination bands $[0^{\circ};-30^{\circ}]$ and $[-30^{\circ};-90^{\circ}]$ concerning the point source search (Aartsen et al., 2020a). As the Galactic Center lies at the boundary of these two declination bands $\delta\sim-29^{\circ}$, we use these two effective areas to calculate the rate of neutrinos according to

with $\Delta t_{\rm obs}$ as the observation time. Based on this estimate, we provide the range of neutrinos. We receive a total number of neutrinos between $2\cdot 10^{-4}-0.08$ for the full 10 years. For the Medium Energy Starting Events (MESE) effective area, published in Aartsen et al. (2015) as supplemental material, we use the zenith angle bin $-0.2\leq\cos(\theta)\leq 0$, where $\theta=90^{\circ}-\delta$ is the zenith angle. Assuming an observation time of 10 years we receive $\sim 0.03$ events⁴⁴4Note that the published MESE analysis is based on 4 years of data, so the expected number of events for the published analysis is even smaller, $\sim 0.012$. Thus, the sources in the Galactic Center itself are not observable with IceCube at this point.. It should be noted that we do not include the diffuse sea of background cosmic rays that propagate into the Galactic Center from outer regions of the Galaxy, so our numbers are rather considered to be a lower limit.

KM3NeT will have a better detection potential due to its location at the opposite side of Earth compared to IceCube. The sensitivity corresponding to the median $90\,\%$ confidence level (C.L.) upper-flux limit for an unbroken power-law distribution proportional to $E^{-2}$ using 6 years of data results in $\Phi^{90\%}_{\nu_{\mu}+\overline{\nu}_{\mu}}=0.297$ eV cm^-2 s^-1 (Aiello et al., 2019). Considering the extended source analysis of KM3NeT by Ambrogi et al. (2018), the minimum detectable flux is reached at $\sim$ 60 TeV. For lower energies, the flux limit increases rapidly. For energies $>$ 60 TeV, the median angular resolution improves from $\sim$ 0.13^∘ up to 0.07^∘ at $10^{5}$ TeV. Because the angular resolution is smaller than the extent of the Pacman region, the upper limit analyzed for the point source Sgr A^∗ can be extrapolated to a more extended region by considering the median angular resolution following

Here, $\Delta\delta_{\mathrm{Pacman}}$ denotes the extent of the region of interest in degrees and $\delta_{\mathrm{MA}}(E)$ the median angular resolution as a function of the energy. The median angular resolution and KM3NeT using 6 years of data and of IceCube using 7 years of data is presented in Adrián-Martínez et al. (2016) and Aartsen et al. (2017), respectively. Figure 21 shows the $90\,\%$ C.L. upper-flux limit of IceCube and the corresponding sensitivity of KM3NeT. The limit calculation follows Ambrogi et al. (2018), who calculate the sensitivites to extended sources, where the Pacman region has an extension of $\sim 0.5^{\circ}$.