MAGIC observations of the diffuse $\gamma$-ray emission in the vicinity of the Galactic Centre

Aims: $\gamma$ rays can be used as a tracer in the search of sources of Galactic cosmic rays (CRs). We present deep observations of the Galactic Centre (GC) region with the MAGIC telescopes, which we use for inferring the underlying CR distribution. Methods: We observed the GC region for ${\approx}100$ hours with the MAGIC telescopes from 2012 to 2017, at high zenith angles (58-70~deg). This implies a larger energy threshold, but also an increased effective collection area compared to low zenith observations. Using new software, we derive instrument response and background models, enabling us to study the diffuse emission in the region. We use pre-existing data of the gas distribution in the GC region to derive the underlying distribution of CRs. Results: We obtain a significant detection for all four model components used to fit our data (Sgr~A*, ``Arc'', G0.9+0.1, and an extended component for the Galactic Ridge). We find that the diffuse component is best described as a power-law with index 2 and an exponential cut-off at around 20~TeV with the significance of the cut-off being only 2~$\sigma$. The derived cosmic-ray profile hints to a peak at the GC position, with a measured profile index of $1.2 \pm 0.3$, supporting the hypothesis of a CR accelerator at the GC. We argue that the measurements of this profile are presently limited by our knowledge of the gas distribution in the GC vicinity.


Introduction
The Galactic Centre (thereafter GC) is one of the most extraordinary regions in the very high energy (VHE, >100 GeV) γ-ray sky, containing a rich variety of sources, capable of accelerating charged particles and thereby producing VHE γ-ray emission (van Eldik 2015; Aharonian et al. 2006a;Collaboration & Aharonian 2006). If those γ rays are produced in hadronic interactions, their sources may also contribute to the overall cosmicray "sea", filling the region with energised charged particles. This sea of cosmic rays (CRs) is believed to be responsible for the diffuse Galactic plane emission detected by EGRET (Hunter et al. 1997) and later studied by Fermi/LAT (Abdo et al. 2010), as well as the central γ-ray ridge detected with all major Imaging Atmospheric Cherenkov Telescopes (IACTs, Collaboration & Aharonian 2006;Archer et al. 2016;Ahnen et al. 2017a). The detected γ rays originate primarily from deep inelastic nuclear interactions of the high-energetic particles with the surrounding gas, which can also be traced via its emission in the radio band (CS emission, Tsuboi et al. 1999;CO emission, Oka et al. 1998). Thus deep γ-ray observations in the TeV energy range can be used to infer the underlying CR density, provided that the gas distribution is reliably measured.
This idea has been applied using observations of the central ≈ ± 1 • region (in latitude) of the Galaxy (corresponding to ≈300 pc at a distance of 8.5 kpc) with the High Energy Stereoscopic System (H.E.S.S.), suggesting an inhomogeneous distribution of CRs around the GC (Aharonian et al. 2006b;H.E.S.S. et al. 2016. Furthermore, those data indicate that the central super-massive black hole of our Galaxy may itself accelerate particles to PeV (10 15 eV) energies.
Observing the GC region with the MAGIC telescopes is only possible at large zenith distances (> 58 • ), implying an increased optical thickness of the atmosphere and larger distance to the shower maximum, which leads to a stronger dilution of the Cherenkov light. On the other hand, for geometric reasons, these data benefit from an increased collection area (≈1 km 2 , Ahnen et al. 2017a) and thus a boost of sensitivity at energies above several TeVs. Ahnen et al. (2017a) have already presented the first part of the MAGIC GC observation campaign in 2012-2015, consisting of ≈70 hours of exposure. Apart from the detection of multi-TeV emission from the Galactic plane, those data suggested the presence of a new VHE source close to the so-called radio "Arc" in the GC vicinity, detected also by other telescopes (Archer et al. 2016;H.E.S.S. et al. 2018). We continued observing the GC during 2015 -2017, adding up to a total exposure time of ≈100 hours. In addition, we employ the recently developed 2D likelihood analysis package SkyPrism (Vovk et al. 2018), which provides a more sensitive analysis of extended sources for MAGIC. This paper is organised as follows. In Sect. 2 we describe the MAGIC observations, the data selection strategy and the data analysis techniques. In Sect. 3 we provide a detailed description of the various analysis steps and report on their results. In Sect. 4 we estimate the biases, stemming from our background modelling and analysis approaches. Finally, Sects. 5 and 6 provide a general discussion and a summary of the obtained results.

MAGIC observations of the Galactic Centre region
The MAGIC (Major Atmospheric Gamma Imaging Cherenkov) telescopes are two 17 m diameter IACTs, located at a distance of 85 m from each other at an altitude of 2200 m a.s.l. at the Roque de los Muchachos Observatory on the Canary Island of La Palma, Spain (28 • N, 18 • W). The telescopes record flashes of Cherenkov light produced by Extensive Air Showers (EAS) initiated in the upper atmosphere by γ-ray photons within a field of view of 1.5 • (in radius). Both telescopes are operated together in a coincidence trigger stereoscopic mode, in which only events simultaneously triggering both telescopes are recorded and analysed (Aleksić et al. 2016). For low zenith distance observations and for E > 220 GeV, the integral sensitivity of MAGIC is (0.66 ± 0.03)% in units of the Crab Nebula flux (C.U.) for 50 hours of observation (Aleksić et al. 2016). At larger zenith angles (above 60 • ), at which the GC is observable at La Palma, the MAGIC collection area reaches 1 km 2 at the energies above 10 TeV (Ahnen et al. 2017a).
The MAGIC observations of the GC region, used for this study, were carried out between April 2012 and May 2017 (typically from March to July). The observation time after quality selection cuts is ≈100 hr. The data were taken in the so-called "Wobble" mode (Fomin et al. 1994), centred on Sgr A* with a pointing offset of 0.4 • .
The acquired data have been processed with the standard MAGIC Analysis and Reconstruction Software (MARS, Zanin et al. 2013). These steps include selection of good data quality periods, the low-level processing of the camera images, the calculation of image parameters and the reconstruction of the energy and incoming direction of the primary γ rays. The last step is performed using the machine learning technique Random Forests (Albert et al. 2008), where Monte Carlo events and real background data are used for training, and which is also used for background suppression through event classification.
At the quality selection step, we removed events corresponding to known temporary hardware issues or recorded during bad weather periods. As the full data sample contains some nights with weak moonlight, we also exclude time periods with sky brightness 2.5 times the typical dark night sky brightness. Data recorded with higher sky brightness lead to biases when analysed with the same pipeline that is used for data that has been recorded during dark nights, but can in principle still be used when accepting a higher energy threshold when using special analysis settings (see Ahnen et al. 2017b). In order to select only good quality data we applied cuts on several instrument and weather related parameters. Those include the mean photomultiplier currents, the event trigger rate, the number of stars detected by the MAGIC star-guider cameras and the cloud-/aerosol-induced absorption in the telescope's field of view. The latter is derived from measurements with an infra-red pyrometer and a LIDAR system Fruck et al. 2014).

Data analysis
The high-level analysis of the acquired data has been performed with the set of utilities available in the MAGIC SkyPrism package (Vovk et al. 2018). SkyPrism provides a set of routines, aimed at 2D spatial likelihood analysis of MAGIC data, which allows us to self-consistently account for sources of complex morphology within the field of view. Furthermore, it also incorporates algorithms for computing the instrument point spread Article number, page 3 of 10 A&A proofs: manuscript no. MAGIC_GalDiffuse function, exposure and background maps, which we intensively use throughout the analysis presented here.

Diffuse background model construction
The wobble observational scheme (Fomin et al. 1994), used here, provides a simultaneous background estimate from comparison of several adjunct (albeit exposed to slightly different sky regions) positions. Still the applicability of this approach is limited by the distance between the wobble pointings and the target, which in our case was d w = 0.4 • . In practice, a source with an extension comparable to or larger than twice the offset distance (0.8 • ), in the same direction, will still (partially) contribute to the background measured in the same camera region. The diffuse emission of the Galactic plane clearly exceeds this limit, so a special treatment is necessary.
Such contamination of the background map with excess γray signal can be avoided if contributions from locations close to known γ-ray sources are excluded during map construction (Vovk et al. 2018). We thus mask out camera regions that correspond to the Galactic plane and use the rest of the camera (free from known/expected sources) to estimate the corresponding background. This approach, implemented in the SkyPrism package, allows us to significantly reduce the background bias.
Due to the remaining limitations of our background modeling technique and the used wobble scheme, we can not completely remove the possible bias in the background estimation from a source as extended as the diffuse Galactic plane emission. Our estimates described in detail in Sect. 4.1 and illustrated in Fig. 7 demonstrate that in the l = [−1.5 • ; 1.5 • ] longitude range along the plane the remaining bias is within 10-20% of the assumed local source luminosity throughout the Galactic plane and does not exceed 30-50% (equivalent to ≈3% of background) in the outskirts. At the same time the flux bias is rather constant across the plane and its variations, averaged in the latitude range b = [−0.2 • ; 0.2 • ] do not exceed ≈1% of the background, which corresponds to 10-20% of the estimated source flux for a highly extended source.

MAGIC view of the Galactic Centre region
The large (> 58 • ) zenith angle GC observations imply an increased energy threshold of ≈1 TeV, though an analysis is also possible at even lower energies, given that the studied source is bright enough (Ahnen et al. 2017a). Thus we have performed the spectral analysis in the 400 GeV -50 TeV energy range; the morphology study of the diffuse emission in the GC region was done above γ-ray energies of 1 TeV.
The sky map (above E = 1 TeV) of the GC vicinity, produced with the described diffuse background estimation scheme, is shown in Fig. 1. The Galactic plane is visible over 2 • across the image. The significance of this detection, computed using the CS radio emission profile (Tsuboi et al. 1999) as an approximation (see Sect. 3.3 and 3.4 for details of the method), results in a ≈17 standard deviations (σ) incompatibility with the null hypothesis of background and point sources only. Other sources visible in the image are Sgr A*, G0.9+0.1, and the so-called "Arc", detected with significances of ≈48 σ, ≈11 σ and ≈6.4 σ respectively, propagating also uncertainties of the background and exposure model. Galactic coordinates at energies above 1 TeV, smeared with a kernel resembling the MAGIC PSF. The pre-trial statistical significance of regions with excess counts is highlighted by green contours at the levels 5σ and 3σ. The (smeared) MAGIC PSF is indicated by 39% and 68% containment contours. The white contours show radio line emission from CS molecules, tracing dense gas (Tsuboi et al. 1999).

Galactic plane brightness scan
The CR distribution profile in the GC surroundings should roughly correspond to the brightness distribution of the detected γ-ray emission. Previous measurements have already shown evidence that the 100 GeV brightness of the Galactic plane is peaking towards the Sgr A*, indicating a concentration of cosmic-rays ( The cosmic-ray distribution around the GC can be inferred by solving the integral equation where S (x, y) is the image plane γ-ray brightness distribution, ρ gas and ρ CR are the number densities of the gas and CRs respectively and A is a factor, which takes into account the protonproton interaction cross-section, distance to observer and additional constants. A proper solution of this equation requires knowledge of the full 3D gas density distribution ρ gas (x, y, z), which is challenging to obtain. Indeed, while in the image plane (x, y) the resolution of radio surveys reaches ≈0.01 • (for instance the CS J = 1 − 0 emission radio survey of Tsuboi et al. 1999) equivalent to 1-2 pc scale, the line-of-sight distance z can hardly be obtained with an accuracy better than several tens of parsecs. For this reason we solve an approximate expression of Eq. 1: which splits the problem into the projected gas (directly inferred from the radio data) and cosmic-ray distributions P gas (x, y) and P CR (x, y) respectively. To avoid degeneracy, we just consider radially symmetric cosmic-ray profiles ρ CR ∝ r −α (with r being the distance from the GC) and their projections onto the image plane. These simplifications result in a certain bias of our measurement, which we quantify in Sect. 4.2.
First we test whether a homogeneous cosmic-ray distribution ρ CR = const is consistent with the MAGIC data. Sgr A* and G0.9 point sources, the "Arc" source, and the diffuse emission model S (x, y) (computed from the CS emission maps with ρ CR = const), it results in χ 2 /d.o. f. ≈ 69/46 degrees of freedom, equivalent to ≈2.4σ disagreement of the data with the model.
In order to investigate if the MAGIC data are in a better agreement with a ρ CR const type cosmic-ray profiles we use a grid search for the optimal value of α using the profile shown in Fig. 2 and estimated the cosmic-ray density P CR (d) = S (d)/P gas (d) as function of the projected off-centre distance d using a full maximum likelihood fit of the measured γ-ray brightness around the GC.
The results of the likelihood-profile scan are shown in Fig. 3, where the density ρ CR is converted to the cosmic-ray energy density w CR using the same procedure as in H.E.S.S. et al. (2016) (formula (2) of the methods section): where w CR is in eV/cm 3 , η N ≈ 1.5 accounts for nuclei heavier than hydrogen, in both CRs and in the target gas. For estimation of the H 2 target mass M based on the CS radio maps we used the procedure described by the authors in section 4.2 of Tsuboi et al. (1999): where the excitation temperature T ex of CS is 30 K, T MB dv is the measured velocity integrated antenna temperature, A is the area of the GC region, µ(H 2 ) is the mass of the hydrogen molecule and X(CS ) = 10 −8 is the relative abundance of CS in H 2 clouds.
In the scan we have tested two different scenarios: (1) the ρ CR ∝ r −α profile dominates the cosmic-ray density in the region and (2) the peaked cosmic-ray profile is found on top of an underlying homogeneous density, so that ρ CR ∝ r −α + const. To account for the background and telescope exposure uncertainties, we have generated 50 random exposure and background maps representing the uncertainty range of our reconstruction method and repeated the scan for each of them. We then averaged the resulting likelihood values in each bin of the scan phase space, which is equivalent to marginalization over these random map representations. The same technique of propagating the uncertainties of the background and exposure models has also been applied throughout the spectral analysis (Sec. 3.4). The resulting averaged values were then used to compute the confidence contours, shown in Fig. 3. In addition to α − w CR combined confidence contours, the right panel of the same figure shows the marginalised uncertainties for the power-law index α for both scenarios (1) and (2), with darker and lighter colours respectively.
As it is illustrated in Fig. 4, also the full likelihood fit to the obtained sky map is not consistent with the ρ CR = const assumption. To perform this fit, we have first split the CS emission map of Tsuboi et al. (1999) (integrated over the radial velocity) into a sequence of concentric rings, with their own normalisation factors. Since the CS emission is highly peaked toward the Galactic plane, this is effectively equivalent to a longitudinal split of the Plane inside the b = [−0.1 • , 0.1 • ] stripe. When computing the normalisations of the rings, we've also included the Sgr A* point source, as well as G0.9+0.1 and the "Arc" source to the fitted model. The resulting cosmic-ray density profile, shown in

Spectral analysis of the detected sources
To compute energy spectra of the detected sources in the MAGIC field of view -including the diffuse Galactic plane emissionwe have used the the SkyPrism package. The spatial model used for the fit includes three point sources (Sgr A*, G0.9+0.1 and the Arc source at RA=17:46:00, Dec=-28:53:00) and the velocity integrated CS map, re-scaled with the ρ CR ∝ r −1.2 bestfit cosmic-ray profile. The fit has been performed in the energy range from 400 GeV to 50 TeV with two methods -energy binwise (7 logarithmic energy bins) and assuming a certain spectral shape model for each of the sources. In the latter case we applied a forward folding procedure considering the energy migration matrix during the fit. All spectra were generated from the general form which can result in a power-law, log-parabola and power-law with cut-off spectral shape depending on the choice of β and E cut . The normalisation energy for all the sources was set to E 0 = 2 TeV, keeping the correlation between the spectral parameters minimal.
To obtain the best fit parameters of the assumed spectral models we perform a maximal Poissonian likelihood fit to the energy-binned MAGIC sky maps. To ensure a good accuracy of the estimated uncertainties, we have additionally used the Markov Chain Monte Carlo (MCMC) sampler emcee on the parameter space (Foreman-Mackey et al. 2013). The uncertainties of the exposure and background models that were derived from Monte Carlo simulations and data regions more than 0.3 • off the Galactic plane, were propagated to the final results by processing 60 random representations through the MCMC sampler and merging the samples. The best-fit values for the detected sources, obtained through this fit, are given in Tab. 1, along with the corresponding errors and detection significances. The obtained spectra (data points and fit results) are shown in Fig. 5. The data points are not the result of spectral unfolding, but spillover corrections based on the energy migration matrix and the fitted spectral shape were applied. The obtained MAGIC spectrum is consistent with the earlier estimate of the Galactic ridge Spectral Energy Density (SED) (H.E.S.S. et al. 2018), as displayed in Fig. 6. A likelihood ratio test comparing the model for the diffuse component with cut-off to a pure power-law results in the ≈2 σ preference for the cut-off for the MAGIC data set. We estimate the systematic uncertainties, arising from uncertainties on the energy and flux normalization scales, following the procedure discussed in Ahnen et al. (2017a) (based on a detailed study by Aleksić et al. 2016). The resulting estimates are indicated by gray arrows in Fig. 5 and Fig. 6, where the vertical arrows indicate the effect of the flux normalization errors at different energies and the horizontal or inclined arrows indicate the effect of the energy scale uncertainty.

Bias from the background modelling
In order to quantify the bias resulting from the background estimation, we used a simplified simulation of the background map based on the initial assumption on the extension and brightness of the sources in the GC region. In this simulation we assume that the true signal measured by MAGIC consists of five contributions: extended gas emission (assumed to be traced by the CS map (Tsuboi et al. 1999)), point-like Sgr A*, the "Arc" source (Archer et al. 2016;Ahnen et al. 2017a;H.E.S.S. et al. 2018), point-like G0.9+0.1 and isotropic background. The relative normalisations of these components are taken from the previous analysis of Ahnen et al. (2017a); we used our best fit results for a cross-check. This composite image is then used to sample the photons in the telescope camera coordinates as a function of pointing azimuth and zenith, following the MAGIC pointing during the GC observations. The resulting event list is supplied to the background estimation routine for a comparison of the reconstructed vs. assumed background.
Based on the results, illustrated in Fig. 7, we expect the bias to stay below 2% in units of background flux nearly everywhere in the sky-maps. This translates to a bias of the measured flux of the diffuse emission of more than 40% in some regions along the edges of the Galactic plane, but less than 30% for the brighter regions along the plane and at the centre. Still, using the CS map as an approximation for the γ-ray emission, the total bias on the integral flux of the Galactic plane component is estimated to be in the range 7-12%.  with the spectral fit results of the sources, detected in the MAGIC field of view. The spectrum type acronyms stand for: PL -powerlaw and PLC -cut-off power-law. Normalisation factor N is given in units of 10 −25 [ph/(cm 2 s eV)]; normalisation energy is set to E 0 = 2 TeV for all the sources. The curvature parameter β = 0 in all cases. For the sources fitted with the power-law model the values of E cut are not given. All uncertainties correspond to a 68% confidence interval.  For the components corresponding to Sgr A*, the "Arc" and the CR/MC component a power-law shape with exponential cut-off has been used while G0.9+0.1 can be described with a simple power-law. The error bars and bands were computed from the MCMC samples and correspond to 68% confidence range. No spectral unfolding has been applied to the data points, but the effect of spillovers due to energy migration has been corrected for, based on the spectral shapes. Gray arrows indicate the size (length) and direction (orientation) of the SED shifts due to the systematical uncertainties in the analysis (see Sect. 3.4 for details).

Bias and uncertainty from the assumed gas distribution model
An accurate modelling of the γ-ray emission from GC region requires detailed knowledge of the gas distribution in three dimensions. The angular resolution of MAGIC is better than 0.1 • , which translates into ≈15 pc at 8.5 kpc distance from the Earth. As a result, MAGIC can map the profile of γ-ray emission at the projected distances of tens of parsecs from the position of the central supermassive black hole (SMBH); a conversion of this profile to the cosmic-ray density distribution then naturally requires that the line-of-sight distances to the gas clouds in the region are known with similar or better accuracy.  This requirement is very difficult to fulfil in practice. At larger distances the locations of the gas clouds are generally inferred from their kinematics, assuming a certain model of gas orbital motion (e.g. Sofue 1995;Nakanishi & Sofue 2003). In the vicinity of the GC, however, this approach can no longer be applied, as the inability to put the source in front or behind the black hole image plane at the scales of tens of parsecs leads to degeneracy in the calculations. The required information -to a certain degree -can be reconstructed using measurements of line-of-sight absorption of the molecular cloud emission, which provides the necessary line-of-sight position estimates (Sawada et al. 2004). Nevertheless, the line-of-sight locations of separate clouds in the GC region can hardly be reconstructed with an accuracy better than ≈50 pc.
In the absence and therefore negligence of the line-of-sight information, Eq. 1 naturally simplifies to Eq. 2, which works only with the projected gas and cosmic-ray densities. Depending on the real distribution ρ gas (x, y, z), the transition Eq. 1→ Eq. 2 may bear an oversimplification, resulting in a biased cosmic-ray profile ρ CR (r).
We do not attempt to reconstruct a fully realistic structure of the central ≈200 pc of our Galaxy, which is rather complex (Ferrière et al. 2007). Still, in order to quantify the associated bias in our analysis, we reconstruct the 3D gas distribution in the GC region based on the measurements of Sawada et al. (2004) and Tsuboi et al. (1999). This reconstruction then allows us to compare the profiles obtained accounting for or neglecting the line-of-sight information.
To perform the reconstruction, we used the fact that the 2D CS gas emission images of Tsuboi et al. (1999) are complemented by the radial velocity v rad information, which already provides the line-of-sight information in an indirect way. A mapping of v rad to the missing z coordinate is found in Sawada et al. (2004), where it is given in the (x, z) projection (in our notation). Hence, the measured intensity of the radio emission I(x, y, v rad ) can be mapped to a full 3D cube I(x, y, z) through a relation v rad (x, y) ↔ v rad (x, z). The obtained cube gives an approximate picture of 3D gas distribution in the innermost ≈400 × 400 × 100 pc of the GC region.
Despite the limitations of this approach, based on this 3D cube one can quantify the biases in the calculations from the previous section. In particular, we can check whether the transition Eq. 1→ Eq. 2 still gives a valid estimate for a realistic gas density distribution.
To verify this, we took the example case of a ρ CR ∝ r −1 profile and computed the projected cosmic-ray density P CR (r) directly and using the P CR (r) = S (r)/P gas (r) relation. The result of this comparison is shown in Fig. 8, where the exact projected shapes for α = {1, 2} profiles are also given.
It is clear from this figure that the inferred profile in the α = 1 case would resemble a steeper one with α ≈ 1.5. This way an α ≈ 1.0 − 1.5 measurement, presented in Sect. 3.3, is likely to correspond to the true α ≈ 1 distribution -accounting for this bias. In order to directly check this assumption, we have performed a fit of inferred profiles of a form ρ CR ∝ r −α + const for α = [0; 2] to the projected cosmic-ray densities, obtained here (Fig. 4) and earlier in H.E. S.S. et al. (2016). The results of this test, shown in Fig. 9, indeed suggest that the true, deprojected cosmic-ray profile has a slope of α = 0.88 +0.16 −0.07 (at 68% confidence level), consistent with the α = 1 assumption.  Fig. 8. Projected cosmic-ray density as a function of the distance from the GC, computed using a 3D reconstruction of the gas distribution in the region. The underlining profile is assumed to have a ρ CR ∝ r −1 shape. Blue line shows the true projection, whereas the orange one depicts the P CR (r) = S (r)/P gas (r) estimate. Dashed and dotted grey lines represent the exact solutions for the α = {1, 2} cases respectively, arbitrarily scaled to fit the image. The calculation was performed for the ±10 pc Galactic latitude slice along the Galactic plane.
In a similar way we can also estimate the effect of the uncertain line-of-sight measurements on the derived value of α. To reconstruct the effect we randomly shifted the positions of the cube grid cells in the z direction and then regenerated the expected cosmic-ray profiles P CR (r). The random shifts were obtained by adding 10-40 Gaussian-distributed displacements with standard deviations from 50 to 200 pc at random positions within the cube. The amplitude of the Gaussians was fixed to 25 parsecs, while their width puts a coherence scale, not allowing the neighbouring bins to have very different shifts. The resulting shifts in Projected 1/r 2 Fig. 10. Projected cosmic-ray density as a function of the distance from the GC, computed using a 3D reconstruction of the gas distribution in the region. The profiles are inferred from the P CR (r) = S (r)/P gas (r) estimate using the random line-of-sight shifts of the underlying gas distribution (see Sect. 4.2 for details). The true profile assumed in the simulation is ρ CR ≈ r −1 . Dashed and dotted red lines represent the exact solutions for the α = {0.5, 1, 1.5, 2} cases respectively, arbitrarily scaled to fit the image.
the cube do not exceed ≈ ± 50 pc, roughly resembling the corresponding measurement uncertainties.
The outcome of this calculation is shown in Fig. 10. The comparison with the exact α = {1, 2} profiles projections indicates that the uncertainty in the line-of-sight position measurements can dramatically change the shape of the measured cosmic-ray profile. Within our simulation setup, assuming an α = 1 cosmic-ray profile, we can generally conclude that a measurement of α ≈ 1.5 is equally likely.
It should also be noted that the main evidence for the centrally peaked cosmic-ray profile comes from the measurements in the central ≈30 − 50 pc, where the uncertainty on the gas content is the highest. Any unaccounted absorption of the molecular radio emission from such dense regions may lead to a general underestimation of the gas content in the direct vicinity of the GC (in addition to the line-of-sight uncertainty), and thereby to a more modest estimate of the cosmic-ray density there.

Discussion
Based on deep VHE γ-ray observations using the MAGIC telescopes we can reconstruct the cosmic-ray distribution profile in the GC region. In this process we use the simplifications outlined in Sect. 3.3, which result from our partial knowledge of the line-of-sight distribution of the target gas material. This implies that measurements of the cosmic-ray distribution in the GC vicinity can only be improved with more accurate gas observation in other (radio, infra-red etc) domains and are not limited by the sensitivity of current γ-ray instruments. Accordingly, even highly sensitive future CTA observations of this region, for an analysis like the one presented here, depend on higher-resolution molecular line measurements.
Our estimations described in Sec. 4.2 indicate that in the absence of such measurements the derived cosmic-ray profile may look sharper than it is in reality. At the same time, the scatter in the profiles, based on the quoted uncertainties in the line-of-sight gas clouds positions, may make any interpretation of the measured shape less reliable. As an example, the α = 1 cosmic-ray profile, simulated in Sect. 4.2, may appear as α between 0.5 and 1.5 depending on the exact locations of the clouds (see Fig. 10).
The existence of a centrally peaked cosmic-ray distribution around GC, revealed by the H.E.S.S. data (H.E.S.S. et al. 2016) and confirmed here, however, seems reliable. The MAGIC data suggest an α = 1.2 ± 0.3 profile and therefore do not contradict the α = 1 scenario, in which CRs diffuse outwards from a central source, especially given the considerations outlined above.
This peaked cosmic-ray profile, together with a hard emission spectrum, was interpreted as a signature of PeV cosmic ray acceleration by SgrA* (H.E. S.S. et al. 2016). It should be noted here that alternative explanations for the enhanced density of CRs at GC exist, ranging from millisecond pulsars (Guépin et al. 2018) to dark matter (Lacroix et al. 2016). It has also been shown by Gaggero et al. (2017) that it is possible to consistently model the diffuse γ-ray emission from GC in the Fermi and VHE emission energy range with particles from the Galactic CR sea, nearly without the need of an addition central CR source. The authors also state that in the presence of such a CR sea, the maximum energy of any excess from a central source is even less certain.
MAGIC measurements of the diffuse γ-ray spectrum over a larger region of the central molecular zone, ≈150 pc in width, favor, on a ≈2σ level, a cut-off in the γ-ray spectrum over a pure power law. The 1σ confidence range for the cut-off energy spans from 10 TeV to 80 TeV (see Sect. 3.4). This corresponds to proton energies of ≈0.1 − 1 PeV, which means that the data are still marginally compatible with the PeVatron scenario.

Conclusions
The presence and proximity of a central super-massive black hole make the GC region a unique laboratory for studying comicray acceleration near black holes. The above presented MAGIC