© The Authors 2025

1 Introduction

Observations of very-high-energy (VHE) cosmic gamma-ray sources are crucial for understanding the extreme processes that accelerate particles to PeV energies in our Galaxy (e.g., Gabici et al. 2019; de Oña Wilhelmi et al. 2024) and beyond (Bose et al. 2022b). However, detecting VHE gamma rays presents significant challenges. Since photon fluxes from cosmic sources follow negative power laws, which are also reduced by propagation in the interstellar medium, satellite-based experiments such as Fermi-Large Area Telescope (Atwood et al. 2009) do not have enough collection area in their detection time frame to measure them. Ground-based experiments, such as imaging atmospheric Cherenkov telescopes (IACTs) and water Cherenkov arrays, offer a much larger collection areas than space-borne detectors by exploiting the showers of secondary particles produced when gamma-ray photons interact with atmospheric nuclei. These secondary particles can be detected over a large area on the ground and discriminated by the cosmic ray-induced background. In brief, IACTs detect the Cherenkov light emitted shortly after charged relativistic secondary particles in the showers pass through a medium. These are captured by the optical systems of IACTs and focused on their nanosecond-sensitive cameras. Alternatively, other experiments detect the light induced by shower charged particles directly in the water of their tanks (see e.g., Fegan 1997; de Naurois & Mazin 2015; Bose et al. 2022a; Sitarek 2022, for a review of gamma-ray detection techniques). Beyond TeV energies, the flux sensitivity of IACTs is also limited by the number of telescopes spread on a given area and by their optical field of view (FoV), due to the large angular size of the high-energy showers at moderate impact distances. Additionally, VHE gamma-ray sources in the Galaxy are often extended (see e.g., Cao et al. 2024; H.E.S.S. Collaboration 2018), further limiting their observability with the current generation of IACTs, which tend to have a rapidly decreasing gamma-ray acceptance with increasing offset angle on a scale of ~1°(e.g., Aleksić et al. 2016).

The Single-Mirror Small-Size Telescope (SST-1M) was developed to probe VHE gamma rays while addressing the usual IACT limitations, such as a small FoV and limited duty cycle due to moonlight. Its optical system features a 4-m diameter spherical mirror dish based on the Davies-Cotton design (Davies & Cotton 1957), maintaining a good image quality over a large FoV, while minimizing optical aberrations (Alispach et al. 2025). That is crucial for observations of extended sources, followups of poorly localized transients¹ and a good reconstruction of VHE gamma-ray showers generating images that can extend over several degrees in the camera. The aperture of the Cherenkov camera is composed of a 3 mm thick Borofloat window, which integrates a narrow-band optical filter composed of dielectric layers and an anti-reflective coating (Alispach et al. 2025). This enables the transmission of Cherenkov light to the pixel sensors while reducing the NSB. The camera comprises 1296 silicon photomultipliers (SiPMs) (Heller et al. 2017), enabling safe operation in the presence of high night sky background (NSB) without human intervention, effectively improving robustness against high light rates and duty cycle compared to cameras equipped with photomultiplier pixel (Anderhub et al. 2013). The degradation of the SiPM performances under high NSB has been studied in (Nagai et al. 2019) for the SiPM of the camera and thus are copped for in the analysis. Moreover, SiPMs are insensitive to the magnetic field and no high-voltage power supply is needed.

Currently, two SST-1M telescope prototypes are operating at the Ondřejov Observatory in the Czech Republic, at an altitude of 510 m a.s.l. The telescopes² are separated by 152.5 meters, and their timestamps are synchronized to nanosecond precision using the White Rabbit timing network (Serrano et al. 2009). They operate in a stereoscopic mode, allowing for precise shower reconstruction when observed from different directions. Both telescopes are currently in the commissioning phase, focusing on optimizing operations and gaining a detailed understanding of all components of the stereoscopic system. A comprehensive description of their design and operation is provided in Alispach et al. (2025).

The Crab Nebula, is a nearby (~2 kpc) remnant of a supernova observed in 1054 A.D., hosting a pulsar with a rotation period of 33 ms. It was the first detected VHE gamma-ray source (Weekes et al. 1989), and exhibits a spectrum extending up to PeV energies (LHAASO Collaboration 2021). While signs of variability or flares were reported at GeV energies (Tavani et al. 2011; Abdo et al. 2011), the Crab Nebula maintains a stable flux at VHE, despite years of monitoring with various instruments (e.g., H.E.S.S. Collaboration 2014). It is therefore usually referred to as the “standard candle” of VHE astronomy and is used to assess the performance of VHE gamma-ray observatories. In this study, we evaluate the mono and stereo performance of a pair of the SST-1M telescopes through observations of the Crab Nebula conducted between September 2023 and March 2024.

In Section 2, we describe the calibration procedure and the Monte Carlo (MC) model of the instrument. Section 3 details the MC simulations. In Section 4, we explain the event reconstruction and the data analysis pipeline. Section 5 is dedicated to the mono and stereo performance of SST-1M. In Section 6, we present the Crab Nebula data sample, focusing on the validation of the MC model and presenting the results on the spectrum and the skymap. Section 8 provides a detailed discussion of systematic uncertainties. Finally, we conclude the paper with a summary and final remarks in Section 9.

2 Calibration and MC model of the telescopes

IACTs collect Cherenkov light produced by secondary charged particles in atmospheric air showers. As already explained in Section 1, by imaging this light in the cameras of IACTs, the direction and energy of the primary particle can be reconstructed. The data reduction process, from the photosensor signals into calibrated images, requires a precise characterization of several factors: optical efficiency, baseline, dark count rate, night sky background, gain, and crosstalk. These parameters are used to construct a telescope model for MC simulations, which, together with the precise model of the atmosphere at the observational site, simulate air shower development and Cherenkov photon detection (Section 3). These MC simulations provide a necessary dataset to train machine learning (ML) algorithms used to estimate the energy and direction of incoming particles and discriminate gamma rays from hadrons (Section 4). The MC simulation is obtained using CORSIKA (Heck et al. 1998) for detailed simulation of extensive air showers initiated by gamma rays and cosmic rays, and sim_telarray (Bernlöhr 2008), which simulates the ray tracing, geometry and response of the detector, as well as the attenuation of Cherenkov light in the atmosphere.

2.1 Calibrations of the SiPM response

The telescopes are remotely operated via a graphical user interface, which allows users to execute calibration procedures and monitor the telescope subsystems during operation. The dark run is performed before the observation after the thermal stabilization of the telescope camera, with the camera lid closed. Another dark run is performed after the end of the observation. The data are used to obtain pixel gain, dark count rate, optical crosstalk, and electronic noise. Analogue signals from the pixels are routed to the crate-based electronics called Digicam, where they are digitized using 250 MHz Flash Analog-to-Digital Converters. The calibration process, based on Alispach et al. (2020) and further developed, provides the necessary conversion coefficients from analog-to-digital converter (ADC) counts to unbiased photoelectrons (p.e.). Within this procedure, the baseline, namely, the average measured signal in 1024 preceding time bins is computed in the Digicam and subtracted.

2.1.1 Photoelectron spectrum

We defined the gain, g [ADC/p.e.], of a pixel, coupled with its readout electronics, as the average number of integrated ADC counts recorded per single avalanche event in the SiPM. During observations, optical crosstalk may occur when a secondary avalanche is triggered within a secondary SiPM micro-cell by photons emitted during the discharge of a neighboring primary cell (Nagai et al. 2019). This effect occurs only between microcells within the same SiPM pixel and, if it is not corrected, leads to an overestimation of the measured p.e. count. The conversion factor from ADC counts to p.e. is given by $g^{\ast }=\frac{1}{\nu g}$ [p.e./ADC], where ν is the averaged p.e. produced per primary avalanche. These parameters are routinely monitored before and after each night of observation through dedicated dark runs as shown in Alispach et al. (2025). During dark runs, spontaneous avalanches are triggered by thermal electrons with a Poisson probability, P_λ. The total number of produced p.e. per primary avalanche, accounting for crosstalk, is well described by a Borel distribution of parameter ν^* (Vinogradov 2012). These events exhibit the same characteristics as photon-triggered avalanches, allowing dark runs to be used for evaluating the SiPM’s response. For each pixel, the p.e. spectrum is constructed by stacking the maximum charge from the integrated signal in n_s = 7 consecutive samples from each waveform. The ADC distribution of dark events is built from approximately 30’000 randomly triggered waveforms of n_T = 50 samples. The resulting histogram is then fitted using the probability density function of the ADC counts (n_ADC), namely, the p.e. spectrum: pes(n_ADC), built by summing each n p.e. component as $$\mathrm{p}\mathrm{e}\mathrm{s}(n_{\mathrm{A}\mathrm{D}\mathrm{C}})=P_{\lambda }(0)A_{\mathrm{e}\mathrm{l}}+(1-P_{\lambda }(0))A_{\mathrm{p}\mathrm{e}},$$(1)

where A_el describe the ADC distribution for waveforms without dark events described by a Gaussian distribution, G, with respective mean μ_el and standard deviation ${\sigma }_{\mathrm{e}\mathrm{l}}:A_{\mathrm{e}\mathrm{l}}=\mathrm{G}({\mu }_{\mathrm{e}\mathrm{l}},{\sigma }_{\mathrm{e}\mathrm{l}})(n_{\mathrm{A}\mathrm{D}\mathrm{C}})$. A_pe represents the sum of each n > 0 p.e., weighted with the point mass function, following the Borel distribution, $$B_{{\nu }^{\ast },n}=\frac{e^{-n{\nu }^{\ast }}(n{\nu }^{\ast }{)}^{n-1}}{n!}.$$(2)

Then, $$A_{\mathrm{p}\mathrm{e}}={\displaystyle \sum _{n=1}^{\mathrm{\infty }}}{B}_{{\nu }^{\ast },n}\mathrm{G}(n_0+ng,\sqrt{{\sigma }_{\mathrm{e}\mathrm{l}}^2+n{\sigma }_{\mathrm{p}\mathrm{e}}^2})(n_{\mathrm{A}\mathrm{D}\mathrm{C}}),$$(3)

where σ_pe represent the standard deviation of ADC counts distribution per single avalanche charge. The averaged p.e. produced per primary avalanche, ν, used in the analysis is corrected for the rate of dark event in the integration window: $\nu =\frac{1}{1-{\nu }^{\ast }}-\lambda \frac{n_s}{n_T}$. $n_0=-\nu \lambda g\frac{n_T}{n_s}$ is the true level of the electronic baseline.

This fitting procedure gives an estimation of the SiPM gain and crosstalk probability with relevant uncertainties for each camera pixel. In Figure 1, we show the averaged distribution of dark ADC counts for all camera pixels.

To reliably describe the response of the SiPM camera to incoming photons, we tuned the sim_telarray MC model parameters to comply with the measured dark run data. In the sim_telarray, the response parameters of individual pixels are drawn from a single population with parameters defined by the MC model. The tuning of the model includes the optical crosstalk probability, electronic noise, and the amplitude of a p.e. pulse in the waveform, which is the parameter used by sim_telarray to transform the detected p.e. into ADC counts. To mimic the thermal rate of avalanches (dark noise) in SiPMs, an NSB value in the simulations was set to 1.8 and 1.5 MHz for SST-1M-1 and SST-1M-2, respectively. The single p.e. pulse shape was determined from the data, stacking more than 50’000 individual p.e. pulses collected during dark runs. To incorporate the above-described crosstalk in simulations, we use the singlephotoelectron spectrum generated from a Poisson branching process for the optical crosstalk smeared by the gain fluctuation of the SiPM. The resulting distribution is a (smeared) Borel distribution. The model have been found to be in good agreement with the experimental probability distributions for dark counts of SiPMs (Vinogradov 2012).

A comparison of the measured and MC-simulated dark runs is shown in Fig. 1, where we present a normalized distribution of waveform-extracted charges of all camera pixels across 11’972 and 100’000 events for the data and simulations, respectively. The first peak, centered slightly above 0 ADC counts, corresponds to a complex interplay between the baseline subtraction, search for the pulse maximum in the noise-only waveform and the electronic noise. The slight disagreement in its shape for SST-1M-2 between the first and second peak can be attributed to the asymmetric distribution of the electronic noise, more important in SST-1M-2. This effect is currently not simulated and will be the subject of further studies. Nevertheless, the high-level analysis is not affected as the measured signal is composed of stacked pulses coming from the second peak and further peaks, which have been aptly modeled for both telescopes.

Fig. 1 Comparison between data and MC of the normalized (to area) ADC charge distributions extracted from the dark runs for the two telescopes. To avoid overlaps of the two curves, the SST-1M-2 distribution is multiplied by 0.1. Each count is calculated from the 50 samples (200 ns) long waveform in the same manner as done for the triggered data. See the explanation in the text. An asymmetric distribution of the electronic noise (more significant in SST-1M-2) can explain the slight disagreement between the first and the second peak.

Fig. 2 Left: relation between muon ring radius and ADC counts under various NSB conditions for SST-1M-1. Right: behavior of measured muon ring charge as a function of the baseline shift induced by the NSB. Blue dotted and orange dashed lines represent the expected behavior given in Nagai et al. (2019) for SST-1M-1 and SST-1M-2, respectively. Filled blue (SST-1M-1) and empty orange (SST-1M-2) markers represent the fitted estimation of the relative ADC counts for 1.15° radius muon images.

2.1.2 Voltage drop induced by NSB

The SST-1M cameras are equipped with SiPM sensors, which offer several advantages over traditional photomultiplier tubes (PMTs). Unlike PMTs, SiPMs can operate under high and continuous illumination without experiencing aging effects or damage. For IACTs, continuous light primarily originates from the NSB, which includes diffuse starlight, artificial light sources, moonlight, and ground reflections. Continuous exposure to light causes a steady current to flow through the SiPM bias circuitry. To prevent this high current from causing damage, a bias resistor (R_bias) is connected in series with the SiPM. The current flowing through R_bias induces a drop of the biasing voltage, which translates naturally into a reduction of the gain (g), the photon detection efficiency (PDE), and optical crosstalk probability. The frontend electronics of the SST-1M is DC coupled so the direct current resulting from the feedback loop can be evaluated by measuring the pixel’s baseline shift. The behavior of SiPM parameters under such conditions has been fully described and experimentally tested for the camera of SST-1M-1 in Nagai et al. (2019).

The two camera prototypes in operation feature different R_bias (R_bias = 10 kΩ for the camera of SST-1M-1 and R_bias = 2.4 kΩ for SST-1M-2) as it was observed after operation of SST-1M-1 that lowering R_bias would decrease the impact of their order to lower the effects of NSB for the second prototype. This results in significantly different behaviors of the two telescopes under the same NSB levels. The baseline shifts induced by the NSB are monitored during observations using 100 Hz clocked interleaved pedestal events. The impact of the voltage drop induced by the NSB on the reconstructed charge can be evaluated using Cherenkov light produced by cosmic muons under different NSB conditions. Figure 2 shows the behavior of the average ADC counts in cleaned muon images with the baseline shift for both cameras (see Section 2.2 for a detailed description of the muon analysis), demonstrating a good agreement between the voltage drop impact on the muon image charge data and the expected electronic behavior that was modeled according to Nagai et al. (2019). This good agreement confirmed that the model can be applied in the data processing, during the calibration step, to correct the image intensity (sum of the integrated pixel charges) as a function of the voltage drop. The voltage drop is derived by measuring the shift between the baselines measured during the dark runs and the one estimated every second during the observation. The extracted intensity is then corrected according to the model to a level that would be the one measured in the absence of NSB.

2.2 Optical efficiency calibration

The optical efficiency of an IACT refers to the effective fraction of Cherenkov light transmitted through the entire optical system before reaching the photon sensors. It is influenced by mirror reflectance, the NSB filtering window transmittance, the reflectance of the lightguides coupled to the SiPMs, the shadowing of the mirror by the telescope structure, and the SiPM PDE (Heller et al. 2017). Since the early 1990s (Fleury & Artemis-Whipple Collaboration 1991), Cherenkov light produced by atmospheric muons has been recognized as an excellent test beam for the estimation of optical efficiency (Gaug et al. 2019). Atmospheric muons are produced in hadronic showers at high altitudes and can be detected if they pass near one of the telescopes. Cherenkov light is emitted along an emission cone following the muon trajectory. If this trajectory aligns with the telescope FoV, the Cherenkov cone of light will sweep large portion of the mirror and be projected on the focal plane as a ring, with a radius that corresponds to the Cherenkov angle. The characteristic shape and consistency of these images make them an ideal tool for calibration and testing the accuracy of the MC model in the simulations.

Muon events are collected during routine observations and processed using a dedicated analysis pipeline. Based on the small timing spread in photon arrivals and stable delay between the trigger and the signal, charge extractions are done using a unique fixed integration window of 28 ns for all pixels. The integrated images are cleaned using the tailcuts algorithm (Lessard et al. 2001). The brightest pixels in the camera can be identified when they are above the so-called image threshold, set at 5 p.e.. Clusters of pixels, called islands, are built by adding neighboring pixels if their intensity is above the boundary threshold set at 4 p.e. Here, islands with two or more neighboring pixels were retained. A circle was fitted to the surviving pixels. We estimated the muon charge by summing the charge in all pixels with centers within 0.15 degrees of the fitted circle. Taking into account both the pixel geometry and the optical PSF, this led to the collection of approximately 90% of the muon signal. Only events with a fitted circle radius between 0.8 and 1.4 degrees were selected for the analysis. To assess the completeness of the muon ring, the signal within the ring region was summed across 12 regions, each covering a 30° angular section. Only images that exhibited a signal above 7 p.e. in at least 10 out of the 12 regions were considered for further analysis. Additionally, the signal outside the selected pixels was also taken into account. Images with an integrated charge exceeding 20 p.e. outside the selected pixels were excluded. These criteria have been designed to filter out shower images and noise fluctuations, ensuring that only muon images with well-defined ring structures are used in the calibration studies. During typical observations, the muon image rates after selection are approximately 0.1 Hz and 0.2 Hz for SST-1M-1 and SST-1M-2, respectively.

A dedicated MC simulation of muon events with parameters listed in Section 3 was used to tune the optical efficiency of the instrument model. Figure 3 shows the relation between the muon ring radius and its integrated intensity. The measured values for both telescopes in the time period of the Crab Nebula observation campaign (see Section 6) are compared with the one obtained with MC simulations, showing the consistency between observational data and the instrument model. The optical efficiency of the telescopes is not constant over time. Gradual degradation of mirror reflectivity, camera window transmittance, and the performance of camera electronic components contribute to a slow but steady decline in the reconstructed muon intensity. Figure 4 shows the evolution of the optical efficiency with time for both telescopes. During the time period of the Crab Nebula observation campaign presented in this work, the degradation of optical efficiency is estimated to be approximately 2%.

2.3 Optical point spread function

The optical point spread function (PSF) is measured with a dedicated CCD camera installed in the center of the mirror dish and a dedicated PSF screen installed on the camera lid (Alispach et al. 2025). In sim_telarray, the PSF can be reproduced by smearing the horizontal and vertical alignment of individual mirror facets (whose optical properties are known from lab measurements), along with ray-tracing of light of a simulated point-like source. The optical PSF was measured for both telescopes in December 2023, providing D₈₀ = 10.4 mm and 12.2 mm for SST-1M-1 and SST-1M-2, respectively, where D₈₀ represents the diameter of a circle containing 80% of the reflected light by the mirrors. This translates to the on-axis PSF of approximately 0.1°. The stiffness of the telescope structure ensures that the PSF does not depend on the telescope pointing direction and is stable in time. The PSF is measured and optimized by pointing at a bright star and adjusting individual mirror facets to reach the smallest PSF in the PSF screen after any optical system intervention, and periodically twice a year (Alispach et al. 2025). The current observation did not show any PSF degradation over time.

Fig. 3 Charge of muon images as a function of the fitted ring radius for SST-1M-1 and SST-1M-2. The solid markers represent observational data collected from September 2023 to April 2024, during the first stereo Crab observation campaign (blue for SST-1M-1 and orange for SST-1M-2). Measured data agree with the MC simulation, shown with open markers and dashed lines. Notable differences between SST-1M-1 and SST-1M-2 arise from differences in their optical designs.

Fig. 4 Evolution with time of the charge in the muon ring pixels for fixed radius of r = 1.2° for SST-1M-1 (blue points) and SST-1M-2 (orange points). Each data point corresponds to one single observation night. The green-shaded area indicates the time period of the Crab Nebula observation campaign presented in this work. The linear fits (dashed for SST-1M-1, dotted for SST-1M-2) enable the estimation of the relative decrease in optical efficiency over time.

2.4 Model of the atmosphere in Ondřejov

A model of the molecular atmosphere, namely, its density profile as a function of altitude, is crucial for a reliable description of particle interactions in air showers, production of Cherenkov light, and its subsequent attenuation along the path to the telescope. For the Ondřejov Observatory, it was extracted from the ECMWF ERA 5 database³. Moreover, precise estimation of the amount of aerosols present in the atmosphere is important when determining the primary particle’s energy, since the incoming Cherenkov light gets attenuated by Mie scattering. The aerosols also contribute to the overall NSB, since the scattering aerosol particles act as a light source. Unfortunately, there is no instrument capable of high-precision measurements of the vertical aerosol optical depth (VAOD) at the Ondřejov Observatory. As such, we made use of the open access results from the Sun photometer located in Košetice⁴ (Blumthaler et al. 1997). This site is located 45 km away from Ondrrejov, but otherwise bears remarkable similarities. It is located at the same height above sea level, in the midst of a hilly landscape and shielded from the nearest village, at a sufficient distance from the nearest town.

The instrument measures the VAOD at a range of wavelengths. We make use of the filter band at 380 nm. This is justified by the transmissivity of the camera and of the optical system being rather flat in the relevant spectral region, but the Cherenkov spectrum increasing sharply at lower wavelengths. The Sun photometer provides only daytime measurements, so for each nighttime measurement, we calculate the VAOD as the mean of the average values of the preceding and following day. It turns out that the VAOD is remarkably stable for this particular dataset of Crab sample measurements and oscillates around the value of 0.05. We estimated that the systematic uncertainties in this sample related to VAOD estimation are approximately of this order, and so we use this value for MC simulations for the whole period. It is worth noting that the observed VAOD is remarkably low; for example, a similar period in 2024 exhibits large variations with a mean value of approximately 0.2. It should be mentioned, however, that the sun photometer measures an integral VAOD, while the shower reconstruction requires a proper measurement of the height-differential VAOD and this must be properly addressed in the atmospheric calibration strategy of IACTs (e.g., Gaug 2017). Such a simplification can lead to systematics on the energy scale, as discussed in Section 8.

With this VAOD value in hand we then calculate the atmosphere transmissivity using the MODTRAN software (Berk et al. 2014) in 50 layers spaced between 0.5 and 100 km, utilizing an average seasonal vertical atmospheric density profile from the ECMWF ERA 5 dataset. It can then serve as an input to the sim_telarray simulation.

3 MC simulations

The MC simulations were produced in two steps. First, particle interactions and Cherenkov light from air showers in CORSIKA v7.7402 (Heck et al. 1998). Second, the attenuation of the light and the response of the SST-1M telescopes in sim_telarray v2021-12-25 (Bernlöhr 2008). To cover different zenith angles at which the Crab Nebula had been observed, five fixed zenith angle regions of 20°, 30°, 40°, 50°, and 60° were simulated. In each of the regions, three sets were simulated, namely the diffuse samples of gamma rays and protons to establish the shower reconstruction, and a point-like sample of gamma rays to directly calculate the SST-1M response to a point-like source. To keep the number of simulations manageable, we chose to simulate for a fixed azimuth at which the sources culminate, that is, the telescopes were looking to the South in all simulations. Incoming directions of primaries for the diffuse samples were randomly thrown around the telescope-axis direction up to 10° opening angle and sampled on the ground in a way to emulate an isotropic flux. The maximum distance up to which the impact points were sampled in the telescope-tilted ground frame was set between 928m - 1273 m for gammas and 1032 m - 1414 m for protons, both evolving with the zenith angle θ as cos^−0.5 θ motivated by Abe et al. (2023b) and confirmed by us to cover the full detection volume. Similarly, for calculation efficiency purposes, the lowest simulated energy evolved with the zenith angle from 200 GeV to 2 TeV for gammas and from 400 GeV to 3 TeV for protons to roughly catch the changing energy threshold of the detection. The maximum simulated energy was 800 TeV for gammas and 1300 TeV for protons to probe the SST-1M sensitivity at the highest energies, and the simulated energy spectrum was $\frac{\mathrm{d}\mathrm{N}}{\mathrm{d}\mathrm{E}}\propto E^{-2}$.

While the SST-1M camera is designed to operate under a wide range of NSB conditions, and the telescopes can also collect data under the full moon, the results presented in the following analysis are limited to dark conditions. The gain, PDE, and optical crosstalk were set following the model described in Nagai et al. (2019). The MC events were corrected in the same way as the real data and used to train the energy and angular regression models, as well as the gamma-hadron classification one. For typical dark nights at the Ondřejov Observatory, the average NSB rate is approximately 100 MHz per pixel. Due to slightly different hardware of the cameras, see Section 2.1.2, the rate of random pulses in pixels caused by NSB was set differently for SST-1M-1 and SST-1M-2 to 94 MHz and 120 MHz, respectively. The left panel of Figure 5 shows the distribution of the averaged standard deviation of the baselines, which are directly impacted by the NSB level, compared to the values used in the MC simulations for both cameras. To stay conservative, the NSB levels in the MC were selected to describe the reconstructed p.e. distributions from pedestal events under the NSB slightly larger than typical during the Crab Nebula observations, depicted in the right part of Figure 5. Despite that the Ondrrejov Observatory is located in an area with a strict public lighting policy, occasional car flashes or the light from nearby villages and towns lead to a substantial NSB variability and a strong dependence on the zenith angle. To estimate the systematic uncertainty of the telescope sensitivities connected to the discrepancy between the NSB level in the data and MC (Section 8) we also performed simulations under the lower NSB conditions, representing the darkest nights of observation, with NSB rates of 55 MHz and 72 MHz for SST-1M-1 and SST-1M-2, respectively.

At the zenith angle of 20°, we simulated 33 millions diffuse gamma-ray showers and 3 millions point-like gammas, each resampled in CORSIKA 20 times, and 26 millions protons, each resampled 100 times. The amount was reduced for other zenith angles to keep roughly the same number of events above the minimum simulated energy. In total, 1.86 billion of diffuse gammas, 260 millions of point-like gammas and 7.4 billions of proton induced events were available for the telescope simulation in sim_telarray and further processing in sst1mpipe (see Section 4).

Besides the simulations necessary for the reconstruction setup and Instrument Response Functions (IRFs) determination, a dedicated set of muon simulations was performed to assess the overall optical throughput of the telescopes described in Section 2.2. It consists of 5 millions of muons directed to the telescope with the random scatter around the telescope axis direction up to 6°, impact distances up to 3 m, and the energy range of 4 GeV - 1 TeV following the power-law spectrum $\frac{\mathrm{d}\mathrm{N}}{\mathrm{d}\mathrm{E}}\propto E^{-2}$.

The full MC production was split into two independent samples for machine learning-based reconstruction (using random forests, RF) of the energy, direction, and classification - the training sample used for RF training, and the testing sample for the estimation of the performance and production of IRFs. In order not to limit the RF performance with the training sample size, we conducted a detailed study of mono and stereo performance, varying the number of training events. We carefully monitored energy, angular resolution, and gamma-hadron separation power, especially for high energies, which are most affected by the lack of statistics. It turned out that the best performance can be reached if the number of mono diffuse gammas and protons in the last arbitrary energy bin (500-800 TeV) is higher than ≈15 000, which we kept roughly constant for all zenith angles (events which survived cleaning; see Sect. 4).

Fig. 5 Left : distribution of the averaged standard deviation of the baselines for SST-1M-1 (dashed line) and SST-1M-2 (solid line). The vertical lines represent the value assumed in the MC model of SST-1M-1 (94 MHz/pixel) and SST-1M-2 (120 MHz/pixel). The vertical axis represents the number of hours with a baseline standard deviation (on the horizontal axis) in the final Crab Nebula data sample. Right: comparison between the MC-predicted and pedestal events-measured p.e. distributions under the slightly higher than typical NSB conditions during the Crab Nebula observations.

4 Event reconstruction

The data and MC shower reconstruction and analysis are performed in sst1mpipe⁵ (Juryšek et al. 2025), which is a standalone Python-based pipeline based on ctapipe (Kosack et al. 2021) functionalities. The basic structure of sst1mpipe is inspired by lstchain (Lopez-Coto et al. 2022) and closely follows the ctapipe data model. For low-level event handling sst1mpipe relies on some functionalities of digicampipe⁶, a pipeline developed for camera commissioning, calibration, and low-level data processing. High-level data product of sst1mpipe (photon lists together with IRFs) follows the modern “gamma astro data formats” (GADF) standards (Nigro et al. 2021), making it compatible with contemporary tools for spectral or spatial data analysis like gammapy (Donath et al. 2021).

4.1 Charge extraction

Raw data recorded by each telescope contains for each triggered event a 12-bit digitized waveform in all 1296 camera pixels with 50 samples and 4 ns intervals. First, the baseline, determined as an average of 1024 samples before the readout, is subtracted. The baseline-subtracted waveforms are then calibrated, namely converted into the number of p.e. using the camera gain described in Section 2. Except for this last step, all the following steps are common for data and MC, this last containing already calibrated waveforms, so that both are processed in the same way. Each calibrated waveform is integrated with an eight-sample window around the sample with the maximum signal using the LocalPeakWindowSum algorithm of ctapipe. An integration correction based on the pulse template is applied to correct the remaining charge outside the integration window. In addition to the integrated charge, the pulse time, defined as the amplitude-weighted average of the time samples within the window, is extracted. Faulty pixels, for which the determination of the camera gain failed (typically ~5/20 out of 1296 pixels for SST-1M-1/2), are flagged and the integrated charge, together with the pulse time are replaced by the average value of their neighboring pixels using NeighborAverage method of ctapipe. Constant monitoring of the pixels is done during observation's and faulty detected pixels are treated the same. We note that in the future, the dead pixels can be corrected at the trigger level, exploiting the fully digital readout and trigger system of the SST-1M camera (Heller et al. 2017).

4.2 Image cleaning and parametrization

To remove noisy pixels containing only the NSB, we adopted the two-stage modification of the standard tailcut image cleaning (Abe et al. 2023b). The default tailcuts were set to (8 p.e., 4 p.e.) for both telescopes. These tail cuts were determined empirically as a trade-off between the fraction of true Cherenkov pixels that survived cleaning on one side, and too many noisy pixels left in the cleaned image on the other. First, in the MC sample and the data we optimized for the reference NSB level by demanding about 1% survival rate of pedestal events (which contain only noise). Then we considered the high variability of NSB conditions in Ondřejov, and set the tailcuts higher (by 1 p.e. to 1.5 p.e.) to prevent a higher NSB from degrading the overall performance.

To further suppress noisy pixels due to stars and planets in the FoV or to a significant difference between the NSB in the MC and the data, we increased the default tailcuts in each pixel to max(8p.e., 〈Q_PED〉 + 2.5σ_Qped), where 〈Q_PED〉 is a mean charge in each pixel over 1000 interleaved pedestal events corresponding to 10 seconds, and σ_Qped stands for its standard deviation (Abe et al. 2023b). The multiplicative factor on standard deviation was set to keep the number of altered pixel thresholds at the level of 1-2% so that its impact on the MC and data agreement is negligible. Further improvements of the cleaning performance was achieved by demanding that the arrival time of the neighbor pixel is within 8 ns of the arrival time of the pixel, which exploits the natural correlation of Cherenkov pulse times. In the next step, cleaned pixels containing signal from the Cherenkov photons were parameterized with an extended set of Hillas parameters (Hillas 1985) to describe the shape, orientation, and position of the shower image (for a full list of image parameters used for the event reconstruction, see Figure 6).

4.3 RF reconstruction

4.3.1 Mono reconstruction

To reconstruct the properties of the primary gamma-ray photon and the shower geometry, the image parameters are fed into RFs (Breiman 2001) trained on the training set of MC events using scikit-learn framework (Pedregosa et al. 2011). The RF classifier is used for gamma-hadron separation (resulting in so-called “gammaness”⁷ parameter for each reconstructed event), and the RF regressor is used for energy reconstruction. To determine the arrival direction of each shower in mono analysis, we adopt the so-called DISP method (Lessard et al. 2001), assuming that the source lies on the main axis of the shower image. An RF regressor is trained to reconstruct the distance of the source from the image centroid (disp norm), together with an RF classifier to determine on which side of the centroid along the image axis the source lies (disp sign).

4.3.2 Stereo events

The Crab Nebula observations were conducted with both telescopes triggered independently. Therefore, in the stereo reconstruction of such a data sample, the first step is to match events resulting from the same shower seen stereoscopically with both telescopes. With the typical mono trigger rate of around 200 Hz and White Rabbit sampling precision of 16 ns, the event matching can be based only on event timestamps with negligible probability of random coincidences. Moreover, event matching is performed on the subset of events that survived cleaning, further suppressing the likelihood of random coincidences (below 0.01% of events). In the subsequent analysis, we consider the event stereoscopic, if, for a given shower seen with SST-1M-1, there is a counterpart seen with SST-1M-2 within 10 μs window. Another possibility is exploiting coincidence event tagging by the Software Array Trigger (Alispach et al. 2025), which is currently being tested.

In the reconstruction of stereo events, we proceed similarly to Abe et al. (2023a). First, the geometrical reconstruction of the shower impact distance (tel_impact_distance) and the height of shower maximum (h_max) is calculated using the axis crossing method (Hofmann et al. 1999). The stereo parameters are then used as additional features for the RF training. The energy regressor and gamma-hadron classifier are trained independently for each telescope and their outputs are averaged using the image intensity as a weighting factor. For the shower direction, the stereo reconstruction follows the approach of the MAGIC analysis and reconstruction software (Aleksic et al. 2016). First, disp norm was reconstructed for both telescopes with the RF regressor. Then, four tentative source positions are calculated, averaging the reconstructed positions for both telescopes for all possible combinations of the head-tail shower orientations. Finally, the combination providing the closest angular distance between the reconstructed coordinates from the two telescopes is adopted.

4.3.3 RF performance

Hyperparameters of the RFs were set to reach a trade-off between a good performance and manageable size and training time for each RF. The RF regressors used for energy and disp norm reconstruction use 150 estimators, maximum tree depth 30, minimal number of samples to split a node 10, and the squared error criterion to measure the quality of each split. The gammahadron classifier, and disp sign classifier in the case of the mono reconstruction, employs 100 estimators, a maximum depth of 100, and 10 samples needed to split a node. The Gini impurity is adopted as a measure of the quality of each split, but only the square root of the total number of features is randomly selected to evaluate the Gini index in each step.

The training parameters and their relative importance⁸ for each model and both mono and stereo reconstructions are shown in Figure 6. All RFs, including those for stereo reconstruction, are trained per-telescope, but we only show feature importance for SST-1M-1, as the difference between the telescopes is very small. The most important of the RF mono and stereo parameters of the gamma-hadron classifiers is the Hillas’ ellipse width, which is related to the higher lateral momentum of secondary particles in the hadronic showers. The energy reconstruction is surprisingly dominated by the length of the ellipse, which is due to its correlation with its intensity and impact parameter. We tested the importance of the parameters when the length is not included between the RF parameters, so that the impact and intensity became the most relevant parameters, while the performance degraded by a few percent. We also note that the training sample is dominated by low-energy showers close to the energy threshold with relatively small images, where the correlation of Hillas intensity and Hillas length is even more profound. Interestingly, the stereoscopic parameters (h_max and impact_dist, see Fig. 6) are not taking over the importance for stereo reconstruction, which may suggest that there is room for further performance improvement in stereo. In the case of direction reconstruction, quite expectedly, the Hillas length and time gradient play an important role in mono reconstruction, while in the stereo reconstruction, the impact point distance to the telescope is the most relevant parameter.

The performance of RFs for hadronic background suppression can be assessed with the receiver operating characteristic (ROC) curve. For each gammaness cut, a position along the ROC curve represents the fraction of true gammas and false gammas (i.e., protons mistakenly reconstructed as gammas) left in the sample. Maximizing the former and minimizing the latter, the performance of the gamma-hadron classifier can be expressed as the integral of the ROC curve - area under the curve (AUC). Figure 7 shows the ROC and AUC for both mono and stereo observation regimes at two different zenith angles, with 30° being the typical zenith angle of the Crab Nebula observation from Ondřejov, and 50° representing the high zenith angle performance for comparison. Including stereo parameters in the training features and a combination of gammaness from both telescopes improves the performance of the classifier as expected. It can be noted that the degradation of the performance with increasing zenith angle is due to the increasing telescope(s)-shower distance, leading to smaller shower images and lower precision of the reconstructed shower geometry.

Fig. 6 Gini importance for different RFs for the SST-1M-1. The features used for the reconstruction: log of the intensity, width, length, width/length ratio (wl), timing slope, skewness, kurtosis, leakage, and coordinates of the shower center of gravity in the FoV (x,y). Height of the shower maximum (h max) and distance of the impact point from the telescope (impact dist) are used only in the stereo reconstruction. Upper left: gamma-hadron classifier. Upper right: energy regressor. Bottom left: disp norm regressor. Bottom right: disp sign classifier (Only used in mono reconstruction).

5 Mono and stereo performance at 510 m a.s.l

The RFs trained on the training MC dataset were used to reconstruct the direction, the energy, and to classify the particle of the point-like and diffuse gammas and diffuse protons at all simulated zenith angles. The resulting reconstructed events were then used to calculate the IRFs (effective area, PSF, energy migration matrix, and background model), applying the same quality and selection cuts as in the real data analysis. The IRFs were computed from the MC testing sample using the pyirf Python package (Dominik et al. 2023; Linhoff et al. 2024). A global cut on the image “intensity” of >45 p.e. was applied to bring the data and MC on the same intensity threshold (see Section 6). A cut on the “leakage2”⁹ <0.7 was also introduced to remove significantly incomplete shower images from the sample. The gamma-hadron separation performance based on parametrized shower images depends on the energy of the primary gamma ray. Low-energy gamma ray and proton showers produce less Cherenkov light and so smaller images, so that the gamma-hadron discrimination capability is challenged by the limited resolution of the optical system and by the pixel size of the camera. Similarly to the approach in Abe et al. (2023b), we introduced an energy-dependent gammaness cut to keep the same fraction (i.e., 60%) of gammas in the final sample in all energy bins.

Fig. 7 ROC curve of the RF gamma-hadron classifier for 30° (solid line) and 50° (dotted line) zenith angles. Line colors represent different regimes of observations.

5.1 Energy threshold and effective area

In IACT, the energy threshold (E_T) is the true energy of primary gamma rays, for which the differential event rate of gamma rays with a given spectral energy distribution reaches its maximum. The left part of Figure 8 shows the rate of point-like gamma rays at 20° zenith angle, weighted on the Crab Nebula spectrum (Aleksic et al. 2016) for different stages of the analysis. As a consequence of a low altitude of the Ondřejov observatory, compared to the other IACT sites, the trigger threshold of SST-1M is relatively high. This is because the amount of emitted Cherenkov light depends on the energy of the primary gamma ray, while the amount of collected light also depends on its attenuation between the emission point and the telescopes. Thus, with decreasing altitude of the telescopes, the distance to the shower maximum at typically 10 km a.s.l. increases, and the trigger efficiency drops. At the trigger level, for mono observation E_T ≈ 0.6 TeV, and it raises to E_T ≈ 1.0 TeV at the analysis level, where the energy-dependent cut on gammaness and global cut on θ are applied. Also, θ represents the angular distance between the true and reconstructed gamma-ray directions. The low-energy events are naturally more affected by the quality cuts on image intensity. The energy threshold is higher in the stereo regime, namely E_T ≈ 0.7 TeV and E_T ≈ 1.3 TeV at the trigger and analysis level, respectively, as the low-energy showers emitting fewer Cherenkov photons have a lower probability of triggering both telescopes for significantly different impact parameters. The total trigger rate is also lower in stereo because the probability that the same shower triggers both telescopes is lower than just one telescope being triggered. The energy threshold also strongly depends on the zenith angle of observation as the distance between the shower and the telescopes increases with the zenith angle. For 60° zenith, the analysislevel energy threshold is 6 TeV and 10 TeV for mono and stereo, respectively.

The effective area is the ratio of reconstructed point-like gamma rays over the number of simulated ones, multiplied by the area on which the simulated showers have been thrown in the frame orthogonal to the telescope pointing. It is shown in bins of true energy in the right panel of Figure 8.

5.2 Energy and angular resolution

The energy resolution is defined as 68% containment around the median value of ΔE/E_True = (E_Reco - E_True)/E_True for reconstructed gammas, while the median itself is the so-called energy bias. Both are shown in Figure 9 as a function of E_True for mono and stereo at two different zenith angles. As expected, the energy resolution in Figure 9 improves with energy, reaching Δ E/E_True,mono ≈ 15% and ΔE/E_True,stereo ≈ 10%. At very high energies, the event reconstruction starts to be limited by the containment of the shower images within the camera FoV, degrading the performance. High zenith angle observations, increasing the average shower-to-telescope distance, partially compensate for this effect at the expense of a higher energy threshold. Stereo reconstruction improves the energy resolution significantly at energies below 10 TeV, where the impact parameter is not well reconstructed in the mono analysis.

The angular resolution is defined as the 68% quantile of the distribution of θ. As pointed out by Abe et al. (2023b), in mono reconstruction, the PSF consists of two components: the central component is made of events with properly reconstructed shower orientation (the disp sign parameter). It is surrounded by a ring of events with a wrongly reconstructed disp sign, which becomes the dominant fraction of the PSF near the energy threshold. This feature is not present in stereo by the nature of the stereoscopic reconstruction. The radius of the ring is about 1.5° , which means that these events do not affect the capability of the instrument to resolve two point-like sources and thus we consider only the sample of properly reconstructed gamma rays to characterize the angular resolution in Figure 10. We note that the IRFs used for skymaps and spectral analysis contain the full PSF. The angular resolution in mono is below −0.20° above the energy threshold for both telescopes. In stereo reconstruction, the angular resolution reaches −0.07°. Figure 10 shows the angular resolution for different zenith angles. The effect of the increasing energy threshold with the zenith angle is clearly visible. At very high energies, on the other hand, the increasing zenith and thus increasing telescope-to-shower distance helps to lower the fraction of events not fully contained in the camera, which leads to better angular resolution.

Fig. 8 Left: differential rate of point-like gamma rays with Crab Nebula spectrum in mono and stereo for 20° zenith angle. Event rates for different stages of the analysis are shown: all triggered events (solid line), all events that survived cleaning (dashed line), events that survived quality cuts (dash-dotted line), and events that survived θ2 and gammaness cuts used in the analysis (dotted line). Right: effective area for 20° zenith angle.

Fig. 9 Energy resolution (left) and energy bias (right) for 30° (solid line) and 50° (dotted line) zenith angles evaluated using on-axis point-like gammas for testing. Line colors represent different regimes of observations. For the sake of clarity, only results for SST-1M-1 are shown in mono as both telescopes are very similar.

5.3 Flux sensitivity

The flux sensitivity is defined as the minimum flux of a source that can be detected in a given time, usually 50 hours in the IACT community. We followed the usual definition of the differential sensitivity (e.g., Abe et al. 2023b,a), where the detection with 5σ statistical significance (Li & Ma 1983) should be achieved in individual energy bins dividing each energy decade into five equal logarithmic intervals, assuming the ratio between the signal and background region sizes, namely, the alpha parameter (Li & Ma 1983), of 0.2. We also required at least ten excess events surviving cuts in each energy bin and the signal to background ratio of at least 5%. The former condition affects mostly the high-energy sensitivity bins, where the number of detected gamma rays is typically low due to the naturally decreasing spectrum of all gamma-ray sources, while the latter is the most important at low energies due to low gamma-hadron separation performance, leading to a low signal-to-background ratio.

The differential sensitivity of SST-1M mono/stereo for different zenith angles shown in Figure 11 was evaluated for on-axis point-like gammas re-weighted using the Crab Nebula spectrum (Aleksic et al. 2016), and diffuse protons were re-weighted on protons + Helium spectrum (Alemanno et al. 2024). To select the region of interest, we applied a global θ cut, keeping −80% of reconstructed point-like gamma-rays. Both cuts on gamma effi-ciency¹⁰ and θ were optimized on MC, to reach the best detection significance for a source with the Crab Nebula spectrum. The stereo observation improves the sensitivity by a factor of −2. It should be noted, however, that the sensitivity at low energies reflects the current triggering scheme of the telescopes, and can be further improved using a shorter coincidence window of the acceptance of the two mono-triggers. This would further suppress NSB and enable lower energy events to be included than is currently possible. As expected, the energy threshold increases at high zenith angles as the layer of the atmosphere is larger, requiring higher energy events to generate detectable showers. On the other hand, at energies above ≈50 TeV, the sensitivity improves with the zenith angle due to the effect of the Cherenkov cone geometry, leading to an increased effective area.

Fig. 10 Angular resolution for 30° (solid line) and 50° (dotted line) zenith angles evaluated on testing on-axis point-like gammas. Line colors represent different regimes of observations. For the sake of clarity, only results for SST-1M-1 are shown in mono, as both telescopes are very similar.

Fig. 11 Differential sensitivity at 5σ confidence level for 30° (solid line) and 50° (dotted line) zenith angles evaluated on testing on-axis point-like gammas and diffuse protons. Line colors represent different regimes of observations. For the sake of clarity, only results for SST-1M-1 are shown in mono, as both telescopes are very similar. We note that in contrast to stereo, the additional condition requiring at least ten excess events does not affect the last energy bins for the 50° zenith angle in mono, as the number of excess events is satisfied by the 5σ condition only (see Section 5.3 for more details).

Fig. 12 Simulated gamma-ray rate from the Crab Nebula and simulated protons weighted by the CR spectrum in the signal region after analysis cuts as a function of the offset at 20° zenith angle. Only events with energies above the ET are taken into account.

5.4 Off-axis performance

Observations of the SST-1M telescopes are taken in the wobble mode (Fomin et al. 1994), where the telescope pointing has a small offset to the true position of the source (see Section 6). However, in the case of observation of extended sources, the offset may be larger to avoid contamination of the signal region with the background. Moreover, in the case of follow-up of a poorly localized alert (e.g., gamma-ray bursts, gravitational waves, or neutrino events), the source may end up located anywhere in the telescope FoV. We therefore study the performance of the SST-1M mono and stereo for different source angular distances from the camera center, using MC simulation of point-like gammas and diffuse protons at 20° zenith angle.

Figure 12 shows the integral rate of simulated gammas and protons above the energy threshold (re-weighted on the Crab Nebula, and CR spectrum, respectively) in the signal region after the analysis cuts for both telescopes. We note that in this analysis, we apply an energy-dependent gamma efficiency cut of 60%, optimized for each offset. Due to the large camera optical FoV (9°) and having all pixels part of the trigger, the acceptance is almost flat up to 2.5°, where it drops by ≈10%, which makes the SST-1M an ideal instrument for observation of extended gamma-ray sources.

The integral sensitivity shown in Figure 13 is given in a fraction of the Crab Nebula flux (Crab units, C.U.). It shows a faster degradation than the acceptance above ~3°, as the background acceptance drops slower than the acceptance for pointlike gammas. The angular resolution integrated over all energies above the energy threshold (Fig. 13, right) slowly degrades with increasing offset, reaching about 5% and 11% degradation at 2.5°, for mono and stereo, respectively.

6 Crab Nebula observation

The data sample used in this study was collected between September 2023 and March 2024. The White Rabbit timing network was installed in October 2023, and all data collected before this date are mono only. Due to the maintenance campaign between September 2023 and January 2024, the September 2023 data were all taken with SST-1M-1 only, while a part of the January 2024 data were taken only with SST-1M-2. The observations were performed in wobble mode. It allows for an estimation of the rate of the background events from the regions with the same offset as the signal region, reflected with respect to the center of the FoV. Utilizing the wobble observation strategy, the background can be determined directly from the same dataset as the signal, and no dedicated observations of the background are needed. Being in the commissioning phase, the observation strategy was not yet optimized during the Crab Nebula data acquisition, and thus part of the data sample was taken with 0.7° wobble offset, while all data taken after December 2023 were taken with 1.4° offset to take advantage of having more background regions and therefore better statistics of the background events.

Each observation run, in which telescopes were tracking fixed equatorial coordinates α and δ (one wobble offset), was typically 20 minutes long. In total 46 hours and 52 hours of data were collected with SST-1M-1 and SST-1M-2, respectively. After event matching, we end up with 33 hours of stereo data in total. The zenith angle in the dataset ranges from 27° up to 53°.

Fig. 13 Integral sensitivity for 50 hours of observation (left), and angular resolution for a source with the Crab Nebula spectrum (right) as a function of the offset at 20° zenith angle. All events with energies above the ET are integrated.

Fig. 14 Left: run-averaged rate of mono events with intensity above 200 p.e., corrected for cos(zenith angle). Center: run-averaged fraction of pedestal events that survived the tailcut cleaning. Right: run-averaged fraction of pixels with increased tailcuts. The solid blue line marks the cuts used in the run selection.

6.1 Run selection

At first, we identified a high-quality dataset taken under good atmospheric conditions, also taking into account the highly variable NSB level at Ondřejov. As atmospheric transparency modulates the event rate, the main selection criterion is the stability of the average event rate per run (after the cleaning). The event rate is also affected by the trigger threshold, which can be changed during the observation to account for NSB variability (Alispach et al. 2025). To avoid this affecting the run selection, we require the event rate to be stable for events with intensity above 200 p.e., which is far enough from the intensity threshold to be unaffected by the trigger threshold in low to moderate NSB conditions. The event rate also typically decreases with increasing zenith angle of the source, which may bias the selection. This has been accounted for by correcting the event rate by a factor of 1/ cos(zenith angle), which yields a good parameterisation of the airmass as a function of the zenith angle over the range of zenith angles in our data sample. Figure 14 (left panel) shows the airmass-corrected, run-averaged event rate for both telescopes. To ensure stable atmospheric conditions over the entire data set, we require the airmass-corrected event rate to be >10 Hz in mono and >9 Hz in stereo to account for a slightly lower rate of stereoscopic events of lower intensities. We also rejected very short runs with a live time <200 s, where the statistics of the cleaned events is insufficient to reliably estimate the averaged event rates.

With increasing NSB level, the capability of cleaning to remove noisy pixels is degraded. This degradation can be observed by monitoring the fraction of pedestal events that survive the cleaning cuts (center of Figure 14), for which we require a survival probability lower than 2% in the final sample. Local NSB fluctuations, for instance, due to bright stars in the FoV, are handled with the adaptive cleaning, which raises the tailcuts locally in the affected pixels. To avoid deviating too much from the MC sample with fixed homogeneous NSB across the full FoV, we require that the average fraction of affected pixels in the entire run to stay below 10%. We also checked for a stricter cut on the fraction of affected pixels of 1%, resulting in the ratio of the Hillas parameter distributions (for shower images after the quality cuts) staying within (0.95, 1.05).

Figure 15 shows the run-averaged distributions of event intensities before and after the run selection compared with MC simulations. Above the intensity threshold (which we define as the intensity where the differential distribution reaches the maximum), the differential event rates of selected runs agree very well. Below the intensity threshold, the differences due to variable trigger settings and fluctuations of the NSB spoil the MC-data agreement. Therefore, in further analysis, we introduce a common cut on the intensity >45 p.e., bringing both telescope data and MC to a common analysis threshold.

After the selection cuts, the final sample used for the analysis is 33 hours and 27 hours for SST-1M-1 and SST-1M-2, respectively. The sample of stereo data after the quality cuts is 25 hours. Due to the negligible dead time of the detector (Heller et al. 2017), we consider the final amount of collected hours equal to the effective observation time.

Fig. 15 Per-run rates of the event intensities (i.e., only events which survived cleaning contribute to the distributions) for SST-1M-1 mono (left) and SST-1M-2 mono (center), and SST-1M-1 stereo (right). The distributions shown are corrected for cos(zenith). Dotted lines mark all runs taken, solid lines mark runs after selection cuts. The dashed line shows the distribution of MC protons re-weighted on the proton + helium spectrum. The solid blue line marks the intensity threshold of 45 p.e. introduced in the higher-level analysis.

6.2 Crab Nebula spectrum

To assess the capabilities of the SST-1M telescopes to reproduce the spectral energy distribution (SED) of a known source, we performed the spectral analysis on the run-selected Crab data set in gammapy v1.0.1 (Donath et al. 2021). The analysis was performed independently for stereoscopic events, as well as for all mono events detected with both telescopes, to evaluate the consistency of the obtained results. The cuts applied on the dataset were exactly those used for performance estimation (Section 5).

The signal region size (cut on θ) was set to the value providing the best significance of the Crab Nebula detection (optimized on MC) - 0.2°, and 0.12° for mono and stereo, respectively, which correspond to about 80% containment of the point-like gammaray events. We note that the angular resolution of the instrument is rather stable in the range of zenith angles analyzed (except for the effect of the energy threshold). As the dataset consists of data collected using the wobble method, we adopt the reflected region method for background estimation, where the background is estimated from regions with the same offsets from the center of the FoV as the signal region, assuming the radial symmetry of the acceptance (Berge et al. 2007). As the position of the signal region, we adopted the Crab Nebula coordinates determined by H.E.S.S. Collaboration (2020): α₂₀₀₀ = 83.62875°, δ₂₀₀₀ = +22.01236°.

We performed forward-folding likelihood fit in the energy range 2.5-50 TeV, assuming a Power-Law (PL) spectral shape for the differential flux: dφ/dE = φ₀(E/E₀)^-Γ, where the reference energy E₀ = 7 TeV was fixed on the value close to the decorrelation energy. The lower bound of the energy range was selected correspondingly to the stereo energy threshold below 40° zenith angle, in which most data were taken. We note that the energy threshold for mono is lower (≈1.5 TeV), but aiming at comparing the results of both observation modes, we want to avoid any bias, given by, for instance, a curvature of the intrinsic Crab Nebula spectrum. To justify the use of a single PL spectral model, we tested for a possible curvature in our dataset assuming an alternative Log-parabola (LP) spectral shape in a form of dφ/dE = φ₀(E/E₀)^{-α-β logE/E₀}, where β represents the curvature. As PL and LP are nested models (LP with β = 0 degenerates to PL), we can use likelihood ratio test to determine whether LP is preferred over PL (Wilks 1938; Protassov et al. 2002), resulting in −2Δ log L < 3 and σ < 2 for both mono and stereo. We note that on the relatively narrow energy range given by the high energy threshold on one side, and relatively short integration time on the other, the curvature is expected to be relatively small. The best-fit PL model is shown in Figure 16, and the bestfitting parameters are listed in Table 1. The model parameters are consistent for all three datasets within the uncertainties. We also performed a maximum likelihood estimation of the Crab Nebula flux using the gammapy function FluxPointsEstimator in 15 logarithmically spaced energy bins in the entire fitting range. In this procedure, gammapy recalculates the flux normalization in each energy bin, assuming a PL model with spectral index fixed on the best-fitting value for the entire dataset.

Fig. 16 Left: SED of the Crab Nebula measured with SST-1M telescopes. Different color bands represent the best-fitting spectra for the SST-1M-1 mono (blue), SST-1M-2 mono (yellow), and SST-1M stereo (green) datasets. Right: stereo flux points derived from the PL spectral model. SEDs from other experiments are shown for comparison in the energy range of their validity (left), and corresponding flux points (right) - MAGIC (Aleksić et al. 2015; MAGIC Collaboration 2020), H.E.S.S. (Aharonian et al. 2024), VERITAS (Meagher & VERITAS Collaboration 2015), LST-1 (Abe et al. 2023b), LHAASO (LHAASO Collaboration 2021), and HAWC (Abeysekara et al. 2019).

6.3 Crab Nebula skymap

The excess and significance distributions of the Crab Nebula observations were computed using the ring background method (Pühlhofer et al. 2003) implemented in the gammapy framework (Donath et al. 2021). A background ring with a radius of 1° and a width of 0.3° was used to estimate background events, 0.3° region around the Crab Nebula was excluded from the background estimation. The radial behavior of background events was estimated using MC simulations of diffuse protons. We fitted a two dimensional, symmetrical Gaussian on the excess distribution to estimate the position of the excess. The best-fit coordinates were found to be ((α₂₀₀₀ = 83.62°, δ₂₀₀₀ = 21.99°) ± 0.02°). This position is consistent with the expected Crab Nebula location and in line with the pointing precision of the instrument measured to be 0.02° (Alispach et al. 2025). The excess distribution, convolved with a disk kernel of 0.05° is shown in Figure 17.

The significance map shown in Figure 17 was calculated using the Li and Ma formula for a wide region and convolved with a disk kernel of 0.11°. The Crab Nebula is masked on the figure to highlight the remarkable homogeneity across the entire 7.5° × 7.5° wide field of view. The one-dimensional distribution of the significance map is shown in Figure 17 and compares the significance distribution across all bins and the OFF regions. The OFF distribution is closely described by a Gaussian with the mean μ = −0.16 and standard deviation σ = 1.02, consistent with theoretical expectations for an unbiased background. This demonstrates the accuracy of our hadronic background modeling using MC simulations and paves the way for observations of bright, large, extended sources.

7 MC model validation

The final sample of the Crab Nebula mono and stereo data was used to validate the MC model of both telescopes. A good agreement between the cosmic-ray-dominated event rate and MC protons (re-weighted on p+He spectrum to account for the CR composition) shown in Figure 15 demonstrates that the MC model provides a valid low-level description of the instruments, including calibration, cleaning, and the atmospheric effects discussed in Section 2.4.

To further demonstrate a good agreement at the analysis level, we used the gamma-like events in the signal region in the data and compared them with the point-like MC simulations, reweighted on the number of events from the Crab Nebula (Aleksic et al. 2016) detected in the effective observation time. The data sample is restricted to the zenith angles between 25° and 35° (representing a majority of the dataset), and simulations with the zenith angle fixed at 30° . We further suppressed the background in the signal region by subtracting the background events, taken from three regions radially symmetrical to the center of the FoV. For all parameter comparisons, we carefully checked the differences between the two wobble offsets (0.7°/1.4°), each compared with the respective off-axis MC. Due to only a negligible difference (<2%) in the acceptance between the two offsets (see Section 5.4), and to increase the statistics in the data sample, we merge them in the distributions shown. We note, however, that in the spectral (Section 6.2) and skymap (Section 6.3) analysis, each offset is handled properly as the full enclosure IRFs are used.

Figure 18 shows the distribution of the squared angular distances of the reconstructed event directions from the center of the Crab Nebula TeV gamma-ray emission (H.E.S.S. Collaboration 2020), the so-called θ² distribution. The same energy-dependent gammaness cut on 60% gamma efficiency as in the final analysis in Section 6 is applied (it was optimized on MC and turns out to yield the best signal-to-noise ratio for a source with the Crab Nebula spectrum). The comparison is shown in three bins of intensity, selected so that the rate of excess mono events in each intensity bin is comparable, demonstrating a good agreement for small shower images close to the intensity threshold as well as the extended shower images with very high intensities. Naturally, the distributions become narrower with increasing intensity, following the expected performance of the arrival direction reconstruction. For the sake of clarity in the case of mono, only results for SST-1M-1 are shown, as they are very similar for SST-1M-2. One can notice a clear improvement in the angular resolution for the stereo analysis. Small discrepancies between MC and data may be due to small mispointing (so far no pointing correction has been applied on the analysis level), or due to a small contamination of the background region with the signal events. We note that Figure 18 compares the absolute rates as no normalization between the data and MC is applied, proving overall agreement of the excess event rates between Crab Nebula data and MC. This is further demonstrated in Figure 19, showing the rate of excess events at the analysis level in different intensity bins integrated in the signal region with a radius of 0.2°, which was carefully selected to prevent the background regions of the same size from overlapping.

A good agreement is reached also in the distribution of gam-maness shown in Figure 20, proving that the Hillas parameters in the data are well reproduced in MC used for RF training. In this analysis, we consider all events in the signal and background regions (no gammaness cut applied) to avoid cropped distributions, which results in both signal and background regions being background-dominated. Distributions of gammaness in both are then subtracted, leaving only the gammaness distribution for the excess event rate. The disagreement for high gammaness close to the intensity threshold in stereo may be due to small NSB differences between MC and data, leading to different fractions of image pixels surviving the tailcuts, and it is a subject of further investigation.

Fig. 17 Left : excess map of the Crab Nebula region within a region of 1.5° × 1.5°. The best-fit position of the excess is marked with a “+” symbol, and the expected position of the Crab Nebula, as measured in H.E.S.S. Collaboration (2020), is indicated by the green “x”. The distribution was convolved with a disk kernel of 0.05°. Center: significance map across the 7.5° × 7.5° field of view centered on Crab Nebula position. The map convolved with a disk kernel with a radius of 0.11°. 0.3° region around the Crab Nebula is masked. Right: one-dimensional distribution of the significance values, including all bins (red) and OFF-region bins (gray). The OFF-region significance distribution is fitted by a Gaussian function (dashed black line).

Fig. 18 θ2 distributions for Crab Nebula excess event rates compared with point-like gamma MC in different bins of Hillas intensities. Top row: SST-1M-1 mono, Bottom row: stereo.

Fig. 19 Excess event rates binned in Hillas intensity for SST-1M-1 (left) and SST-1M-2 (right).

8 Systematic uncertainties

Event reconstruction in the IACT technique relies on MC simulations, which describe the particle interactions in the atmosphere, the production and subsequent propagation of the Cherenkov light, and finally, the response of the telescope optics and camera electronics. All these parameters are known with only limited precision and are often variable in time, and therefore contribute to the total systematic uncertainties on the reconstructed quantities. In this section, we estimate three classes of systematic uncertainties on the energy scale, flux normalization, and spectral index. Results of this study are summarized in Table 2. We note that as the statistical uncertanties of the resulting Crab Nebula SED are relatively high due to the limited size of the data sample, some of the systematic estimates are not conclusive and will be subject to a future study.

8.1 Energy scale

The absolute energy scale of the instrument is affected by any effect that contributes to the uncertainty of the image “intensity”.

For example, if the reconstructed Hillas intensities in MC are over-estimated (i.e., for a given true energy, the Hillas intensity in data is lower than in MC), some of the real events would not survive the intensity cut, leading to an underestimation of the flux, especially at low energies. The overall effect of the miscalibrated energy scale is given by the shape of the source spectrum and the effective area, and thus is quite complex to determine (see, e.g., Aleksic et al. 2016).

The Earth’s atmosphere is a critical part of the IACT technique, enabling particle showers to develop, and presenting an attenuation medium for emitted Cherenkov photons. Using MC simulated with a fixed model of the atmosphere inevitably leads to systematic uncertainty if used on data taken during variable conditions. In Section 2.4, we describe how the average VAOD used in MC was determined. To estimate the effect of its variability on the light scale, we trained an RF regressor with atmospheric transmissivity corresponding to VAOD = 0.2, and then reconstructed MC with different simulated VAOD. It turned out that a change of 0.1 roughly corresponds to a 10% change in the light scale, and correspondingly, a 10% bias in reconstructed energy. This is in line with what is naively expected since the Mie attenuation causes a I/I₀ ∝ e^-VAOD = e^−0.1 ≈ 0.9 shift in detected flux of Cherenkov photons. Our estimated systematic uncertainty of 0.05 in the measured VAOD (for the Crab Nebula dataset with exceptionally low VAOD and its variations, see Section 2.4) therefore corresponds approximately to a 5% uncertainty on the light scale.

Seasonal variations of atmospheric molecular profiles (and thus variations of the index of refraction and other relevant quantities) can have a sizable impact on the density of Cherenkov photons detected at the ground (Bernlöhr 2000). Since we use only a single representative profile for simulations covering the whole period of the Crab measurements, this represents another potential source of systematic error. To estimate this effect, we constructed three average seasonal profiles out of weekly ECMWF data, corresponding to winter, summer, and an intermediate period. We then produce three sets of simulations, one for each profile. We let the RF regressor train on the intermediate profile and then let it reconstruct the MC for all three cases.

The systematic shift in reconstructed energy and the light scale between a fixed intermediate atmosphere and a winter or summer atmosphere is roughly 1-3%, slightly depending on the energy.

The amount of detected light is also affected by the total optical efficiency of the telescope, considering mirror reflectivity, transmissivity of the camera window, and shadowing of the telescope structure. This effect is taken into account by analysis of the muon rings (Section 2.2), which allows for the tuning of the optical efficiency in MC. The systematic uncertainties induced by the discrepancies between the instrument and MC optical efficiencies are evaluated through the precision of the fit between the muon radius and the muon image charge (see Fig. 3), which is used to tune the MC optical efficiency. The selection criteria of muon events do not significantly affect the agreement between MC and data. However, for SST-1M-2, the slope of the linear fit differs between MC and data, resulting in an additional 5% uncertainty. The systematic uncertainty on the shower charge due to discrepancies between MC and data is estimated to be 4% for SST-1M-1 and 6% for SST-1M-2.

Moreover, slow degradation of the optical components of the telescopes leads to about 5% loss of optical efficiency in a year, while the MC are tuned to the average value. During the course of the Crab observation presented in this study, the optical efficiency of both telescopes decreased by about 2%, and thus we conclude the total systematics on the optical efficiency to be 5/6% for SST-1M-1/2.

The NSB-dependent effects described in Section 2.1.2 can introduce significant systematic uncertainties in the energy scale, reaching up to 20% for SST-1M-1 and 5% for SST-1M-2 if not corrected. We note that such a large difference between the two telescopes is driven by their different susceptibility to the NSB (see Section 2.1.2). However, the dependence of the electronic response on the baseline shift is accounted for and corrected in calibration. The correction factor is evaluated on the analysis level using the baseline shift determined from pedestal events taken with a frequency of 100 Hz. The precision of the correction is assessed using the muon data, where the dispersion of the average muon charge across different NSB bins after correction remains below 2%, consistent with statistical uncertainties.

Digicam readout is characterized with a negligible dead time (Heller et al. 2017). However, in case of extreme trigger rates on the level of a few GHz, which may briefly occur due to a flash from a car passing along a road near the observation site, or under extreme NSB conditions (Moon in the FoV, very Large Zenith Angle observations in Ondrejov), the camera server event buffer may reach saturation. This can lead to a reduction in the rate of stored events. The Crab Nebula sample used in this study was cleaned for these effects, and only runs with low-to-moderate NSB conditions were used in the analysis.

The reconstructed charge is given by the integral of the waveform on a predefined integration window width. We test the difference in the number of reconstructed p.e. varying the size of the integration window, even though this effect is mitigated by the integration correction (Section 4.1). We find the difference in the reconstructed pixel charge <4% for both telescopes.

Hillas intensity is one of the most important features in RF reconstruction of the energy of the primary gamma-ray photon. The number of pixels, that survive cleaning depends on the NSB level. Moreover, the amount of reconstructed charge in the individual pixels is affected by the NSB-dependent PDE, gain, and optical crosstalk, and thus any discrepancy in NSB between MC and data may lead to an offset in reconstructed energy. Although the impact of these effects is reduced by the adaptive cleaning thresholds, see Section 4.2, and voltage drop correction based on muon-ring analysis, Section 2.1.2, we tested for this by running a dedicated MC production with a very low NSB level on which we trained a set of RFs used to reconstruct testing data from our main MC production. We compared the results with those obtained if the RFs trained on MC with the same NSB level were used, resulting in 2% systematics in the energy scale.

The total systematic error on the energy scale of the Crab Nebula dataset is below 10%. We note that for periods with higher VAOD variations when no VAOD monitoring device is present on site, the energy scale systematics are dominated by its uncertainty, and thus can be higher. To evaluate the effect of the energy scale systematics on the measured SED, we perform a set of spectral analyses of the Crab Nebula data sample, each with a true energy axis in the IRFs scaled in the range of ±10% for both telescopes. We found that the spectral indices are not significantly affected by the energy scaling (ΔΓ < 0.3%), while the maximum variation of the flux normalization is Δφ₀ = ±18% for mono and Δφ₀ = ± 19% for stereo, effectively due to the shift of the SED along the energy axis. We note that using the scale invariance property of PL results in similar systematics in φ₀, while no systematic in Γ is expected.

Fig. 20 Comparison of gammaness distribution for MC simulations and Crab Nebula excess events in different intensity bins. Top row: mono SST-1M-1, Bottom row: stereo (binned in the SST-1M-1 intensity).

8.2 Flux normalization

The number of excess detected counts, and thus the source flux, depend on the estimated background counts. In the wobble observation mode, the number of background events in the source region is estimated from one or more regions of the same size at the same distance from the camera center as the source, assuming radially symmetrical acceptance. This method is therefore naturally affected by any inhomogeneity in the distribution of the background events in the FoV, leading to a systematic error in the estimated flux. Small inhomogeneities are usually caused by stars present in the FoV, locally increasing event rate, or by dead pixels, which decrease the acceptance in some parts of the camera. On top of that, small asymmetry of acceptance is expected in the case of stereo observations with only two telescopes (e.g., Aleksic et al. 2016). The latter is mitigated by wobbling, where the source and the background region positions are swapped every run. To evaluate the systematics related to background inhomogeneity, we compared the background estimated in two off-source regions for the full Crab Nebula sample, at a distance of 1.4° from the center (and equidistant to the Crab Nebula). Unfortunately, having a relatively low event rate due to a high energy threshold, the total number of background events in our sample is not large enough to make a conclusive statement about the background systematics. After applying relatively loose cuts (gamma efficiency = 90% and θ < 0.3°), we end up with 3.5 ± 2.0%, 2.3 ± 2.0%, and 2.1 ± 5.0% difference in the number of events between the two regions, for SST-1M-1, SST-1M-2, and in stereo, respectively¹¹. For stereo the difference is consistent with statistical uncertainty. We note that the effect of background uncertainty on the flux normalization depends on the signal-to-noise ratio, and while it may not be important for the Crab Nebula data, it becomes significant for faint sources.

Variable NSB in the data effectively affects the waveform variances in individual pixels, leading to differences in extracted Hillas parameters, which may result in variable acceptance for gamma rays. One may expect a general performance degradation with increasing NSB level. Even though this effect is partially mitigated with adaptive cleaning (see Section 4.2), we performed a dedicated MC study, carefully testing what is the effective area of gamma rays for both telescopes at different zenith angles and under the variable NSB conditions in the Crab Nebula sample. We found rather energy independent 10% change of the effective area. We note that this level of systematics is only valid for the presented Crab Nebula sample in this paper and may differ for different atmospheric conditions. For analysis of a different data sample spanning over a large range of NSB, one may consider using dedicated MC tuned on data in several NSB bins. The effect of simulated NSB on the intensity threshold in our data sample turns out to be at the level of a few p.e., below the intensity threshold adopted in the final analysis (45 p.e.).

To evaluate what is the effect of using RFs trained on MC with fixed level of NSB on the gamma-ray acceptance, we extended the study described in Section 8.1, and calculated the relative change in the effective area, if RFs trained on low NSB MC are used to reconstruct the baseline testing MC. We found the difference lower than 3% (energy independent) for all bins in the zenith angle, for both mono and stereo.

We combined both NSB-related factors that affect the effective area and performed a spectral analysis of the Crab Nebula sample with the effective area changed by ±10%. We confirmed that the relative difference in flux normalization is Δφ₀ = ±10% (for both mono and stereo), not changing the spectral index as expected for energy-independent effective area variation.

Due to small remaining discrepancies between MC and data, the reconstructed flux depends on the cuts applied in the analysis (tighter cuts result in underestimated fluxes). We test this effect on the final Crab Nebula SED, varying the gammaness and direction cuts in the range of ±10% the fraction of gamma-like events left in the sample after the cut. We varied the gamma efficiency (energy-dependent gammaness cut) applied on the data and corresponding IRFs between 50-70%, and θ cut between 0.17°-0.30° and 0.10°-0.18°, for mono and stereo, respectively, which correspond to the 70-90% containment, averaged over all energy bins considered in the analysis. All other parameters of the spectral analysis were left the same as described in Section 6.2. We found that the effect on flux normalization is Δφ_0,stereo = ±6%, Δφ_0,1 = ±5%, and Δφ_0,2 = ±7%.

8.3 Spectral index

Uncertainty on the background counts is the largest in low energy bins, where the number of excess gamma rays is usually smaller than the background, and thus the signal-to-noise ratio is quite low. For high energies, the situation is usually the opposite, and therefore the background systematics discussed in Section 8.2 results also in an uncertainty in the spectral index ΔΓ. Following Aleksic et al. (2016), ΔΓ can be approximately estimated as: $$\mathrm{\Delta }\mathrm{\Gamma }=2\frac{\sqrt{(5\mathrm{\%}/\mathrm{S}\mathrm{B}{\mathrm{R}}_{\mathrm{L}\mathrm{E}}{)}^2+5\mathrm{\%}/\mathrm{S}\mathrm{B}{\mathrm{R}}_{\mathrm{H}\mathrm{E}}{)}^2}}{\log (E_{\mathrm{m}\mathrm{a}\mathrm{x}}/E_{\mathrm{m}\mathrm{i}\mathrm{n}})},$$(4)

where we consider reconstruction of a spectrum in the energy range (E_min, E_max), with two energy bins only, having signal to background ratio SBR_LE/HE. 5% in the formula comes from a relatively conservative estimate of systematics in the background for both mono and stereo modes of operation (Section 8.2). For the Crab Nebula data sample, with high SBR in both energy bins, ΔΓ_mono = 0.04, ΔΓ_stereo = 0.01. We note that this is an estimate relevant only for the Crab Nebula data sample presented in this study, as 5% uncertainty on the background counts is only a limit given by relatively small statistics.

Using the results of the analysis performed in Section 8.2, we can evaluate how the spectral index is affected by the analysis cuts. Varying the gamma efficiency and the size of the signal region, we found that for stereo and SST-1M-1 in mono, the spectral index change within 3% (ΔΓ₁ = +0.08/-0.05, ΔΓ₂ = +0.03/-0.08); whereas for SST-1M-2 in mono, it is within 6% (ΔΓ₂ = +0.0/-0.2).

9 Summary and conclusions

We evaluated the low-altitude performance of the SST-1M telescopes in mono and stereo modes using MC simulations and Crab Nebula data taken during commissioning from September 2023 to March 2024. The MC model of the telescopes was carefully tuned to accurately represent both telescopes, and the atmospheric conditions at the Ondřejov observatory. Due to the dependence of the gain, PDE, and the optical crosstalk on the NSB level, special care was taken with the telescope calibration and tuning of the simulated NSB level to reflect low to moderate NSB conditions in Ondřejov.

We obtained the trigger threshold (for a source with the spectrum of the Crab Nebula) at low zenith angles to be E_T,mono = 0.6 TeV, and E_T,stereo = 0.7 TeV, rising up to E_T,mono = 1 TeV, and E_T,stereo = 1.3 TeV at the analysis level. At the low zenith angle, the energy resolution in mono is about 20% at low energies near the energy threshold and reaches about 15% at higher energies. In stereo, the energy resolution is between 10-15%. The angular resolution is better than 0.18° above the energy threshold in mono and reaches 0.10° in stereo. The flux sensitivity in stereo shows an improvement of about a factor of 2 compared to the mono reconstruction. The integral sensitivity above the energy threshold in stereo reaches about 7% C.U. in 50 hours. We also show the zenith angle dependence of these quantities, showing an improvement at large energies with increasing zenith angle at the expense of increasing energy threshold.

Additionally, we performed an MC study of the off-axis performance, which shows a remarkably flat acceptance up to about 2.5° offset, where it drops by 10%. The angular resolution and the integral sensitivity show similar behavior, making the SST-1M an ideal instrument for multi-TeV observations of extended gamma-ray sources or for the follow-up of poorly localized transients.

We used the low-zenith-angle Crab Nebula data to demonstrate a good MC-data agreement, showing that the real observation matches the expected performance. The derived Crab Nebula SED is in good agreement with the results of other observatories within the reported uncertainties. We carefully studied the main sources of systematic uncertainties and evaluated their impact on the measured SED using dedicated MC simulations. The total systematic uncertainties were found to be <10% in the energy scale, 23% in the flux normalization, and about 3-6% in the spectral index. We note that the most important effect on the flux normalization systematics is the uncertainty in the energy scale, dominated by the uncertainty in VAOD and optical efficiency. While the former can be mitigated with the adoption of on-site atmospheric monitoring of VAOD or the so-called Cherenkov transparency correction (Hahn et al. 2014; Stefanik et al. 2019), improvements to the muon analysis have to be investigated to make improvements to the latter. The second largest contribution comes from the variable NSB in the data, leading to an effective area mismatch with the MC simulated with a fixed NSB level. Improvements in the analysis methods to reduce these effects will be the subject of future studies.

Acknowledgements

This publication was created as part of the projects funded in Poland by the Minister of Science based on agreements number 2024/WK/03 and DIR/WK/2017/12. The construction, calibration, software control and support for operation of the SST-1M cameras is supported by SNF (grants CRSII2_141877, 20FL21_154221, CRSII2_160830, _166913), by the Boninchi Foundation, and by the Université de Genève, Faculté de Sciences, Département de Physique Nucléaire et Corpusculaire. The Czech partner institutions acknowledge support of the infrastructure and research projects by Ministry of Education, Youth and Sports of the Czech Republic (MEYS) and the European Union funds (EU), MEYS LM2023047, EU/MEYS CZ.02.01.01/00/22_008/0004632, CZ.02.01.01/00/22_010/0008598, Horizon Europe MSCA COFUND Physics for Future 101081515, and Czech Science Foundation, GACR 23-05827S. Author contribution. J. Juryšek: paper project coordination, MC performance, and validation, Crab Nebula SED, study of the systematics uncertainties. T. Tavernier: Calibration, muon analysis and voltage drop correction, optical efficiency, Crab Nebula skymap. V. Novotny: MC production, low-level MC-data tuning. J. Blazek: VAOD and atmospheric effect. D. Mandât, M. Pech: telescope operation and data acquisition. All these authors contributed to the paper drafting. The rest of the authors contributed to one or more of the following: design, construction, maintenance, SW development, and discussion and approval of the draft.

References

All Tables

All Figures

Fig. 1 Comparison between data and MC of the normalized (to area) ADC charge distributions extracted from the dark runs for the two telescopes. To avoid overlaps of the two curves, the SST-1M-2 distribution is multiplied by 0.1. Each count is calculated from the 50 samples (200 ns) long waveform in the same manner as done for the triggered data. See the explanation in the text. An asymmetric distribution of the electronic noise (more significant in SST-1M-2) can explain the slight disagreement between the first and the second peak.

In the text

Fig. 2 Left: relation between muon ring radius and ADC counts under various NSB conditions for SST-1M-1. Right: behavior of measured muon ring charge as a function of the baseline shift induced by the NSB. Blue dotted and orange dashed lines represent the expected behavior given in Nagai et al. (2019) for SST-1M-1 and SST-1M-2, respectively. Filled blue (SST-1M-1) and empty orange (SST-1M-2) markers represent the fitted estimation of the relative ADC counts for 1.15° radius muon images.

In the text

Fig. 3 Charge of muon images as a function of the fitted ring radius for SST-1M-1 and SST-1M-2. The solid markers represent observational data collected from September 2023 to April 2024, during the first stereo Crab observation campaign (blue for SST-1M-1 and orange for SST-1M-2). Measured data agree with the MC simulation, shown with open markers and dashed lines. Notable differences between SST-1M-1 and SST-1M-2 arise from differences in their optical designs.

In the text

Fig. 4 Evolution with time of the charge in the muon ring pixels for fixed radius of r = 1.2° for SST-1M-1 (blue points) and SST-1M-2 (orange points). Each data point corresponds to one single observation night. The green-shaded area indicates the time period of the Crab Nebula observation campaign presented in this work. The linear fits (dashed for SST-1M-1, dotted for SST-1M-2) enable the estimation of the relative decrease in optical efficiency over time.

In the text

Fig. 5 Left : distribution of the averaged standard deviation of the baselines for SST-1M-1 (dashed line) and SST-1M-2 (solid line). The vertical lines represent the value assumed in the MC model of SST-1M-1 (94 MHz/pixel) and SST-1M-2 (120 MHz/pixel). The vertical axis represents the number of hours with a baseline standard deviation (on the horizontal axis) in the final Crab Nebula data sample. Right: comparison between the MC-predicted and pedestal events-measured p.e. distributions under the slightly higher than typical NSB conditions during the Crab Nebula observations.

In the text

Fig. 6 Gini importance for different RFs for the SST-1M-1. The features used for the reconstruction: log of the intensity, width, length, width/length ratio (wl), timing slope, skewness, kurtosis, leakage, and coordinates of the shower center of gravity in the FoV (x,y). Height of the shower maximum (h max) and distance of the impact point from the telescope (impact dist) are used only in the stereo reconstruction. Upper left: gamma-hadron classifier. Upper right: energy regressor. Bottom left: disp norm regressor. Bottom right: disp sign classifier (Only used in mono reconstruction).

In the text

	Fig. 7 ROC curve of the RF gamma-hadron classifier for 30° (solid line) and 50° (dotted line) zenith angles. Line colors represent different regimes of observations.
In the text

Fig. 8 Left: differential rate of point-like gamma rays with Crab Nebula spectrum in mono and stereo for 20° zenith angle. Event rates for different stages of the analysis are shown: all triggered events (solid line), all events that survived cleaning (dashed line), events that survived quality cuts (dash-dotted line), and events that survived θ2 and gammaness cuts used in the analysis (dotted line). Right: effective area for 20° zenith angle.

In the text

	Fig. 9 Energy resolution (left) and energy bias (right) for 30° (solid line) and 50° (dotted line) zenith angles evaluated using on-axis point-like gammas for testing. Line colors represent different regimes of observations. For the sake of clarity, only results for SST-1M-1 are shown in mono as both telescopes are very similar.
In the text

	Fig. 10 Angular resolution for 30° (solid line) and 50° (dotted line) zenith angles evaluated on testing on-axis point-like gammas. Line colors represent different regimes of observations. For the sake of clarity, only results for SST-1M-1 are shown in mono, as both telescopes are very similar.
In the text

Fig. 11 Differential sensitivity at 5σ confidence level for 30° (solid line) and 50° (dotted line) zenith angles evaluated on testing on-axis point-like gammas and diffuse protons. Line colors represent different regimes of observations. For the sake of clarity, only results for SST-1M-1 are shown in mono, as both telescopes are very similar. We note that in contrast to stereo, the additional condition requiring at least ten excess events does not affect the last energy bins for the 50° zenith angle in mono, as the number of excess events is satisfied by the 5σ condition only (see Section 5.3 for more details).

In the text

	Fig. 12 Simulated gamma-ray rate from the Crab Nebula and simulated protons weighted by the CR spectrum in the signal region after analysis cuts as a function of the offset at 20° zenith angle. Only events with energies above the ET are taken into account.
In the text

	Fig. 13 Integral sensitivity for 50 hours of observation (left), and angular resolution for a source with the Crab Nebula spectrum (right) as a function of the offset at 20° zenith angle. All events with energies above the ET are integrated.
In the text

	Fig. 14 Left: run-averaged rate of mono events with intensity above 200 p.e., corrected for cos(zenith angle). Center: run-averaged fraction of pedestal events that survived the tailcut cleaning. Right: run-averaged fraction of pixels with increased tailcuts. The solid blue line marks the cuts used in the run selection.
In the text

Fig. 15 Per-run rates of the event intensities (i.e., only events which survived cleaning contribute to the distributions) for SST-1M-1 mono (left) and SST-1M-2 mono (center), and SST-1M-1 stereo (right). The distributions shown are corrected for cos(zenith). Dotted lines mark all runs taken, solid lines mark runs after selection cuts. The dashed line shows the distribution of MC protons re-weighted on the proton + helium spectrum. The solid blue line marks the intensity threshold of 45 p.e. introduced in the higher-level analysis.

In the text

Fig. 16 Left: SED of the Crab Nebula measured with SST-1M telescopes. Different color bands represent the best-fitting spectra for the SST-1M-1 mono (blue), SST-1M-2 mono (yellow), and SST-1M stereo (green) datasets. Right: stereo flux points derived from the PL spectral model. SEDs from other experiments are shown for comparison in the energy range of their validity (left), and corresponding flux points (right) - MAGIC (Aleksić et al. 2015; MAGIC Collaboration 2020), H.E.S.S. (Aharonian et al. 2024), VERITAS (Meagher & VERITAS Collaboration 2015), LST-1 (Abe et al. 2023b), LHAASO (LHAASO Collaboration 2021), and HAWC (Abeysekara et al. 2019).

In the text

Fig. 17 Left : excess map of the Crab Nebula region within a region of 1.5° × 1.5°. The best-fit position of the excess is marked with a “+” symbol, and the expected position of the Crab Nebula, as measured in H.E.S.S. Collaboration (2020), is indicated by the green “x”. The distribution was convolved with a disk kernel of 0.05°. Center: significance map across the 7.5° × 7.5° field of view centered on Crab Nebula position. The map convolved with a disk kernel with a radius of 0.11°. 0.3° region around the Crab Nebula is masked. Right: one-dimensional distribution of the significance values, including all bins (red) and OFF-region bins (gray). The OFF-region significance distribution is fitted by a Gaussian function (dashed black line).

In the text

	Fig. 18 θ2 distributions for Crab Nebula excess event rates compared with point-like gamma MC in different bins of Hillas intensities. Top row: SST-1M-1 mono, Bottom row: stereo.
In the text

	Fig. 19 Excess event rates binned in Hillas intensity for SST-1M-1 (left) and SST-1M-2 (right).
In the text

	Fig. 20 Comparison of gammaness distribution for MC simulations and Crab Nebula excess events in different intensity bins. Top row: mono SST-1M-1, Bottom row: stereo (binned in the SST-1M-1 intensity).
In the text