MAGIC very large zenith angle observations of the Crab Nebula up to 100 TeV

Aims. We measure the Crab Nebula γ -ray spectral energy distribution in the ∼ 100 TeV energy domain and test the validity of existing leptonic emission models at these high energies. Methods. We used the novel very large zenith angle observations with the MAGIC telescope system to increase the collection area above 10 TeV. We also developed an auxiliary procedure of monitoring atmospheric transmission in order to assure proper calibration of the accumulated data. This employs recording optical images of the stellar ﬁeld next to the source position, which provides a better than 10% accuracy for the transmission measurements. Results. We demonstrate that MAGIC very large zenith angle observations yield a collection area larger than a square kilometer. In only ∼ 56 h of observations, we detect the γ -ray emission from the Crab Nebula up to 100 TeV, thus providing the highest energy measurement of this source to date with Imaging Atmospheric Cherenkov Telescopes. Comparing accumulated and archival MAGIC and Fermi /LAT data with some of the existing emission models, we ﬁnd that none of them provides an accurate description of the 1 GeV to 100 TeV γ -ray signal.


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 A&A 635, A158 (2020) 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 13 eV. The SSC part of the emission at these energies is produced in the deep Klein-Nishina regime E γ bkg E e /(m 2 e c 4 ) 1, with E γ bkg 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 2 . 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. Notes. The observational time is given separately for the rise and set of the source on the horizon. Observation time is given in hours, after the data selection cuts.

MAGIC very large zenith angle observations of the Crab Nebula
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).

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.
A158, page 2 of 10 MAGIC Collaboration: MAGIC VLZA observations of the Crab Nebula up to 100 TeV

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 , 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.

Fermi/LAT data analysis
In this work we made use of the publicly available Fermi/LAT Pass 8 data set 1 . We used the Fermi Science Tools package 2 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 lowenergy 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.

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 2 (compared to ∼0.1 km 2 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 the curvature 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 angle observations (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 3 , 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 10 1 10 2 Energy, TeV  Dashed colored lines (blue, green, and orange, respectively) denote the archival best fits to the Crab Nebula spectra from Aharonian et al. (2004Aharonian et al. ( , 2006, and Aleksić et al. (2015). Data points from Aharonian et al. (2004Aharonian et al. ( , 2006, Abeysekara et al. (2019) and Amenomori et al. (2019) are also shown for comparison. 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 h of 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 observations from 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 dN/dE = f 0 (E/E 0 ) (α+β log 10 (E/E 0 )) , 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(Aleksić et al. , 2016a, we find this fit to be in poor agreement with the data (χ 2 ≈ 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 0 = 1 TeV we find α = −2.48 ± 0.03, β = −0.23 ± 0.01, and f 0 = (2.95 ± 0.27) × 10 −23 ph/(cm 2 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 10 3 10 2 10 1 10 0 10 1 10 2 Energy, TeV It should be noted though that these high χ 2 values seem to be dominated by the point-to-point systematics, which is not accounted for here.

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 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.  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. A158, page 7 of 10 A&A 635, A158 (2020) 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.
The additional mispointing yields a broader PSF than predicted by Monte Carlo simulation, thus affecting the true event containment within 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-weighting this 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 onand 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-75 deg where no significant difference is present between the real and simulated Crab Nebula VLZA signals.

B.4. Energy bias and resolution
We estimate the 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 with dedicated 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.
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 follows a simple law: Here c is the number of background-subtracted CCD counts, c 0 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 0 can be determined from Eq. (B.1) from measured CCD counts at different zenith angles from the selected star. Knowing c 0 , the average absorption coefficient α during the subsequent observational sessions can be estimated as α = − log (c/c 0 )/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 m EAS air . The resulting absorption can then be estimated as τ data = exp (−αm EAS air (z)). 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 0 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 4 CCD counts and the resulting uncertainty is less than 1%.
The uncertainty on the calibration constants c 0 was computed from several c 0 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%.
The calculation of the exact transmission correction for a specific air shower from such stellar light measurements depends on the assumed distribution of the absorber in the atmosphere, which induces an additional uncertainty. Even though for mild A158, page 9 of 10 A&A 635, A158 (2020) Recorded CCD counts c from the star HD 35708 as a function of air mass, converted to the stellar magnitudes m = m 0 − log 2.512 c with an arbitrary reference point m 0 . 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 c 0 for the particular filter used. The RMS of these extrapolations sets the c 0 uncertainty and is ≈0.04 mag for this calibration sample.
( 20%) light absorption the transmission estimates are generally accurate to within 3-5%, our estimates suggest that deviations of up to 10% are still possible in 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  Aleksić et al. (2012). The values affecting the telescope energy scale and flux normalisation are labelled "ES" and "FN" correspondingly.
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.