© MAGIC Collaboration 2020

1 Introduction

The Crab Nebula broad-band emission is usually interpreted in the framework of leptonic models. The radio to megaelectronvolt (MeV) gamma-ray emission is attributed to synchrotron radiation of energetic electrons in the 120−150 μG nebula magnetic field. At higher energies, giga- to teraelectronvolt (GeV to TeV) emission is linked to the inverse Compton (IC) scattering on the synchrotron, infrared (IR) and cosmic microwave background (CMB) photons. The morphology of the Nebula, revealed by the optical and X-ray data, is non-trivial. Nonetheless, its broad-band spectral energy distribution (SED) is reasonably well described even by one-zone models, involving diverse electron distributions (Meyer et al. 2010) and/or propagation and cooling effects (Martín et al. 2012; Fraschetti & Pohl 2017).

At the same time, most of the proposed models fail to describe the details of the Crab Nebula SED (Aleksić et al. 2015). The highest model-to-data deviations lie in the kiloelectronvolt (keV) to MeV range, where the SED softens, and GeV to TeV range, where the IC peak appears broader than suggested by several models. In addition to IC-related emission, bremsstrahlung and proton–proton interactions may also contribute to the GeV–TeV emission, if emitting electrons are at least partially confined in the filaments of the nebula, filled with the ionised gas (Atoyan & Aharonian 1996, see also Sect. 5).

One way to resolve this degeneracy is offered by observations at the highest energies above several tens of TeV, where at least the bremsstrahlung process gives a subdominant contribution to the Nebula emission. The dominant emission at those energies is due to the combination of the synchrotron self-compton (SSC) and IC/CMB emission of electrons with energies ≳ 10¹³ eV. The SSC part of the emission at these energies is produced in the deep Klein-Nishina regime , with being the energy of the background photons, E_e that of electrons, m_e the electron mass, and c the speed of light. This means that the SSC spectrum traces that of the underlying electron population. This way, the apparent changes in the synchrotron spectrum at keV-MeV energies should also manifest themselves in the ~ 10−100 TeV energy band. The absence of the corresponding spectral changes at these energies would indicate the subdominant nature of the leptonic SSC emission at the highest energies, emitted by the Nebula in favour of other competing mechanisms.

Observations at energies ≳10 TeV are usually associated with low event count rates from astrophysical sources. The collection area A_eff of Imaging Atmospheric Cherenkov Telescopes (IACTs) is determined by the size of the Cherenkov light cone from the γ-ray induced extended air showers (EAS). For vertical observations the collection area of a single telescope is ~ 0.05 km². The collection area can be increased using a larger number of telescopes, like the forthcoming Cherenkov Telescope Array (CTA) observatory. Alternatively, a similar effect can be achieved by using observations at higher zenith angles (ZA; Sommers & Elbert 1987) (see also Konopelko et al. 1999; Ahnen et al. 2017). This observation mode leads to an increase in the Cherenkov pool size due to the larger distances to the showers. At the same time the reduced photon density on the ground shifts the energy threshold of the telescope to significantly higher energies. Technical details of the novel very large zenith angle (VLZA) observation technique can be found in Mirzoyan et al. (2018).

In this paper we present the results of the Crab Nebula observation at VLZAs (> 70°) with the MAGIC telescopes and discuss them in the context of other multi-wavelength data of this source.

2 MAGIC very large zenith angle observations of the Crab Nebula

2.1 The MAGIC telescopes

The Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescopes are a system of two 17 m diameter IACTs, located at an altitude of 2200 m above sea level at the Roque de los Muchachos Observatory on the Canary Island of La Palma, Spain (28°N, 18°W).

The telescopes are used to image flashes of Cherenkov light produced by the charged component of EAS initiated in the upper atmosphere by gamma-ray photons with energies ≳ 30 GeV. Both telescopes are nominally operated together in coincidence mode (also called stereoscopic mode), in which only events simultaneously triggering both telescopes are recorded and analysed (Aleksić et al. 2016a). For low zenith angle (< 30°) 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 h of observations (Aleksić et al. 2016a).

2.2 Observations

The data sample presented here was accumulated from December 2014 to November 2018 in the ZA range 70°− 80° and comprises ≈56 h of good-quality data (after the initial data selection as described below; ≈ 88 h before the selection), taken in the so-called wobble mode (Fomin et al. 1994). The summary of observational time per year is given in Table 1.

Two configurations can be used for VLZA observations of the Crab Nebula by MAGIC: during the source rise and during the set above the horizon. These configurations give somewhat different sensitivities of the MAGIC stereoscopic system to the incoming γ-ray flux due to the varying projected inter-telescope distance, as seen from the direction of the source (stereo baseline). The two MAGIC telescopes are located north-east and south-west from the system centre, thus providing a larger baseline in the north–north-west and south–south-east directions. With the declination of ≈ 22°, the Crab Nebula rises at ≈77° and sets at ≈283° azimuths (counting from the north), i.e. north–north-east and north–north-west correspondingly. Here both “rise” and “set” configurations were used with most of the data taken in the set direction due to the larger stereo baseline.

Observations in the VLZA regime correspond to shower distances of ≳ 50−100 km from the telescopes as opposed to ~10 km at lower zenith angles ≲30°. As such, these measurements are subject to increased light attenuation due to the scattering and absorption in the atmosphere. The standard MAGIC way to account for this effect (i.e. using a dedicated micro-LIDAR system; Fruck et al. 2014) only allows us to probe the atmospheric absorption at distances ≲ 20 km. To ensure an appropriate control over the wavelength-dependent atmosphere attenuation we took additional contemporaneous images of the stellar field next to the Crab Nebula with red, green, and blue filters which allows the total atmospheric transmission to be monitored with an accuracy better than 10%. The details of this procedure are given in Appendix B.5 (see also Mirzoyan et al. 2018).

We did not change the optical focusing of the telescopes (usually set to a distance of 10 km) as our tests with the Monte Carlo simulation did not indicate any significant performance improvement by doing so.

3 Data analysis

3.1 MAGIC data analysis

The acquired data are reduced with the standard MAGIC Analysis and Reconstruction Software (MARS; Zanin et al. 2013). We first remove events detected during adverse weather conditions and those corresponding to the known temporary hardware issues. Due to the nature of VLZA observations, several usual data cuts are no longer efficient. These include the presence of the clouds in the telescope’s field of view measured with an infrared pyrometer system (Gaug et al. 2014), and the number of stars detected by the MAGIC star-guider cameras. The corresponding measurements were used to cross-check the applied event selection. The latter was based on the cuts on the mean currents of photomultipliers, the event trigger rate, and the LIDAR transmission at the maximum accessible range.

We use the standard MAGIC MARS routines to reconstruct the initial direction and impact distance with respect to the MAGIC telescopes for the recorded EAS images. These were augmented with contemporaneous atmospheric transmission monitoring and corrections, as explained in Appendix B.5.

To reconstruct the energy of the EAS initiating particle, three different methods were used: (a) standard MARS procedure based on a look-up table (LUT, Aleksić et al. 2012) created from the MAGIC Monte Carlo (MC) simulations, (b) random forest (RF) multivariate analysis, and (c) neural network (NN) regression. The applied cosmic-ray background suppression was based on the classification scheme implemented with both RF (Albert et al. 2008; Aleksić et al. 2012) and NN. All these techniques yield consistent results. To derive the results shown below, LUT energy estimation and RF event classification techniques were employed.

3.2 Fermi/LAT data analysis

In this work we made use of the publicly available Fermi/LAT Pass 8 data set¹. We used the Fermi Science Tools package² v9r33p0 for data processing, retaining only the “Source” (P8R2_SOURCE_V6) class events, registered until June 2017 within the 75° zenith angle; we did not apply the ROI-based zenith angle cut. The photons selected from the 20° region around the Crab Nebula position were also required to lie within the Crab Pulsar 0.60-0.82 phase range where the nebula emission dominates (Abdo et al. 2010). We did not apply an additional gating of the Crab Nebula flares (Abdo et al. 2011; Buehler et al. 2012). These flares mostly affect the low-energy synchrotron emission of the Nebula, and due to their ≲ 500 MeV cut-off energy are practically undetectable above ~3 GeV energy, and even in the 0.3−1 GeV energy range contributions of the flares to the average Nebula flux in 9 yr is limited to 5− 10% due to their short duration. The fluxes of all the sources in the 300 MeV–510 GeV energy range in the selected region were estimated from the joint likelihood fit. Given the Crab Nebula brightness in the Fermi/LAT energy range, the fitted model included the diffuse background components (namely iso_P8R2_SOURCE_V6_v06.txt and gll_iem_v06.fits) and 20 brightest objects from the 3FGL catalogue (Acero et al. 2015) within 28 degrees from the source of interest. The Crab Nebula model itself comprised a single power law; using more complex multi-component models is not required with the narrow energy bins (five per decade) used here.

Fig. 1 MAGIC collection area, estimated for an observational sample in the zenith angle range 70°− 80° with Monte Carlo simulations (see also Mirzoyan et al. 2018), and for comparison the collection area (for 20° zenith angle observations in so-called Production 3 layout) of the future CTA array (see Sect. 4 for details). Also shown is the MAGIC collection area at lower zenith angles from Aleksić et al. (2016a).

4 Results

During the VLZA data taking, the low-energy threshold after the data selection cuts quickly increases from ~ 1 TeV at zenith angle of 70° to ~ 10 TeV when approaching 80°. The collection area at energies above 70 TeV quickly reaches approximately 2 km² (compared to ~ 0.1 km² for low zenith angle observations, Aleksić et al. 2016a), leading to an unprecedented gamma-ray collection area.

To estimate the MAGIC performance for the acquired VLZA data, we used a dedicated Monte Carlo simulation, describing the MAGIC observations of gamma-ray induced air showers in the zenith angle range 70°− 80°. This simulation was performed with the Corsika code (Heck et al. 1998), modified to include the MAGIC specific output. It also included thecurvature of the Earth’s atmosphere to properly describe the increasing air column density during the near-horizon observations. The rest of the simulation procedure was performed the same way as for the lower zenith angleobservations (e.g. Aleksić et al. 2016a). The resulting collection area estimated after the data selection cuts is shown in Fig. 1. For comparison, shown is the expected collection area of the CTA³, which is currently under construction.

It should be noted, that despite of the increase in the collection area, MAGIC VLZA performance is impacted by the limited reconstruction of the shower parameters, resulting from their remoteness (≳ 50−100 km) and correspondingly smaller image size in the telescope camera. In the case of MAGIC, the measured images sizes are decreased to 3–4 pixels for most of the detected showers. Though the number of excess events in our data sample changes with energy as expected given the collection area, we find that these small images degrade the performance of the cosmic-ray background suppression technique we employ. In addition we also note an approximately twofold degradation of the energy and angular resolution, compared to observations at small zenith angles (Aleksić et al. 2016a). Given that this is a novel technique for IACT observations, we anticipate that a certain improvements can be achieved with dedicated, optimised analysis. Still we find that processing of the VLZA data with the standard MAGIC MARS tools at present allows us to perform the interesting studies shown below.

In 56 hof observations, the Crab Nebula signal at estimated energies above 30 TeV was detected at a ≈ 6.5σ Li & Ma (1983) significance level. Despite the increased energy threshold, the spectrum could be reconstructed down to the energy of ~ 1 TeV.

In order to reconstruct the Crab Nebula SED at energies above 1 TeV, we made use of all three energy estimation methods outlined in Sect. 3.1. We applied the background rejection with both the standard MARS routines and dedicated NN, each time adjusting the cuts so as to maintain 90% of Monte Carlo gamma rays in each energy bin.

The SED of the Crab Nebula up to ~100 TeV, obtained with the LUT energy estimation method, is shown in Fig. 2, which shows that the previous HEGRA (Aharonian et al. 2004) spectrum, produced with about 400 h of data, is within ≲ 20% of the MAGIC results. A comparison of these data with the previous lower energy measurements, including the lower ZA MAGIC observationsfrom Aleksić et al. (2015), is given in Fig. 3. Our data do not support the indications for the ~ 30 TeV high-energy cut-off, suggested earlier (Aharonian et al. 2006).

In Fig. 3 we also fit the Crab Nebula SED above 60 GeV with the log-parabola function , also using the archival MAGIC data (Aleksić et al. 2015) in addition to the VLZA measurements. Accounting for 15% systematics (17% for VLZA) in the MAGIC energy scale and 11% (20% for VLZA) on the flux scale (Aleksić et al. 2015, 2016a), we find this fit to be in poor agreement with the data (χ² ≈ 33 over 14 degrees of freedom), indicating that the log-parabola does not provide a good match to the Crab Nebula spectral shape over the entire 60 GeV–100 TeV energy range. Otherwise, the best-fit parameters are similar to those found in Aleksić et al. (2015); for fixed E₀ =1 TeV we find α = −2.48 ± 0.03, β = −0.23 ± 0.01, and f₀ = (2.95 ± 0.27) × 10⁻²³ ph∕(cm² sec eV) (all uncertainties correspond to 1σ confidence range).

As Fig. 3 shows, the overall shape of the GeV to 100 TeV γ-ray emission of the Crab Nebula can be reasonably described within the framework of existing theoretical models like Meyer et al. (2010), Martín et al. (2012), and Fraschetti & Pohl (2017). At the same time, these models do not reproduce the gradual softening of the inverse Compton (IC) emission peak at multi-TeV energies. A simultaneous fit of Fermi/LAT, and archival and VLZA MAGIC data to the best-fit model curves yields relatively high χ² values: 183.0∕26 d.o.f. for Meyer et al. (2010), 77.5∕26 d.o.f. for Martín et al. (2012), and 140.0∕26 d.o.f. for Fraschetti & Pohl (2017) (as presented in Aleksić et al. 2015). It should be noted though that these high χ² values seem to be dominated by the point-to-point systematics, which is not accounted for here.

Fig. 2 Spectral energy distribution of the Crab Nebula, obtained with MAGIC VLZA observations (red). Results of the spectral unfolding (data points) and of the forward-folding (band) procedures are shown. Dashed colored lines (blue, green, and orange, respectively) denote the archival best fits to the Crab Nebula spectra from Aharonian et al. (2004, 2006), and Aleksić et al. (2015). Data points from Aharonian et al. (2004, 2006), Abeysekara et al. (2019) and Amenomori et al. (2019) are also shown for comparison.

Fig. 3 Crab Nebula spectrum obtained here compared to the lower energy measurements with MAGIC (Aleksić et al. 2015) and Fermi/ LAT. The blue shaded band denotes the 68% confidence flux range, estimated from the fit with the log-parabola function. Solid lines of different colours show several leptonic models from Meyer et al. (2010), Aleksić et al. (2015) and Fraschetti & Pohl (2017), previously constructed for the Crab Nebula. Dashed lines denote predictions for hadronic contribution from Bednarek & Protheroe (1997) (BP), and Amato et al. (2003) (AGB). Shown is the bulk wind Lorentz factor (Γ), and is the effective target material density increase over its mean value.

5 Discussion

The novel method of VLZA observations with the MAGIC telescopes allows us to detect gamma rays up to hundreds of TeV in a few tens of hours. It can be efficiently used to search for astrophysical sources accelerating particles to PeV energies. This observation technique requires simultaneous measurements of the atmospheric transparency, careful studies of systematics, and properly tailored MC data. To this end we have developed an auxiliary atmospheric transmission measurement procedure employing contemporaneous measurements of stellar light from the region next to the γ-ray target.

As discussed earlier, the remoteness of the showers from the telescope during VLZA observations has an impact on the instrument performance. At present this impacts our ability to infer the highest-energy flux from the soft-spectrum sources like Crab Nebula. Along with other differences in the VLZA shower development (Neronov et al. 2016), this indicates that a revision of the analysis technique is needed to fully reveal the full potential of VLZA data taking.

Still, the obtained VLZA Crab Nebula data at a few tens of TeV agree well with the earlier, lower zenith angle measurements with HEGRA (Aharonian et al. 2004) that were obtained in observations that were approximately eight times longer. VLZA measurements presented here also support the source spectrum extension up to 100 TeV and likely beyond. These results are consistent with findings of HAWC (Abeysekara et al. 2019) and Tibet ASγ collaborations (Amenomori et al. 2019).

Generally, the available multi-wavelength Crab Nebula flux measurements can be explained within the framework of leptonic models; however, they do not provide much flexibility in the spectral shape of the inverse Compton emission part (Atoyan & Aharonian 1996). Testing the source flux ratios in the energy bins above 0.3 TeV, we found that overall they are consistent with the leptonic framework expectation, following the simplified argumentation in Atoyan & Aharonian (1996). However, a direct comparison of several models to the combined Fermi/LAT and MAGIC data suggests that they do not provide an adequate description of the Nebula flux in the energy range up to 100 TeV, yielding significant flux residuals in the 1–3 TeV energy range. This is consistent with the Aleksić et al. (2015) conclusions that the IC peak is broader than expected from the modelling.

At the same time, several theoretical studies have suggested that the highest energy emission of the Crab Nebula may have at least a partial contribution from the interaction of accelerated particles with the ambient medium (e.g. Atoyan & Aharonian 1996; Bednarek & Protheroe 1997; Bednarek & Bartosik 2003; Amato et al. 2003). Such interactions may be intensified if particles are trapped in the Nebula filaments, leading to a approximately twentyfold effective target density increase (Atoyan & Aharonian 1996). The extension of the Crab Nebula synchrotron emission to ~ 100 MeV energies implies the presence of PeV electrons. This suggests the hadrons in the Nebula could also be accelerated to similar energies, given their generally lower losses via synchrotron radiation. Energetic protons would naturally contribute to the γ-ray flux at TeV energies via p−p interactions.Several predictions from such models are plotted with dashed and dotted lines in Fig. 3.

As can be seen from Fig. 3, the MAGIC data disfavour a putative significant hadronic contribution to the measured TeV flux. This in turn implies that the accelerated nuclei constitute a minor fraction of the pulsar wind power and/or do not generally interact with the overdense Nebula filaments. Further VLZA observations can be used to refine this statement via a more accurate spectral shape estimation at the highest energies.

Overall, the VLZA observation technique extends MAGIC sensitivity to the highest energies, allowing us to search for galactic PeVatrons in the pre-CTA era. With an appropriate adaptation, this technique may also be used by CTA itself to further boost its sensitivity to the highest energy γ-ray events.

Acknowledgements

We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF and MPG, the Italian INFN and INAF, the Swiss National Fund SNF, the ERDF under the Spanish MINECO (FPA2015-69818-P, FPA2012-36668, FPA2015-68378-P, FPA2015-69210-C6-2-R, FPA2015-69210-C6-4-R, FPA2015-69210-C6-6-R, AYA2015-71042-P, AYA2016-76012-C3-1-P, ESP2015-71662-C2-2-P, FPA2017-90566-REDC), the Indian Department of Atomic Energy, the Japanese JSPS and MEXT, the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-153/28.08.2018 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia “Severo Ochoa” SEV-2016-0588 and SEV-2015-0548, and Unidad de Excelencia “María de Maeztu” MDM-2014-0369, by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project 13.12.1.3.02, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, the Polish National Research Centre grant UMO-2016/22/M/ST9/00382 and by the Brazilian MCTIC, CNPq and FAPERJ. This research has made use of the CTA instrument response functions provided by the CTA Consortium and Observatory, see http://www.cta-observatory.org/science/cta-performance/ (version prod3b-v1) for more details.

Appendix A Artificial neural network tools for data analysis

Currently, the multivariate analysis methods based on artificial intelligence are being extensively developed and widely used. In the ground-based Cherenkov gamma-ray technique, characterised by a very small level of events constituting the source signal, the application of artificial neural networks (ANNs) has shown very good performance (Maneva et al. 2003).

Both energy reconstruction and gamma-hadron classification for VLZA analysis were checked by applying two different neural network tools. The first is based on JETNET Fortran package (Peterson et al. 1994) implemented for ROOT (Brun & Rademakers 1997) and MARS via C++ wrapper. The second uses the modern deep learning Tensor Flow⁴ libraries implemented in the KERAS package (Chollet et al. 2015). In both analyses ANN results were added to the standard for MAGIC ROOT output files, so that the entire analysis program chain of MAGIC could be applied. We used feed-forwardalgorithms with back-propagation minimisation. The network architecture consists of an input layer, three hidden layers, and an output layer. As input we used Hillas parameters for both telescopes as well as the results of stereoscopic reconstruction (e.g. EAS impact distance and height of the shower maximum). The performance of the two tools is similar and comparable to that of the MAGIC standard RF energy estimator. The advantage of Tensor Flow library is that it is several times faster than JETNET code, enabling better optimisation of the network architecture.

Appendix B Systematic uncertainties in MAGIC VLZA measurements

The systematic uncertainties associated with the MAGIC VLZA observations largely overlap with those derived at lower zenith angles (Aleksić et al. 2016b,a) with regard to the telescope hardware performance. At the same time, the increased distance to EAS at large zenith angles affects the performance of the reconstruction techniques applied. Potentially, this makes them more sensitive to the otherwise small discrepancies between the MAGIC MCs and real VLZA EASs. Below we quantify the systematics contributions specific to VLZA data taking.

B.1 Pointing accuracy

Increased gravitational loads during the VLZA observations, caused by the nearly horizontal orientation of the telescopes, result in bending of the telescope structure and camera support arch. This effect is largely compensated by an active mirror control system and contemporaneous observations of positions of a number of reference stars next to the MAGIC field of view (Aleksić et al. 2016b,a).

To evaluate the effect from the remaining telescope mispointing, we follow the approach taken in Aleksić et al. (2016a) and compare the true sky position of the Crab Nebula to that derived from the data on several different epochs of observations. Due to the lower count rates during the VLZA observations, this comparison is not possible on a nightly basis. We thus combine the data in sets spanning one or more months to properly determine the fitted source position in the skymaps. From these maps we determine the mispointing as the difference between the nominal and the reconstructed source positions.

The mispointing data versus time can be found in Fig. B.1. A two-dimensional plot of the mispointings in RA/Dec coordinates can be found in Fig. B.2. We can conclude that the systematic uncertainty on the reconstructed source position is ≈0.04°. This value is larger than the one reported in Aleksić et al. (2016a), which can be explained by the increased mechanical stress on the telescope structure when observing near the horizon in the VLZA regime. We note that this value is ≈ 4 times smaller than the MAGIC point spread function (PSF) during such observations.

Fig. B.1 Time evolution of mispointing of Crab Nebula using the VLZA method. Each data point represents one or more months of taken data. The solid line represents the results from the analysis using the complete Crab VLZA data set, and the dashed lines represent the error on that value.

Fig. B.2 Mispointing in RA/Dec of each datapoint, as in Fig. B.1. The circles represent specific angular distances, as seen in the legend. The data point where all of the Crab VLZA data are analysed is represented in red.

Furthermore, an additional estimate of the remaining systematic uncertainty on the event direction reconstruction can be obtained from the extension of the reconstructed Crab Nebula total angular profile with respect to the Monte Carlo estimated PSF. The MAGIC PSF shape is well described by a King function (Aleksić et al. 2016a; Da Vela et al. 2018), whereas the putative additional mispointing that is random in nature would appear as an additional smearing of this profile. The fit of the VLZA Crab Nebula data above 10 TeV suggests that such additional mispointing does not exceed 0.04° at the 90% confidence level, as shown in Fig. B.3.

Fig. B.3 Scan of an additional extension of the Crab Nebula angular profile above 10 TeV photon energy above the MC predicted PSF. The Y -axis gives the resulting χ2 value of the fit with 16 degrees of freedom. The extension is expressed as the σ of the additional Gaussian component. The maximum allowed extension provides the measure of the remaining MAGIC mispointing.

The additional mispointing yields a broader PSF than predicted by Monte Carlo simulation, thus affecting the true event containmentwithin a given off-source angle cut. The resulting effect depends on the original (energy dependent) PSF width, and thus changes with the energy. Using the MAGIC VLZA simulations we estimate the impact of this mispointing to be ≲ 4% at ~ 3 TeV and ≲ 8% at ~ 30 TeV energies.

B.2 Zenith angle dependence of the instrument response

The rapid growth of the air mass with the zenith angle in the 70°− 80° range results in a gradual change of the MAGIC response with respect to the incoming EAS. To account for this effect we split the MAGIC MC sample into 100 bins in cosine of zenith angle (in the 0°− 90° range; such binning roughly follows the growth of the air mass). The instrument response functions were then computed by re-weightingthis sample with the zenith angle distribution of the accumulated data, as shown in Fig. B.4.

B.3 MC to data comparison

The development of extended air showers at zenith angles above 70° proceeds primarily in the rarefied layers of the upper atmosphere and at distances of ≳ 50−100 km from the observer. These conditions lead to certain peculiarities in the shower evolution (dependent on the nature of the primary particle; see Neronov et al. 2016), which means that VLZA observations may be associated with a larger MC to data discrepancy than that derived from lower zenith angle data (Aleksić et al. 2016a).

In order to verify this, we compared the distribution of the basic EAS Hillas parameters Size, Length, and Width (Hillas 1985) as well as the height of the maximum shower development in MAGIC Monte Carlo and real event samples. For the latter we used the excess distributions of the same parameters in the on- and off-source regions, derived with loose event selection cuts. This comparison is shown in Fig. B.5 for events in the zenith angle range 70–75deg where no significant difference is present between the real and simulated Crab Nebula VLZA signals.

Fig. B.4 Zenith angle distribution of the accumulated VLZA Crab Nebula observational sample, used to re-weight MAGIC response functions.

B.4 Energy bias and resolution

We estimatethe MAGIC energy resolution and bias comparing the true MC simulated event energies to those obtained with our energy reconstruction algorithms. In order to parametrise the accuracy of the reconstruction we fit a Gaussian to the scaled (E_est − E_true)∕E_true distribution of the estimated energies E_est in narrow bins of E_true. The mean of this distribution is taken as a measure of bias, whereas its width corresponds to the energy resolution of the applied reconstruction procedure.

The resulting energy estimate is subject to uncertainties in the overall MAGIC light throughput, which are estimated to be ≲ 15% (Aleksić et al. 2016a). To estimate the possible impact on the VLZA energy reconstruction, we apply an additional scaling of the amount of light in the VLZA MC sample by ±15%. These scaled MCs are then processed as if no light scale had been applied. The bias and resolution resulting from them (as a function of the true event energy) are given in Fig. B.6.

As can be seen, the overall energy bias varies in the range [−20%;+15%], which gives an estimate of the instrumental MAGIC energy scale uncertainty in the VLZA regime. Since the total amount of light recorded from EAS plays a major role in the event energy reconstruction, it is worth noting that the resulting energy bias is almost directly proportional to the assumed light scale.

B.5 Atmospheric transmission

Atmospheric transmission directly affects the amount of EAS light that reaches the telescope camera. The uncertainties on the transmission contribute to the overall uncertainty of the VLZA energy scale, as discussed above.

To estimate the wavelength-dependent atmospheric transmission, we image the stellar field next to Crab Nebula withdedicated CCD cameras, located close to the centres of the MAGIC-I and MAGIC-II reflector dishes. The images were taken every 90 s with the 90-s exposure, cyclically changing the colour filters from red (λ_mean ~ 640 nm) to green (λ_mean ~ 530 nm) to blue (λ_mean ~ 450 nm). The acquired images were flat-fielded and cleaned of hot pixels and dark current contribution. Then counts from selected bright stars were estimated as a difference of counts from the circular region around the star and the background counts from an annular region of a larger diameter.

Fig. B.5 Comparison of basic EAS reconstruction parameters between the MC simulated (blue) and real (red) event samples, recorded in the 70–75 deg zenith angle range. Size, Length, and Width are the so-called Hillas parameters (Hillas 1985), whereas the shower maximum height is reconstructed from the MAGIC stereo data by a standard analysis pipeline in MARS. Deviation of real distributions from MC does not exceed the 2σ confidence level.

Fig. B.6 Left: bias in the VLZA energy estimation computed for different scaling values of the MAGIC light throughput system. It is estimated as the mean of the Gaussian distribution of the estimated MC energies and is given relative to the assumed true energy in the simulation. Right: VLZA relative energy resolution for different scaling values of the MAGIC light throughput system computed as the sigma of the Gaussian distribution of the estimated MC energies.

In order to calibrate this aperture photometry procedure, an additional imaging of this stellar field was performed on several clean nights, when atmosphere absorption showed no variation with height. During such nights light absorption in each colour filter followsa simple law: (B.1)

Here c is the number of background-subtracted CCD counts, c₀ is the number of counts before absorption, α ≈const is the specific absorption coefficient, and m_air is the air mass at a given zenith angle z. The constant c₀ can be determined from Eq. (B.1) from measured CCD counts at different zenith angles from the selected star. Knowing c₀, the average absorption coefficient α during the subsequent observational sessions can be estimated as α = −log(c∕c₀)∕m_air(z).

Contemporaneous imaging of those selected stars during the VLZA data taking allows us to estimate atmospheric transmission for EAS with a temporal resolution of 1.5−3 min. The height of each shower maximum, estimated as a part of the standard data reduction in MARS, is used to compute the line-of-sight distance to the shower maximum and derive the corresponding value of the air mass . The resulting absorption can then be estimated as . The ratio of this latter value to the absorption assumed in the MAGIC detector Monte Carlo simulations (for the same zenith angle and shower distance) defines the relative light scale s = τ_data∕τ_MC in each of the colour filters, which – after a convolution with the Cherenkov emission spectrum – is finally used to correct the estimated event energies or select time intervals with good transmission.

Though the stellar light measurements described above provide a simple and reliable way to estimate the total atmospheric transmission, they are subject to inaccuracies due to the uncertainties in the derived calibration constants c₀ and uncertainties in the measured CCD counts during the observations. We minimised the latter by choosing the camera exposure time so that the reference stars get ≳3 × 10⁴ CCD counts and the resulting uncertainty is less than 1%.

The uncertainty on the calibration constants c₀ was computed from several estimates, taken on nights with stable atmosphere transparency (< 10% deviations from Eq. (B.1) law), as shown in Fig. B.7. The standard deviation of these estimates suggests that the calibration constants for the reference stars are determined with an accuracy ≲ 5%.

Fig. B.7 Recorded CCD counts c from the star HD 35708 as a function of air mass, converted to the stellar magnitudes m = m0 − log2.512c with an arbitrary reference point m0. These counts were obtained with the green filter during the CCD imaging of the Crab Nebula region on nights with stable atmosphere transparency. Solid lines correspond to the fits of each night (see Eq. (B.1)) assuming α = const; their extrapolation to the zero air mass defines the reference calibration constant c0 for the particular filter used. The RMS of these extrapolations sets the c0 uncertainty and is ≈0.04 mag for this calibration sample.

The calculation of the exact transmission correction for a specific air shower from such stellar light measurements depends on the assumeddistribution of the absorber in the atmosphere, which induces an additional uncertainty. Even though for mild (≲ 20%) light absorption the transmission estimates are generally accurate to within 3–5%, our estimates suggest that deviations of up to 10% are still possiblein some cases. We thus conservatively use the 10% value as an estimate of the atmosphere transmission systematics.

B.6 Total systematic uncertainty

To estimate the total systematic uncertainty in the VLZA case, we also account for non-VLZA specific sources of MAGIC systematics,reviewed in Aleksić et al. (2012); the resulting list of contributions is given in Table B.1. Combining them we find that the telescope’s energy scale is determined with an accuracy of 19% at low energies (~ 3 TeV) and 17% at medium energies (~30 TeV). This is comparable to the MAGIC energy scale systematics estimated from the muon analysis, worsened by the larger uncertainty in the atmosphere transmission due to VLZA conditions. The uncertainty on the reconstructed flux normalisation (excluding the effect of the energy scale) is 14 and 20%, respectively.

References

All Tables

All Figures

Fig. 1 MAGIC collection area, estimated for an observational sample in the zenith angle range 70°− 80° with Monte Carlo simulations (see also Mirzoyan et al. 2018), and for comparison the collection area (for 20° zenith angle observations in so-called Production 3 layout) of the future CTA array (see Sect. 4 for details). Also shown is the MAGIC collection area at lower zenith angles from Aleksić et al. (2016a).

In the text

Fig. 2 Spectral energy distribution of the Crab Nebula, obtained with MAGIC VLZA observations (red). Results of the spectral unfolding (data points) and of the forward-folding (band) procedures are shown. Dashed colored lines (blue, green, and orange, respectively) denote the archival best fits to the Crab Nebula spectra from Aharonian et al. (2004, 2006), and Aleksić et al. (2015). Data points from Aharonian et al. (2004, 2006), Abeysekara et al. (2019) and Amenomori et al. (2019) are also shown for comparison.

In the text

Fig. 3 Crab Nebula spectrum obtained here compared to the lower energy measurements with MAGIC (Aleksić et al. 2015) and Fermi/ LAT. The blue shaded band denotes the 68% confidence flux range, estimated from the fit with the log-parabola function. Solid lines of different colours show several leptonic models from Meyer et al. (2010), Aleksić et al. (2015) and Fraschetti & Pohl (2017), previously constructed for the Crab Nebula. Dashed lines denote predictions for hadronic contribution from Bednarek & Protheroe (1997) (BP), and Amato et al. (2003) (AGB). Shown is the bulk wind Lorentz factor (Γ), and is the effective target material density increase over its mean value.

In the text

	Fig. B.1 Time evolution of mispointing of Crab Nebula using the VLZA method. Each data point represents one or more months of taken data. The solid line represents the results from the analysis using the complete Crab VLZA data set, and the dashed lines represent the error on that value.
In the text

	Fig. B.2 Mispointing in RA/Dec of each datapoint, as in Fig. B.1. The circles represent specific angular distances, as seen in the legend. The data point where all of the Crab VLZA data are analysed is represented in red.
In the text

	Fig. B.3 Scan of an additional extension of the Crab Nebula angular profile above 10 TeV photon energy above the MC predicted PSF. The Y -axis gives the resulting χ2 value of the fit with 16 degrees of freedom. The extension is expressed as the σ of the additional Gaussian component. The maximum allowed extension provides the measure of the remaining MAGIC mispointing.
In the text

	Fig. B.4 Zenith angle distribution of the accumulated VLZA Crab Nebula observational sample, used to re-weight MAGIC response functions.
In the text

Fig. B.5 Comparison of basic EAS reconstruction parameters between the MC simulated (blue) and real (red) event samples, recorded in the 70–75 deg zenith angle range. Size, Length, and Width are the so-called Hillas parameters (Hillas 1985), whereas the shower maximum height is reconstructed from the MAGIC stereo data by a standard analysis pipeline in MARS. Deviation of real distributions from MC does not exceed the 2σ confidence level.

In the text

Fig. B.6 Left: bias in the VLZA energy estimation computed for different scaling values of the MAGIC light throughput system. It is estimated as the mean of the Gaussian distribution of the estimated MC energies and is given relative to the assumed true energy in the simulation. Right: VLZA relative energy resolution for different scaling values of the MAGIC light throughput system computed as the sigma of the Gaussian distribution of the estimated MC energies.

In the text

Fig. B.7 Recorded CCD counts c from the star HD 35708 as a function of air mass, converted to the stellar magnitudes m = m0 − log2.512c with an arbitrary reference point m0. These counts were obtained with the green filter during the CCD imaging of the Crab Nebula region on nights with stable atmosphere transparency. Solid lines correspond to the fits of each night (see Eq. (B.1)) assuming α = const; their extrapolation to the zero air mass defines the reference calibration constant c0 for the particular filter used. The RMS of these extrapolations sets the c0 uncertainty and is ≈0.04 mag for this calibration sample.

In the text