V-shaped Cherenkov images of magnetically-separated gamma-rays

TL;DR

This paper investigates V-shaped Cherenkov images caused by magnetically-separated gamma-ray events, proposing a new geometrical energy estimation method and analyzing its potential to improve CTA measurements.

Contribution

It introduces a novel geometrical method for energy estimation of V-shaped Cherenkov images and evaluates its effectiveness alongside classical methods.

Findings

V-shaped images result from geomagnetic deflection of early particle pairs

The new geometrical method provides an independent energy estimate

Combining methods can help validate systematic uncertainties

Abstract

Cherenkov Telescope Array (CTA) is an upcoming instrument that will start a new generation of atmospheric Cherenkov telescopes. CTA is expected not only to provide an unprecedented sensitivity in the tens of GeV to hundreds of TeV range, but also to considerably improve the systematic uncertainties of the measurements. We study the images registered by Cherenkov telescopes from low energy gamma rays with its first interaction in the upper parts of the atmosphere. The images show a characteristic separation due to the deflection of the first e- e+ pair in the Geomagnetic Field. We evaluate the performance of the standard stereoscopic analysis for such events. We derive also a novel method for energy estimation of V-shaped events based purely on geometrical properties of the image. We investigate the potential of combining the classical energy estimation and the novel method for independent validation of the systematic shifts in the energy scale of Cherenkov telescopes and discuss the limitations of such analysis.

Equations

$$d_1=\frac{6.4[\mathrm{km}]}{\cos \mathrm{ZD}}\ln \left(1+\frac{X_0\cos \mathrm{ZD}}{1340[\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}]\exp (-0.155{\mathrm{H}}_1/\mathrm{km})}\right),$$

$${\Psi }_{+-}\approx 2d_1/R_L=0.1^{\circ }(d_1/4\mathrm{km})(E_0/100\mathrm{GeV}{)}^{-1}(B_{\perp }/0.4\mathrm{G}),$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

For the last 30 years the Cherenkov telescopes have been the main driver of studies of galactic and extragalactic TeV gamma-ray sources. The IACT (Imagining Air Cherenkov Telescopes) technique is based on the measurement of Cherenkov photons produced in the atmosphere by charged relativistic particles forming an Extensive Air Shower (EAS) initiated by the primary gamma-ray. The two-dimensional angular distribution of Cherenkov light appears on the telescope camera as the shower image. The parametrization (Hillas, 1985) of the image allows us to estimate the energy and direction of the primary particle, and to exclude a large fraction of the dominating hadronic background. The rapid progress of both the hardware developments and new analysis methods resulted in great improvement in the sensitivity of Cherenkov telescopes (see e.g. Fig. 17 of Aleksić et al., 2016b). Currently, three major IACT instruments are in operation: H.E.S.S. (Aharonian et al., 2004), MAGIC (Aleksić et al., 2016a) and VERITAS (Weekes et al., 2002). Already with the current generation of IACTs it is possible to obtain for bright flares precision measurements of nightly fluxes with statistical uncertainty of only a few per cent (see e.g. Aharonian et al., 2007; Furniss et al., 2015; Ahnen et al., 2016). Nevertheless, the Cherenkov telescopes are limited by relatively large systematic uncertainties, usually the dominating one is the uncertainty in the energy scale. It is caused mainly by the uncertainty in the atmospheric production and absorption of the Cherenkov photons as well as lack of a “test beam” for IACT instruments. Despite various calibration methods (inter-telescope calibration: Hofmann, 2003, muon analysis: Vacanti et al., 1994, flux comparisons at different thresholds: Aleksić et al., 2016b, LIDAR corrections: Fruck & Gaug, 2015; Devin et al., 2019) it is still at the level of $\sim 15\%$ (see e.g. Aharonian et al., 2006; Abeysekara et al., 2015; Aleksić et al., 2016b), much larger than of satellite experiments (c.f. 2% energy scale systematic uncertainty estimated for Fermi-LAT, Ackermann et al., 2012b). Such an uncertainty for a typical gamma-ray source can easily produce a $\sim 30\%$ systematic error in the absolute normalization of the flux, with a more severe effect for steep spectrum sources and the energies close to the threshold.

The upcoming Cherenkov Telescope Array (CTA) (Actis et al., 2011; Acharya et al., 2013) will start a next generation of large IACT arrays. CTA is designed to study gamma-ray sources from a few tens of GeV to hundreds of TeV with unprecedented sensitivity. To reach this goal it will be composed of sub-arrays of telescopes of different sizes: SST, MST and LST (small, medium and large sized telescope, respectively). To fully exploit the scientific potential of CTA, the increase in the statistical accuracy has to be matched also by efforts to lower systematic uncertainty of the measurements. In fact, various studies are being performed to lower the systematic uncertainty of the CTA: using atmospheric monitoring (Gaug et al., 2015), muon calibration (Brown et al., 2015), cosmic ray electron spectrum (Parsons et al., 2016), inter-telescope calibration (Mitchell et al., 2016).

An interesting method to verify the energy calibration of ground instruments exploits the magnetic field of the Earth (i.e. Geomagnetic Field, GF) to measure the deflection of charged cosmic rays and the resulting shift of the shadow produced by the Moon or the Sun. Even while such a method proved successful for surface arrays of detectors like MILAGRO or HAWC (Christopher, 2011; Abeysekara et al., 2013), its application to IACTs is difficult due to strong background of moonlight. However, the Moon shadow measurements were performed with CLUE experiment, an array of Cherenkov telescopes sensitive only in UV range where the atmospheric absorption is high (Alexandreas et al., 1998). Moreover, so far GF has been considered as a nuisance to the IACT technique, due to its effect on the images of the showers (see e.g. Bowden et al., 1992; Commichau et al., 2008; Szanecki et al., 2013).

In this work we investigate the feasibility of a novel method for verifying the energy calibration of CTA. The method is based on a search of characteristic images created by separation of the first $e^{-}e^{+}$ pair in the upper parts of the atmosphere by GF. In Section 2 we perform basic analytical estimations to evaluate properties of such events. In Section 3 we describe the used simulation and analysis pipeline. Results of the simulations are presented in Section 4. The paper is concluded with discussion in Section 5.

2 Magnetically separated events

In order to evaluate what properties such events would have let us consider a $E_{0}\sim 100$ GeV gamma ray entering the atmosphere at a low zenith angle. On average after $\sim 9/7$ radiation lengths, i.e. at the height of $H_{0}\sim 22$ km it will interact with atmospheric nuclei producing $e^{-}$ and $e^{+}$ with energies $E_{-}=f_{e^{-}}\,E_{0}$ and $E_{+}=(1-f_{e^{-}})\,E_{0}$ respectively. Initially $e^{-}$ and $e^{+}$ will be travelling nearly in the same direction as the primary gamma ray (the angle between the direction of $e^{\pm}$ and primary gamma rays will be of the order of arcsecond). However as they move through the atmosphere they will be more and more deflected by the GF, loosing at the same time energy by production of multiple Bremsstrahlung gamma rays. To estimate a characteristic deflection we can assume that the deflection occurs on a distance, $d_{1}$, corresponding to one radiation length, over which $e^{+}e^{-}$ will lose most of its energy. For the further calculations we will assume a simple exponential atmosphere with thickness $X=1340\cdot e^{-0.155H/\mathrm{km}}/\cos\mathrm{ZD}\mathrm{\,[g\,cm^{-2}]}$111For close to zenith observations, this formula describes the La Palma atmosphere model used in standard CTA simulations with accuracy better than 5% between height of 6 and 35 km., where $\mathrm{ZD}$ is the zenith angle of the shower. We obtain:

[TABLE]

where $X_{0}=37\mathrm{[g\,cm^{-}2]}$ is the radiation length in air, and $H_{1}$ is the height of the first interaction. For $f_{e^{-}}\approx 0.5$ we obtain the total deflection angle between $e^{+}$ and $e^{-}$ of:

[TABLE]

where $R_{L}$ is the Larmor radius of $e^{+}$/$e^{-}$, and $B_{\perp}$ is the component of GF perpendicular to the shower axis. Such a separation is within reach of the angular resolution of the CTA telescopes. In fact, the separation can be even more pronounced if the first interaction occurs higher in the atmosphere (resulting in longer $d_{1}$), and also taking into account that the primary $e^{+}$ and $e^{-}$ will loose energy during the deflection. Nevertheless, as most of the $e^{+}e^{-}$ energy is lost early, the Eq. 2 is still a reasonable approximation of the effective deflection. The two subshowers started by $e^{+}$ and $e^{-}$ can be then in principle angularly separated. The angular separation of the two subshowers also shifts the positions of the Cherenkov photons on the ground, which can be roughly estimated as $\Psi_{+-}H_{1}/\cos\mathrm{ZD}\approx 40$ m for the above scaling values. Such a separation is a factor of a few below the radius of the shower lightpool ($\sim 120$ m). Hence we conclude that usually the same telescope would be able to see both subshowers, however propagating in a slightly different direction. Therefore such events observed under favourable conditions would produce characteristic V-shaped images. The difference of such V-shaped images to the magnetically broadened events studied e.g. by Bowden et al., 1992; Commichau et al., 2008; Szanecki et al., 2013 should be stressed. These authors focus on the global effect of broadening of the image, mostly caused by the deflection of copious, lower energy $e^{+}e^{-}$ pairs produced in far generations in GF over shorter distances. Due to large number of individual $e^{+}e^{-}$ pairs with different energies and different distances from the shower axis, such broadening would not cause a clear separation of the image in two components, which naturally occurs in the scenario presented here. Note however, that each of the two sub-images in V-shaped events would be broadened by this effect, which is automatically taken into account in our simulations.

In Fig. 1 we show an example of an event with $e^{+}$ and $e^{-}$ subshowers separated as seen by 4 LST telescopes.

In each pixel we plot the number of photoelectrons (phe) reconstructed from the signal recorded in that pixel. Night sky background and electronic noise is simulated independently in $e^{+}$ and $e^{-}$ and total images. The light produced in $e^{+}$ and $e^{-}$ subcascade is clearly shifted in the direction perpendicular to the magnetic field vector (i.e. for the simulated geometry in the horizontal direction in the camera). The total image (i.e. both the $e^{+}$ and $e^{-}$ subcascades together) of the shower in T4 (and to a lesser extent also in T2 and T3) is V-shaped. An even clearer separation can be achieved for lower energy showers and for showers that have started earlier in the atmosphere (see Fig. 2). Note however that showers starting early in the atmosphere will produce dimmer images that are more difficult to reconstruct.

3 Simulations and data analysis

We perform simulations of gamma rays using CORSIKA version 7.5 (Heck et al., 1998). We have simulated an array of 4 LST telescopes located in the La Palma site. We simulate gamma rays with zenith angle of $20^{\circ}$ and azimuth corresponding to shower coming from geomagnetic North. Therefore the angle between the direction of the shower and GF is $72^{\circ}$, nearly maximizing the GF effects. The response of the telescopes was simulated using sim_telarray (Bernlöhr, 2008) using settings of the so-called CTA Prod 3 MC (Acharyya et al., 2019). Both programs were modified to track additional information of gamma-ray induced showers. Namely the particle type (most commonly $e^{+}$ and $e^{-}$ for a primary gamma ray shower) of the particles produced in the first interaction is propagated to the further generations of the cascade starting from those particles, and to the Cherenkov photons produced by the latter. Therefore, for each emitted Cherenkov photon we know if it originated in the subcascade caused by $e^{+}$ or $e^{-}$ from the first interaction. The generated MC simulations are described in more detail in Sitarek et al. (2018b). The extraction of signal amplitudes from simulated waveforms, image cleaning (i.e. keeping only pixels with significant signals above NSB fluctuations) and its parametrization, gamma/hadron separation and classical stereo reconstruction and energy estimation is done using MARS/Chimp chain (Zanin et al., 2013; Aleksić et al., 2016b; Sitarek et al., 2018a). Since the deflection of the two subshowers by the GF is roughly inversely proportional to the energy of the primary gamma we simulate only the LST subarray, where the strongest effect is expected.

4 Results

In this section we present results of two types of simulations. In the first, dedicated simulations, we exploit the extra knowledge of separation of shower into $e^{+}$ and $e^{-}$ subshowers to study how the separation into two subshowers affects the image. Next we apply the V-shaped events selection cuts to the second, standard set of simulations, in which the composition of the image from $e^{+}$ and $e^{-}$ subshowers is unknown.

4.1 Comparison of $e^{+}$ and $e^{-}$ images

In order to investigate the separation effect of the $e^{+}$ and $e^{-}$ subcascades we analyze images composed by Cherenkov photons from only one of the two subcascades produced in the first interaction. We clean separately images produced by $e^{+}$ and $e^{-}$ subcascades, calculate their Hillas parameters and exclude those in which at least one of the subcascades did not produce at least 50 phe signal. We define $\Delta\mathrm{COG_{X}}$ as the separation distance of the $e^{+}$ and $e^{-}$ subshower images, measured perpendicularly to the direction of the magnetic field projected onto the camera. For the simulated here case of showers coming from geomagnetic North the positive value of $\Delta\mathrm{COG_{X}}$ means that the images are separated in the direction consistent with deflection in GF. In Fig. 3 we investigate how $\Delta\mathrm{COG_{X}}$ depends on the energy of the primary gamma ray and its height of the first interaction.

The image to image fluctuations are large, however a global shift of the distribution is clearly seen. At TeV energies the separation is small, $\sim 0.05^{\circ}$. It increases going down in energy, to $\sim 0.1^{\circ}$ at about 150 GeV, a consistent value with the result of the order of magnitude calculations of Eq. 2. Interestingly, at the lowest energies the average shift is getting smaller again (see the left panel of Fig. 3). This can be understood in terms of the dependence of the image separation on the height of the first interaction (the right panel of Fig. 3). Only those of the lower energy showers that are generated deeper in the atmosphere will allow enough light to be gathered by the telescopes for the triggering of the event and its consecutive reconstruction. Such events produced deeper in the atmosphere will be affected by a smaller deflection in the GF. Limiting the calculations to only events starting high in the atmosphere ($>25$ km a.s.l., which at energy of $\sim 100$ GeV corresponds to roughly half of the events) shows a monotonic behaviour of the separation angle down to energies of  30 GeV (see the middle panel of Fig. 3).

4.2 Selection of the v-images and rates

Since the images of low energy showers are composed of a small number of pixels, searching for two parts of an image by e.g. a fit with two two-dimensional Gaussian distribution, which has a large number of free parameters can lead to unstable, non-converging fits. Therefore, in order to preselect interesting images for more detailed analysis we apply, telescope-wise, a fast analytic “double Hillas” procedure.

First we calculate the Hillas parameters of the full image. Next, the image is divided by a line crossing the source position on the camera and the COG of the image. The pixels forming the image are divided into two groups depending on which side of the line their center is, and the Hillas parameters are calculated independently for each of the groups. If the image is a classical ellipse such a procedure will approximately cut it along its main axis resulting in two parts of the image separated by $\Delta\mathrm{COG_{X}}\lesssim 0.1^{\circ}$ (of the order of the Width parameter of the image). The average angle between their main axes can be estimated as a ratio of $\Delta\mathrm{COG_{X}}$ to Dist (the angular distance of the COG of the image from the direction of the primary gamma ray), i.e. $\Delta\delta\sim\arcsin(0.1^{\circ}/1^{\circ})\approx 6^{\circ}$. On the other hand, in the case of images like T4 in Fig. 1 or Fig. 2 double-Hillas analysis reveals the two components and results in a larger separation of COGs and main axes of the two parts of the image. To reject fake islands in the image due to night sky background remaining after the cleaning and also to limit the contribution of the hadronic background, we keep only the events in which both parts of the image separately are pointing towards the true source position. We achieve it by cutting in the so-called Miss parameter, i.e. the distance between the main axis of the subshower image and the direction of the source.

In order to check the viability of the procedure described above and to tune the values of the cuts in $\Delta\mathrm{COG_{X}}$ and $\Delta\delta$ in Fig. 5 we compare the $\Delta\delta$ vs $\Delta\mathrm{COG_{X}}$ plot obtained from special MCs where $e^{+}$ and $e^{-}$ subcascades were processed separately (top panel) with the one obtained from regular MCs in which double-Hillas analysis was used.

Note that in the latter case the subimages from $e^{+}$ and $e^{-}$ subcascades cannot be distinguished, i.e. only absolute value of $\Delta\mathrm{COG_{X}}$ can be derived. As shown already in Fig. 3, due to the deflection in GF, $\Delta\mathrm{COG_{X}}$ is shifted towards positive values. The most clear V-shaped events are located in the top right part of the top panel of Fig. 5 (i.e. large $\Delta\delta$ and large positive $\Delta\mathrm{COG_{X}}$). The events with large $\Delta\delta$ but large negative $\Delta\mathrm{COG_{X}}$ are events in which large fluctuations of the shower mimick the separation due to GF deflection. As the number of events with $\Delta\delta>20^{\circ}$ and $\Delta\mathrm{COG_{X}}<-0.15^{\circ}$ is much smaller than the number of events for which $\Delta\delta>20^{\circ}$ and $\Delta\mathrm{COG_{X}}>0.15^{\circ}$, we conclude that the background of such mimicking events will be small. The bottom panel of Fig. 5 obtained with double-Hillas analysis of full images show a similar shape as the folded along the $\Delta\mathrm{COG_{X}}=0$ axis upper panel. As expected the classical elliptical events are located at small values of $\Delta\mathrm{COG_{X}}$ and $\Delta\delta$, while the more interesting, V-shaped events will produce larger values of $\Delta\mathrm{COG_{X}}$ and $\Delta\delta$. We define image as V-shaped if it fulfills the following conditions: $Miss_{Img1}<0.15^{\circ}$, $Miss_{Img2}<0.15^{\circ}$, $\Delta\mathrm{COG_{X}}>0.15^{\circ}$ and $\Delta\delta>20^{\circ}$. $Miss_{Img1}$ and $Miss_{Img2}$ are the $Miss$ parameters calculated independently for the two separated by the double-Hillas analysis parts of the image.

In Fig. 6 we compare the distribution of the true impact parameter (calculated with respect to each telescope) for V-shaped and not-V-shaped images.

For V-shaped images the distribution is shifted to lower values and peaks around 60 m. It is caused by the geometrical effect that the same separation angle $\Psi_{+-}$ between the $e^{+}$ and $e^{-}$ subshowers for different impact parameters will result in different $\Delta\delta$ affecting also $\Delta\mathrm{COG_{X}}$ angle. For example if the impact parameter is zero, the two subshowers will propagate in opposite direction on the camera ($\Delta\delta=180^{\circ}$), maximizing also $\Delta\mathrm{COG_{X}}$.

In Fig. 7 we show the rates of such preselected events in which at least one/two images are classified as V-shaped.

The rates are normalized to a source with a Crab-like spectrum (including its curvature following Aleksić et al., 2015). The lowest energies produce too dim images for both components to be clearly visible. On the contrary, the highest energies produce large images in which the $e^{+}$ and $e^{-}$ subcascades are not-separatable. Therefore we limit the analysis to events with estimated energy between 30 and 300 GeV. At 100 GeV about 1/10 of the events show at least one image with separation of the two components in the double Hillas analysis larger then the applied cuts. At the same energy about 1/100 of the events show a clear indication of being V-shaped in at least two individual telescope images. For observations of a source with a Crab like spectrum the total rate of the V-shaped events surviving the selection criteria in at least two images is 50 per hour. For the further analysis we consider events with V-shape images present in at least two LST telescopes.

4.3 Toy model

According to Eq. 2 the deflection of the two images depends on the energy of the primary gamma ray. The classical energy estimation is based mainly on the amount of Cherenkov light for the reconstructed stereoscopic parameters (the impact parameter and the height of the shower maximum). Hence, comparing the classical energy estimation of those events with the developed here independent geometrical energy estimation has a potential for validation of the total light throughput of the telescopes and this in turn validates the energy calibration. Unfortunately $\Psi_{+-}$ depends also on the height of the first interaction and on the separation of the primary energy into $e^{+}e^{-}$. For individual events those are poorly known, however if a large sample of V-shaped events is gathered, this problem can be dealt with on a statistical basis.

In order to apply such a calibration procedure to Cherenkov telescopes, precise enough knowledge of the GF is needed. In fact, the standard simulation program of EAS, CORSIKA, has a somewhat simplified treatment of the GF, neglecting the height dependence. Nevertheless, according to the WMM2015 model222https://www.ngdc.noaa.gov/geomag/WMM/ the change of the GF components in our range of interest, i.e. between 10 and 40 km a.s.l. is only about 1.5% for the La Palma site.

In order to test the potential of V-shaped events to calibrate/validate the energy scale of LST we developed a simple geometrical model of such events. The model assumes a knowledge on the direction of the shower $\vec{v}_{0}$ (azimuth and zenith angles: $\phi,\theta$), and local GF, $\vec{B}$, and location of all the telescopes (including relative height differences). In addition to this, it has 5 free parameters: the core position on the ground $(x_{0},y_{0})$, the height of the first interaction $H_{0}$, energy of the gamma ray $E_{0}$, and the fraction $f_{e^{-}}$ of the energy taken by the first generation $e^{-}$. All of the parameters (except of $f_{e^{-}}$) can be left free to vary in the fitting procedure or it can be estimated from other methods. In particular the core position can be estimated using stereoscopic reconstruction of the shower (Hofmann et al., 1999). The height of the first interaction, even while poorly constrained, can be estimated with the closest pixel method (Sitarek et al., 2018b) or with a model analysis (de Naurois & Rolland, 2009; Parsons $\&$ Hinton, 2014).

Using $\vec{v}_{0}$, $(x_{0},y_{0})$ and $H_{0}$ we calculate first the point in space, $P_{0}$, where the first interaction happens. Starting from $P_{0}$ the first generation $e^{+}$ and $e^{-}$ pair diverges. We assume that each of the leptons traverses one radiation length $d_{1}$, which geometrical length is determined from $H_{0}$ and $\theta$ using a simple exponential profile of the atmosphere. The deflection angle of $e^{+}$ and $e^{-}$ is computed as $\Psi_{+}=d_{1}/R_{L,+}$ and $\Psi_{-}=d_{1}/R_{L,-}$ respectively, where $R_{L,+}$ and $R_{L,-}$ are their Larmor radii computed using energies $(1-f_{e^{-}})E_{0}$ and $f_{e^{-}}E_{0}$ and perpendicular part of $\vec{B}$. We calculate the directions of the deflection of the $e^{+}$ and $e^{-}$ subshowers as the unit vectors $\vec{n}_{\pm}=\pm\frac{\vec{v}_{0}\times\vec{B}}{|\vec{v}_{0}\times\vec{B}|}$. The direction of the $e^{\pm}$ subshowers is then computed as $\vec{s}_{\pm}=\frac{\vec{v}_{0}+\Psi_{\pm}\vec{n}_{\pm}}{|\vec{v}_{0}+\Psi_{\pm}\vec{n}_{\pm}|}$. The starting points of the $e^{\pm}$ subshowers are then assumed as $P_{\pm}=P_{0}+0.5(\vec{s}_{\pm}+\vec{v}_{0})d_{1}$. Each subshower propagating from point $P_{\pm}$ along $\vec{s}_{\pm}$ direction will then produce a line in the coordinates system on the Cherenkov telescope camera. Therefore, for a given set of starting parameters we obtain 2 lines per each telescope.

It should be noted, that the toy model does not take into account several effects. In particular the stochastic behaviour of Bremsstrahlung, energy losses of the first $e^{\pm}$ pair affecting the deflection, lateral distribution of the electrons in the subshowers and selection bias of the events triggered by the telescopes can modify the resulting image.

4.4 Standard stereoscopic reconstruction for V-shaped events

In order to evaluate how the standard stereo reconstruction deals with the V-shaped events we compare the basic performance parameters for the events with (a) two or more, (b) one and (c) no V-shaped images. Note that the simulations were performed for showers coming from North, where the effect of GF worsens the performance of the Cherenkov telescopes (see e.g. Hassan et al., 2017). In Fig. 8 we show the comparison of angular and energy resolution for different classes of events.

The energy resolution we calculate as a second central moment of the distribution of $(E_{est}-E_{true})/E_{true}$, i.e. including non-Gaussian tails. It should be noted that for those comparisons we do not apply any additional quality cuts (e.g. in number or telescopes used in the reconstruction), nor cuts in the gamma/hadron separation parameter. Such cuts would improve the absolute values of the performance parameters. For the calculation of the energy resolution we apply however a cut in the reconstructed position of the source. The performance parameters for events with some of the images V-shaped are comparable to events with classical elliptical images in all the telescopes. This can be understood taking into account that the most clear V-shaped images occur if the impact parameter of the shower to the telescope is small (compare also Fig. 6) and both subcascades develop on opposite sides of the telescope. In such a case the event is normally relatively well contained within the array of 4 LST telescopes and there are enough non-V images for efficient reconstruction. Moreover, V-shaped images will contribute with lower weight in the energy reconstruction due to their unusual values of Hillas parameters (c.f. the Chimp energy estimation procedure explained in Sitarek et al., 2018a).

In Fig. 9 we show the distribution of the gamma/hadron separation parameter Hadronness for different classes of events.

V-shaped events when parameterized as Hillas ellipses have larger Width, and in case of $\Delta\delta>90^{\circ}$ also larger Length making them more hadron-like. Therefore the distributions for events with V-shaped images are shifted towards larger Hadronness values. Nevertheless the effect is small even for the case of two V-shaped images in one event. The value of Hadronness cuts that keeps 90%, 80% and 60% of the events with no V-shaped images keeps also respectively 87%, 74% and 46% of events with at least two V-shaped images. The fact, that both the angular resolution and the hadronness distribution are not much affected by occurance of the V-shaped images, allows one to use the standard analysis methods for the hadronic background rejection. The remaining background can be estimated using also standard analysis methods, in particular using OFF region observed simultanously with the source of interest.

We check also if the reconstruction of the V-shaped events can be improved using the toy model described in Section 4.3. As an example we perform a joint fit to all the images of an V-shaped event by fixing the $E_{0}$ parameter to the estimated energy of the shower. In the minimization procedure $x_{0}$ and $y_{0}$ are given the starting values obtained from the standard stereoscopic reconstruction and are left free to vary. The $H_{1}$ and $f_{e^{-}}$ are treated as nuisance parameters. In Fig. 10 we compare the impact reconstruction of the classical stereoscopic reconstruction with the one obtained from the toy model fitting procedure for V-shaped events.

In the fitting procedure we either fix the arrival direction of the gamma ray to its true value, or let it free in the fitting procedure.

Since the toy model fitting procedure is using prior knowledge of the arrival direction of the shower its resolution of the impact reconstruction is nominally much better (the RMS of the distribution is 34 m compared to 71 m). If the true direction of the shower is not known and let free to vary, in the toy model fitting procedure the obtained RMS of the distribution is 66 m, much closer to the value obtained from the standard reconstruction. Applying to events without V-shaped images a similar procedure, fitting however instead of two separate subshowers a single shower moving along the nominal direction of the source allows us to obtain a similar improvement in the resolution of the impact parameter when exploiting the knowledge of the true source direction. The obtained RMS of the distribution in this case is 77 m for standard reconstruction, 76 m for geometrical fitting when the source position is left free to vary and 31 m when the source position is set to the nominal value.

4.5 Geometrical energy estimation of V-shaped events

For each image that survived the selection criteria described in Section 4.2 we fit the toy model described in Section 4.3 simultaneously to all images. We minimize using MINUIT2333https://root.cern.ch/root/html534/guides/minuit2/Minuit2.pdf the sum of squared distances of all the pixels in all the telescopes surviving the cleaning to the closest (for a given pixel) of the two lines given by the toy model.

To test the method we fix the core location and height of the first interaction to the true values obtained from the MCs leaving free only the energy of the primary gamma ray and the energy division parameter $f_{e^{-}}$. In Fig. 11 we plot the correlation of the energy obtained from such a fitting procedure to the true energy of the event.

While there is a large spread, the correlation is visible between the true and fitted energy. The energy resolution of the method calculated as the standard deviation of the distribution (computed in the range shown in Fig. 11) is 43%.

To test if this alternative method of the energy estimation (purely geometrical, without using the total light yield seen by the telescope) can be exploited to validate the energy scale of Cherenkov telescope we make a distribution of the ratio of the energy estimation from the stereoscopic reconstruction ($E_{est}$) to the one obtained from the above-described fitting procedure ($E_{sep}$). We perform the study on three sets of MCs. The first, reference Set 1 has nominal total light throughput reflectivity values. In the second and third we have artificially modified the light throughput by $\pm 15\%$ to simulate a systematic uncertainty in the light scale of the telescopes. Such a modification can simulate imperfections in the knowledge of the light collection parameters of the telescope (e.g. quantum efficiency of light sensors, reflectivity of the mirrors) or in the transmission of the atmosphere. All three sets are reconstructed using the nominal (Set 1) MCs. This causes a systematic bias in the energy estimation of Set 2 and Set 3. In Fig. 12 we compare the distribution of $E_{est}/E_{sep}$ for the three sets of MCs.

While the shift of the energy scale is partially reflected in the distribution, the effect is very small. The shift of the energy scale by $\pm 15\%$ results in a shift by only $\pm 3\%$. Such a shift is unfortunately too small for practical application of the method, which will be burdened by a number of its own systematic errors and statistical uncertainties (limited by the rate of V-shaped events and residual hadronic background). The fact that the shift of the $E_{est}/E_{sep}$ proxy is smaller than the shift of $E_{est}$ alone is most probably caused by the bias caused by the height of the first interaction discussed in Section 4.1. Consider a case in which the light throughput of the telescope is overestimated resulting in underestimation of $E_{est}$. In such a case mostly the events that start high in the atmosphere and those that have larger GF separation will be lost as they do not produce enough light to survive the trigger and reconstruction. Since the toy model cannot take into accout the trigger bias the average separations at a given energy given by the toy model are underestimated. Therefore, the fitting procedure counteracts it by favouring lower energies. As a result the effect on $E_{est}/E_{sep}$ is smaller than on $E_{est}$ alone.

5 Discussion and conclusions

We investigated the propagation of the $e^{+}e^{-}$ pair produced in the first interaction of a $\gamma$ ray with atmospheric nucleus. In the energy range most interesting for the LST telescopes (30–300 GeV) the deflection of the first $e^{+}e^{-}$ pair produces a shift of the two components by $\sim 0.1^{\circ}$. If the first interaction occurred high in the atmosphere, and the impact factor of the shower is of the order of a few tens of meter the resulting image can have a very clear V-shaped signature. We derived a simple method to identify the V-shaped images candidates. We investigated how the Chimp/MARS analysis chain of CTA reconstructs events which pass selection criteria for V-shaped images. Despite non-regular shape, the angular resolution as well as energy resolution are similar to the ones obtained from regular, elliptical images. The gamma/hadron separation parameter is slightly shifted to more hadron-like values in case at least two images show V-shaped characteristics. The toy model fitting procedure can be used for more precise estimation of the impact parameter of V-shaped events, however the improvement is mostly related to the assumed nominal position of the gamma-ray source. Similar improvement can be also achieved by applying a simplified geometrical fitting method to classical elliptical images for point-like sources by exploiting the knownledge of the nominal source position. The toy model fitting procedure can be used to perform an alternative estimate of the energy of V-shaped events, not dependent on the total light yield of the telescopes. We investigated if the comparison of such an estimate of the energy with the classical one obtained from stereoscopic reconstruction of the shower can be used to validate the energy scale of the Cherenkov telescope. The change of the telescope reflectivity leaves an imprint on the distribution of the ratio of both energy estimations. However, in the presented analysis the effect has too small magnitude (most probably due to strong dependence of the deflection on the height of the first interaction) that limits the practical usage of the method in the presented here form.

The performance of the proposed here method might be different at larger zenith angles. The first interaction would happen at higher altitude, resulting in longer path of deflection of $e^{+}$ and $e^{-}$ before the development of their subcascades. This would in turn produce more separated sub-images. The observations at higher zenith angle however result in stronger absorption of light in the atmosphere and effectively higher threshold. This would enhance the observation bias on the distribution of the height of the first interaction seen in Fig. 3 worsening the performance of the method. Similarly, the MST telescopes can be used together with LSTs to reconstruct the two subshowers in a more precise way. However, it should be noted that since in V-shaped images the $e^{+}$ and $e^{-}$ subcascades are separated on the camera, the energy threshold for observation of such events should be about a factor two larger then nominal energy threshold for classical elliptical images. As the most pronounced V-shaped events occured for LSTs for energies about 100 GeV, they would be at or below the analysis threshold for MSTs. Finaly, the performance of the method might be improved by using a more refined image fitting algorithm, possibly using templates of the two subshowers (see e.g. de Naurois & Rolland, 2009; Parsons $\&$ Hinton, 2014).

Acknowledgements

This work is supported by the grant through the Polish Narodowe Centrum Nauki No. 2015/19/D/ST9/00616. DS is supported by the Narodowe Centrum Nauki grant No. UMO-2016/22/M/ST9/00583. We would like to thank CTA Consortium and MAGIC Collaboration for allowing us to use their software.

This paper has gone through internal review by the CTA Consortium. In particular we would like to thank Abelardo Moralejo, Konrad Bernlöhr, Bruno Khélifi, Jose Luis Contreras and Johan Bregeon for helpful discussions.