Using Muon Rings for the Calibration of the Cherenkov Telescope Array: A   Systematic Review of the Method and its Potential Accuracy

Markus Gaug; Steven Fegan; Alison Mitchell; Maria-Concetta; Maccarone; Teresa Mineo; Akira Okumura

arXiv:1907.04375·astro-ph.IM·July 11, 2019

Using Muon Rings for the Calibration of the Cherenkov Telescope Array: A Systematic Review of the Method and its Potential Accuracy

Markus Gaug, Steven Fegan, Alison Mitchell, Maria-Concetta, Maccarone, Teresa Mineo, Akira Okumura

PDF

TL;DR

This paper reviews the muon ring calibration method for Cherenkov telescopes, evaluates its potential accuracy for the CTA, and proposes improvements to address systematic effects and enhance calibration precision.

Contribution

It provides a comprehensive review of muon ring calibration, identifies new systematic effects, and suggests analytical methods for improved accuracy in CTA applications.

Findings

01

Identified additional systematic effects affecting calibration accuracy.

02

Proposed new analytical methods to incorporate systematic effects.

03

Estimated muon data rates for calibration and monitoring tasks.

Abstract

The analysis of ring images produced by muons in an Imaging Atmospheric Cherenkov Telescope (IACT) provides a powerful and precise method to calibrate the IACT optical throughput and monitor its optical point-spread function (PSF). First proposed by the Whipple collaboration in the early 90's, this method has been refined by the so-called second generation of IACT experiments: H.E.S.S., MAGIC and VERITAS. We review here the progress made with these instruments and investigate the applicability of the method as the primary throughput calibration method for the different telescope types forming the future Cherenkov Telescope Array (CTA). We find several additional systematic effects not yet taken into account by previous authors and propose several new analytical methods to include these in the analysis. Slight modifications in hardware and analysis need to be made to ensure that such a…

Tables7

Table 1. Table 1: Current CTA telescope designs (Canestrari et al., 2013 ; Aguilar et al., 2014 ; Ambrosi, 2014 ; Inome et al., 2014 ; Schlenstedt, 2014 ; Catalano et al., 2014 ; Glicenstein et al., 2014 ; Puehlhofer et al., 2015 ; Rulten et al., 2016 ; Glicenstein and Shayduk, M. for the NectarCam Collaboration and the CTA Consortium, 2017 ; Mazin et al., 2017 ; Sol et al., 2017 ; White, Richard and Schoorlemmer, Harm, 2017 ; Samarai et al., 2017 ) . The geometrical mirror area has been corrected for the effect of shadowing by the camera or the secondary mirror for the case of the Schwarzschild-Couder design for on-axis incidence of light, but not for camera support structures and other material elements, like e.g. baffles. ∗ Averaged over all azimuth angles, individual values vary around this number because of the hexagonal geometry of the tessellated mirror facets. + flat-to-flat distance of a hexagonal light funnel in front of the photon detector without taking into account the thickness of the plastic and the intended gap between funnels, otherwise 0.002 ∘ –0.003 ∘ less . ∗ ∗ Flat-to-flat distance of a square shaped photon detector without taking into account the ∼ 0.2 similar-to absent 0.2 \sim 0.2 mm gap between pixels . † Two MST camera designs are currently developed which differ by 0.2 ∘ in their respective camera fields-of-view. ‡ The MST optics have been modified with respect to the “classical” Davies-Cotton design such that the radius of curvature of the primary reflector does not correspond exactly to the focal length of the individual mirror facets, but instead to 1.2 times their focal length. This setup improves the timing performance of the reflector with respect to the “classical design”.

		Geometrical	Number	Reflector		Camera		Number	Pixel
	Optics Type	Mirror Area	(Primary)	Aperture		FOV		Camera	FOV
		$A_{geom}$	Mirror Facets	$2 R$		$Ω$		Pixels	$ω$
		(m²)	$N_{mirrors}$	(m)		(deg)		$N_{pixels}$	(deg)
LST	Tessellated Parabolic	378	198	23.	6^∗	4.	5	1855	0.	10⁺
MST	Modified Davies-Cotton^‡	88.5	86	12.	5^∗	7.	8^†	$\sim$ 1800^†	0.	175⁺
SST (ASTRI)	Schwarzschild-Couder	9.2	18	4.	1^∗	10.	5	2368	0.	19^∗∗
SST (GCT)	Schwarzschild-Couder	9.4	6	4.	0	8.	3	2048	0.	15^∗∗
SST-1M	Davies-Cotton	9.4	18	3.	8^∗	9.	1	1296	0.	24⁺

Table 2. Table 2: Published resolutions and accuracies of the estimated parameters impact distance and optical bandwidth , using the different methods presented.

Method	${\hat{ρ}}_{R}$		${\hat{B}}_{μ}$		Reference
	Resolution	Accuracy	Resolution	Accuracy
	(%)	(%)	(%)	(%)
One-dimensional fits
Histogramming	7	n.a.	100	3–12	Jiang et al. (1993); Rose (1995) (Whipple) Shayduk et al. (2003); Goebel et al. (2005) (MAGIC)
Smoothing and renormalizing	2	n.a.	5	2	Bolz (2004) (H.E.S.S.)
cog-based calculations	n.a.	n.a.	9	n.a.	Fegan and Vassiliev (2007) (VERITAS)
Two-dimensional fits
	n.a.	n.a.	15	10	Rovero et al. (1996) (Whipple)
	n.a.	n.a.	5		Chalme-Calvet et al. (2014) (H.E.S.S.)
One-dimensional fits including p.e. statistics
Minimize $χ^{2}$ , Eq. 40	8–13	$>$ 3	50	n.a.	Iacoucci (1998) (CAT) Leroy (2004) (H.E.S.S.)
Two-dimensional fits including p.e. statistics
Minimize $- 2 \ln ℒ$ , Eq. 41	n.a.	n.a.	16	tbd	Mitchell (2016) (H.E.S.S.)

Table 3. Table 3: Approximate molecular transmission and optical bandwidths for observed Cherenkov light from muons for the different telescope types. The first column shows the range of possible maximum emission heights (depending on the muon impact distance, for full rings), the second column the molecular transmission for an averaged Cherenkov light ray from muons at an impact distance of ρ R = 0.5 subscript 𝜌 𝑅 0.5 \rho_{R}=0.5 , and the last column the obtained optical bandwidths for a generic telescope. The uncertainties stem from varying the impact distance of the muons. All values have been derived for an observatory altitude of 2200 m a.s.l.

Telescope	$h_{\max}$	$t_{μ, mol} (300 nm)$	$B_{μ}$
	(m)		(eV)
LST (Hamamatsu PMT)	[550,1100]	0.968	0.	6348 $\pm$ 0.0006
MST (Hamamatsu PMT)	[280,570]	0.982	0.	6401 $\pm$ 0.0006
SST	[95,190]	0.997	0.	6406 $\pm$ 0.0002

Table 4. Table 4: Aerosol transmission t μ , aer subscript 𝑡 𝜇 aer t_{\mathrm{\mu,aer}} at 300 nm wavelength for an averaged Cherenkov light ray from muons at an impact distance of ρ R = 0.5 subscript 𝜌 𝑅 0.5 \rho_{R}=0.5 , for the different telescope types and average and extreme aerosol conditions (see text for details). The obtained optical bandwidth B μ subscript 𝐵 𝜇 B_{\mu} has been calculated for a generic telescope and average molecular transmission. All values have an uncertainty of Δ B μ ≲ 0.001 less-than-or-similar-to Δ subscript 𝐵 𝜇 0.001 \Delta B_{\mu}\lesssim 0.001 due to the residual impact distance dependency. Those entries which deviate by more than 2% from the average are marked as bold. All values have been derived for an observatory altitude of 2200 m a.s.l.

Telescope	$t_{μ, aer}$ (300 nm)	$t_{μ, aer}$ (300 nm)	$t_{μ, aer}$ (300 nm)	$t_{μ, aer}$ (300 nm)
	Average	Extremely	Extremely	Extremely
	(1)	Low	High (no calima)	High (calima)
LST (PMT)	0.977	0.994	0.897	0.968
MST (PMT)	0.987	0.997	0.938	0.984
SST (SiPM)	0.998	0.999	0.988	0.998
	$B_{μ}$	$B_{μ}$	$B_{μ}$	$B_{μ}$
	average	extremely	extremely	extremely
		low	high (no calima)	high (calima)
	(eV)	(eV)	(eV)	(eV)
LST (PMT)	0.624	0.632	0.592	0.614
MST (PMT)	0.634	0.638	0.614	0.630
SST (SiPM)	0.639	0.640	0.636	0.639

Table 5. Table 5: Requirements for the spot size encircling 80% of the light reflected on to the camera ( θ 80 subscript 𝜃 80 \theta_{\mathrm{80}} ). Note that the rms point spread Δ Δ \Delta (Equations 75 and 76 ) encircles only 39.3% of the light in case of a 2-dimensional normal distribution of light. In that case θ 80 ≈ 1.8 ⋅ Δ subscript 𝜃 80 ⋅ 1.8 Δ \theta_{\mathrm{80}}\approx 1.8\cdot\Delta .

Telescope	Requirement	Application Range
Type	on $θ_{80}$	(from Camera Center)
LST	$< {0.11}^{\circ}$	$< {1.2}^{\circ}$
MST	$< {0.18}^{\circ}$	$< {2.8}^{\circ}$
SST	$< {0.25}^{\circ}$	$< {3.2}^{\circ}$

Table 6. Table 6: List of maximally achievable rates of useful muon images for calibration for the proposed telescopes of CTA and corresponding monitoring time scales (see text for details).

	LST		MST		SST
	$ϑ = 0^{\circ}$	$ϑ = 60^{\circ}$	$ϑ = 0^{\circ}$	$ϑ = 60^{\circ}$	$ϑ = 0^{\circ}$	$ϑ = 60^{\circ}$
Maximum impact distance $ρ_{\max}$ (m)	11.8		4.9		2.0
Maximum muon angle $υ_{\max}$ (deg)	0.8		2.5		3.1
Determination of optical bandwidth $B_{μ}$
$R_{\max}$ (Hz)	1.8	0.5	5	1.4	0.8	0.2
$T_{monitor}$ for 1% resolution	14 s	48 s	5 s	18 s	28 min	18 min
Determination of optical bandwidth (per pixel) $B_{pix}$
$R_{\max}$ (Hz)	0.07	0.02	0.10	0.03	0.02	0.006
$T_{monitor}$ for 3% resolution	6 min	21 min	4 min	13 min	2.3 h	7.6 h
Determination of time offset (per pixel) $Δ t_{pix}$
$R_{\max}$ (Hz)	0.1	0.03	0.2	0.04	0.03	0.008
$T_{monitor}$	10 s	36 s	6 s	23 s	36 s	2.2 min
Determination of optical PSF
$R_{\max}$ (Hz)	5	1.6	8	2.5	2.4	0.7
$T_{monitor}$ for 5% resolution	9 s	25 s	5 s	16 s	17 s	1 min

Table 7. Table 7: Summary of systematic uncertainties in percent for CTA. The bold numbers correspond to those which are sizeably larger for CTA than for traditional IACTs, due to the large field-of-view, or new technology employed. The uncertainties marked in blue require additional treatment beyond the analysis of Vacanti et al. ( 1994 ) to be applicable. The numbers in light gray do not apply for standard operation of a given telescope.

Item	LSTs	MSTs	SSTs	Comments
Instrumental part $U_{0}$	%	%	%
1 Determination of average $R$	$<$ 0.5	$<$ 0.5	$<$ 0.5	Account for hexagonal mirrors (Fig. 7), can be calibrated with MC simulations.
Reconstr. Cherenkov angle $θ_{c}$
2 All ring broadening effects	$<$ 0.5	$<$ 0.3	$<$ 0.1	Residual biases.
3 Optical aberration	$<$ 2	$<$ 1	$<$ 1	Assuming plate-scale calibration, otherwise: 4% (LST), 2% (MST and SST).
4 Finite camera focus	$<$ 0.5	$<$ 0.5	$<$ 0.5	Assuming offline correction, otherwise: 8% (LST), 4% (MST), 2% (SST).
5 Bending in geomagnetic field	$<$ 1	$<$ 0.5	$<$ 0.01	Assuming averageing pointings parallel and perpendicular to geomagn. field, otherwise: 2% (LST), 1% (MST).
6 Cherenkov light emission between first and secondary mirror	0	0	$< 0.5$	Assume a suitable correction or cut, outlined in Sect. 5.1.10.
Reconstr. impact distance and $E_{0} (ρ_{R})$
7 Neglect muon inclination	$<$ 0.1	$< 0.1$	$< 0.2$
8 Overall shadow	$<$ 1	$<$ 1	$< 𝟑$	Requires MC calibration using simulations.
9 Changing shadows with inclination angle	$< 𝟏$	$< 𝟐$	$< 𝟏$	Assume careful treatment of effect of displaced quadratic camera shape (Sect. 5.3.2)
Atmospheric transmission $T$
10 Molecular part $T_{mol}$	$<$ 0.1	$<$ 0.1	$<$ 0.04	Assuming correction (Sect. 5.2.2), otherwise $<$ 3% (LST), $<$ 2% (MST), $<$ 0.3% (SST).
11 Aerosol part $T_{aer}$	$<$ 2	$<$ 1	$<$ 0.1	Assume exclusion of bad nights (Sect. 5.2.2), or correction from atmospheric monitoring data.
Reconstructed image size $N_{tot}$
12 Mono trigger biases	0	$<$ 1	$<$ 2	Stereo trigger assumed for the LST, but can be corrected using mono runs.
13 Stereo trigger bias	$<$ 2	$<$ 2	$<$ 5	Assuming that a bias calibration using mono runs is carried out regularly, otherwise $3 - 10$ %.
14 Pulse leakage	$<$ 1	$<$ 1	$<$ 2	Requires signal estimators of sufficient length.
15 Signal estimator biases	$<$ 0.1	$<$ 0.1	$<$ 0.1	Requires un-biased signal estimators, otherwise (8–15)% (Sect. 3.1 and 5.3.1)
16 Image selection biases	$<$ 0.1	$<$ 0.1	$<$ 3	Requires un-biased analysis cuts and image cleaning algorithms optimized for muons.
17 Pixel baselines	$<$ 0.1	$<$ 0.2	$<$ 1	Assuming the baselines are controlled to $<$ 0.3 p.e., as required for CTA
18 Non-active pixels	$<$ 0.5	$<$ 0.5	$<$ 0.5	Requires advanced impact distance fits and less than 10% broken pixels.
Translation muon to gamma efficiency $ε_{μ} \to ε_{γ}$
19 Chromaticity of degradation	$<$ 2	$<$ 1.5	$<$ 1	Assume spectral cutoff by camera protection window, additional monitoring of chromaticity of degradation otherwise: (10–15)% for PMT-based cameras.
20 Dependencies on camera pixel incidence angles	$<$ 0.5	$<$ 0.5	$<$ 1	Assume corrections (Section 5.1.9), otherwise 3–5%.
21 Mis-focused mirrors	$<$ 0.5	$<$ 1	0	Assume only one mirror mis-aligned.
Total	$<$ 4–5	$<$ 3–4	$<$ 5–6	Assuming un-correlated uncertainties
(only items 9–20:)	$<$ 4	$<$ 3.5	$<$ 5	(Assuming previous comparison with simulations)

Equations301

\frac{d ^{2} N _{c}}{d x d ϕ} = \frac{α}{h c} \cdot \int_{ϵ_{1}}^{ϵ_{2}} (1 - \frac{1}{β _{μ}^{2} ( x ) n ^{2} ( ϵ , x )}) d ϵ,

\frac{d ^{2} N _{c}}{d x d ϕ} = \frac{α}{h c} \cdot \int_{ϵ_{1}}^{ϵ_{2}} (1 - \frac{1}{β _{μ}^{2} ( x ) n ^{2} ( ϵ , x )}) d ϵ,

cos θ_{c} (x, ϵ) = \frac{1}{β _{μ} ( x ) \cdot n ( ϵ , x )} .

cos θ_{c} (x, ϵ) = \frac{1}{β _{μ} ( x ) \cdot n ( ϵ , x )} .

E_{t} (x, ϵ) = \frac{m _{μ} \cdot c ^{2}}{1 - 1/ n ( ϵ , x ) ^{2}},

E_{t} (x, ϵ) = \frac{m _{μ} \cdot c ^{2}}{1 - 1/ n ( ϵ , x ) ^{2}},

θ_{c}

θ_{c}

with

θ_{\infty}

\frac{d ^{3} N _{c}}{d x d ϵ d ϕ} ≃ \frac{α}{h c} \cdot sin^{2} θ_{c},

\frac{d ^{3} N _{c}}{d x d ϵ d ϕ} ≃ \frac{α}{h c} \cdot sin^{2} θ_{c},

\frac{d N _{pe}}{d ϕ} (ρ_{R}, ϕ_{0})

\frac{d N _{pe}}{d ϕ} (ρ_{R}, ϕ_{0})

= \frac{α}{2 h c} \cdot sin (2 θ_{c}) \cdot D (ρ, ϕ) \cdot Θ (ρ, υ, ϕ) \cdot B_{μ} .

D (ρ, ϕ)

D (ρ, ϕ)

and ρ_{R} := ρ / R,

\displaystyle\qquad\,\quad\vec{\rho}~{}=~{}\begin{pmatrix}\rho_{x}\\ \rho_{y}\end{pmatrix}~{} =~{}\rho_{R}\cdot\begin{pmatrix}\cos(\phi_{0}{\color[rgb]{0,0,0}+\pi})\\ \sin(\phi_{0}{\color[rgb]{0,0,0}+\pi})\end{pmatrix}

υ = (υ_{x} υ_{y}) = υ \cdot (cos ψ sin ψ) .

Θ (ρ, υ, ϕ)

Θ (ρ, υ, ϕ)

B_{μ}

B_{μ}

ξ_{det} (ϵ) = ξ_{mirr} (ϵ) \cdot ξ_{window} (ϵ) \cdot ξ_{concentrator} (ϵ) \cdot ξ_{pde} (ϵ),

ξ_{det} (ϵ) = ξ_{mirr} (ϵ) \cdot ξ_{window} (ϵ) \cdot ξ_{concentrator} (ϵ) \cdot ξ_{pde} (ϵ),

B_{γ} (h, d, ϑ_{tel}) = \int t_{γ} (h, d, ϑ_{tel}, ϵ) \cdot ξ_{det} (ϵ) d ϵ .

B_{γ} (h, d, ϑ_{tel}) = \int t_{γ} (h, d, ϑ_{tel}, ϵ) \cdot ξ_{det} (ϵ) d ϵ .

B_{γ} (h, d, ϑ_{tel})

B_{γ} (h, d, ϑ_{tel})

with

C_{μ - γ} (h, d, ϑ_{tel})

{\color[rgb]{0,0,0}|\sin(\phi-\phi_{0})|<1/\rho_{R}\quad\mathrm{and}\quad\cos(\phi-\phi_{0})>0}~{}~{}.

{\color[rgb]{0,0,0}|\sin(\phi-\phi_{0})|<1/\rho_{R}\quad\mathrm{and}\quad\cos(\phi-\phi_{0})>0}~{}~{}.

h_{max} = [R \cdot cot θ_{c}, 2 R \cdot cot θ_{c}] ≃ [R / 2 ϵ, R \cdot 2/ ϵ],

h_{max} = [R \cdot cot θ_{c}, 2 R \cdot cot θ_{c}] ≃ [R / 2 ϵ, R \cdot 2/ ϵ],

Q_{pe} (ϕ ∣ ρ, ϕ_{0}, ω)

Q_{pe} (ϕ ∣ ρ, ϕ_{0}, ω)

≃ 1.0 \times 1 0^{2} \cdot (\frac{B _{μ}}{eV}) \cdot (\frac{ω}{deg}) \cdot (\frac{D ( ρ , ϕ - ϕ _{0} )}{m}),

N_{pix}

N_{pix}

Q_{tot} (θ_{c}, ρ)

Q_{tot} (θ_{c}, ρ)

\displaystyle\qquad\mathrm{with:}~{}\Phi=\left\{\begin{array}[]{lcl}\mathrm{arcsin}(1/\rho_{R})&~{}~{}\mathrm{for:~{}}&\rho_{R}>1\\[5.69046pt] \pi/2&~{}~{}\mathrm{for:}{}&\rho_{R}\leq 1\end{array}\right.

\approx U_{0} \cdot θ_{c} \cdot E_{0} (ρ_{R})

\displaystyle\qquad\mathrm{with:}\left\{\begin{array}[]{lcl}U_{0}&=&R\cdot\frac{\displaystyle\alpha}{\displaystyle\hslash c}\cdot B_{\mu}\\[7.11317pt] E_{0}(\rho_{R})&=&\frac{2}{\pi}\int\limits_{0}^{\Phi}\Theta(\rho,\upsilon)\cdot\sqrt{1-\rho_{R}^{2}\sin^{2}\phi}~{}\mathrm{d}\phi\quad,\end{array}\right.

\frac{δ E ( ρ _{R} )}{E ( ρ _{R} )}

\frac{δ E ( ρ _{R} )}{E ( ρ _{R} )}

= \frac{K ( ρ _{R} ) / E ( ρ _{R} ) - 1}{ρ _{R}} \cdot δ ρ_{R}

= \frac{δ ρ _{R}}{F ( ρ _{R} )}

\frac{δ B _{μ}}{B _{μ}} \approx (\frac{δ Q _{tot}}{Q _{tot}})^{2} + (\frac{δ θ _{c}}{θ _{c}})^{2} + (\frac{δ R}{R})^{2} + (\frac{δ ρ _{R}}{F ( ρ _{R} )})^{2},

\frac{δ B _{μ}}{B _{μ}} \approx (\frac{δ Q _{tot}}{Q _{tot}})^{2} + (\frac{δ θ _{c}}{θ _{c}})^{2} + (\frac{δ R}{R})^{2} + (\frac{δ ρ _{R}}{F ( ρ _{R} )})^{2},

Δ θ_{c} = χ^{2} / N_{tot} .

Δ θ_{c} = χ^{2} / N_{tot} .

Δ ϕ = N_{smooth} \overline{ω} = N_{smooth} \frac{ω}{θ _{c}}

Δ ϕ = N_{smooth} \overline{ω} = N_{smooth} \frac{ω}{θ _{c}}

Q_{pe, i} (ϕ_{i} ∣ θ_{c}, ρ, ϕ_{0}, B_{μ})

Q_{pe, i} (ϕ_{i} ∣ θ_{c}, ρ, ϕ_{0}, B_{μ})

\equiv \frac{1}{4} ϕ \geq (ϕ_{i} - 2 \overline{ω}) \sum ϕ \leq (ϕ_{i} + 2 \overline{ω}) Q_{pe} (ϕ ∣ θ_{c}, ρ, ϕ_{0}, B_{μ})

\frac{∣ d ∣}{θ _{c}} = \frac{( cog _{x} - υ _{x} ) ^{2} + ( cog _{y} - υ _{y} ) ^{2}}{θ _{c}} = \frac{1}{2} \frac{ρ _{R}}{E _{0} ( ρ _{R} )} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Using Muon Rings for the Calibration of the Cherenkov Telescope Array: A Systematic Review of the Method and its Potential Accuracy

M. Gaug

Unitat de Física de les Radiacions, Departament de Física, and CERES-IEEC, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Spain

S. Fegan

Laboratoire Leprince-Ringuet, Ecole Polytechnique, CNRS/IN2P3, F-91128 Palaiseau, France

A.M.W. Mitchell

Physik Institut, Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland

M.C. Maccarone

Istituto di Astrofisica Spaziale e Fisica Cosmica di Palermo, INAF, Via Ugo La Malfa 153, I-90146 Palermo, Italy

T. Mineo

Istituto di Astrofisica Spaziale e Fisica Cosmica di Palermo, INAF, Via Ugo La Malfa 153, I-90146 Palermo, Italy

A. Okumura

Institute for Space–Earth Environmental Research and Kobayashi–Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan

Abstract

The analysis of ring images produced by muons in an Imaging Atmospheric Cherenkov Telescope (IACT) provides a powerful and precise method to calibrate the IACT optical throughput and monitor its optical point-spread function (PSF). First proposed by the Whipple collaboration in the early 90’s, this method has been refined by the so-called second generation of IACT experiments: H.E.S.S., MAGIC and VERITAS. We review here the progress made with these instruments and investigate the applicability of the method as the primary throughput calibration method for the different telescope types forming the future Cherenkov Telescope Array (CTA). We find several additional systematic effects not yet taken into account by previous authors and propose several new analytical methods to include these in the analysis. Slight modifications in hardware and analysis need to be made to ensure that such a calibration works as accurately as required for the CTA. We derive analytic estimates for the expected muon data rates for optical throughput calibration, camera pixel flat-fielding and monitoring of the optical PSF. The achievable statistical and systematic uncertainties of the method are also assessed.

gamma rays: general, methods: data analysis, astroparticle physics, atmospheric effects, telescopes

††journal: ApJS

1 Introduction

Muon ring calibration for Imaging Atmospheric Cherenkov Telescopes (IACTs) was first suggested by Hillas and Patterson (1990) and Rowell et al. (1991) and further elaborated in detail by Vacanti et al. (1994). In addition to the Whipple telescope (Fleury and Artemis-Whipple Collaboration, 1991; Jiang et al., 1993; Rose, 1995; Rovero et al., 1996), the method is used to calibrate the optical throughput of practically all currently operating IACTs, including HEGRA (Pühlhofer et al., 2003), H.E.S.S. (Guy, 2003; Leroy et al., 2003; Chalme-Calvet et al., 2014) and MAGIC (Shayduk et al., 2003; Goebel et al., 2005; Meyer et al., 2005), and used as a cross-calibration method for VERITAS (Humensky et al., 2005; Hanna, 2008).

Muon rings were also used to calibrate the optical point-spread function (PSF) of the MAGIC telescopes (Goebel et al., 2005; Meyer et al., 2005). In-depth studies of the muon method have been carried out in several PhD theses, dedicated to the optical throughput calibration of the Whipple telescope (Vacanti, 1993), the HEGRA telescopes (Bolz, 1999), the H.E.S.S. telescopes (Bolz, 2004; Leroy, 2004; Mitchell, 2016), the CAT telescope (Iacoucci, 1998) and the monitoring of the optical PSF of MAGIC (Garczarczyk, 2006). Detailed muon calibration procedures have been described in internal reports of the MAGIC (Meyer, 2005) and VERITAS (Fegan and Vassiliev, 2007) collaborations.

The Cherenkov Telescope Array (CTA) (Actis et al., 2011; Acharya et al., 2013) is the next generation gamma-ray observatory comprising more than one hundred IACTs distributed over two sites (CTA-South at Cerro Armazones near Paranal, Chile, and CTA-North at “Roque de los Muchachos” Observatory (ORM) on the Canary Island of La Palma, Spain, both at around 2200 m a.s.l.) that will provide significant performance enhancement with respect to current facilities (Bernlöhr et al., 2013)111see also https://www.cta-observatory.org/science/cta-performance/.

CTA aims to increase the number of known very high energy (VHE; $E>10~{}\mathrm{GeV}$ ) gamma-ray sources by a factor of ten and to make measurements of their spectral and morphological properties with unprecedented precision (The Cherenkov Telescope Array Consortium: Acharya et al., 2019), for which a well-defined calibration approach, included as an integral part of the CTA design, is paramount (Gaug et al., 2014; Maccarone et al., 2017). In order to meet the required accuracy of 10% for the determination of the global energy scale, the Cherenkov photon transmission through all detector components and subsequent conversion to photoelectrons must be known to better than 5% accuracy, given that the systematic uncertainties due to the imperfect knowledge of the atmospheric conditions cannot be reduced to smaller than 8% within reasonable efforts. Cherenkov light emitted by local muons can serve as a continuously available and well-understood calibration light source, the intensity of which is known to the accuracy with which the Cherenkov angle can be determined (Patrignani, 2017). Because muon calibration can be carried out at the same time as the telescopes perform science observations, the muon calibration method has been chosen as the first option for monitoring of the optical bandwidth of all CTA telescopes (Gaug et al., 2014).

The various telescope types and sizes within CTA (Actis et al., 2011; Acharya et al., 2013) are listed in Table 1 (see The CTA Observatory, 2018)222These include the Large-Sized-Telescope (LST), the Medium-Sized-Telescope (MST) and the Small-Sized-Telescope (SST).. This article aims to address whether the muon calibration method will fulfill the precision and accuracy requirements for each telescope throughout the $>$ 30 years of lifetime of CTA. In this paper, we critically review the method, while in a second follow-up paper, we will compare the performance of the muon calibration between telescopes and across different algorithms with dedicated simulations. Moreover, we investigate the accuracy that muons can achieve for the monitoring of the optical PSF and the camera pixel calibration of each telescope.

We will introduce the main concepts for muon calibration in Section 2. The selection and reconstruction of muon rings are explained in Section 3. The precise reconstruction of their impact distance, ring width, and the telescope’s optical throughput is described in Section 4. Systematic effects are investigated in Sections 5–7, and expected muon image rates are calculated in Section 8. We discuss the findings of the previous sections in Section 9 and draw conclusions in Section 10.

2 Review of the Muon-Calibration Method for IACT’s

Single muons form naturally as part of hadronic air showers and generate Cherenkov light. Most muons detected by IACTs are minimum ionizing particles with a long lifetime.

Muons hence have a large penetration depth, traveling in straight lines toward the ground with Cherenkov radiation emitted under constant opening angle in a cone around the muon direction of travel. As IACTs image in angular space, the Cherenkov light from *local333Local muons are considered those that either hit the telescope mirror directly or come very close to it.

muons traveling parallel (or almost parallel) to the telescope’s optical axis forms a ring-shaped image in the camera. The ring appears complete provided that the muon passes through the primary reflector of the telescope, such that Cherenkov light is reflected from all azimuthal angles around the muon path.

Figure 1 recalls the geometry of the system and introduces the parameters: the impact distance $\rho$ , the inclination angle $\upsilon$ , and the Cherenkov angle $\theta_{\mathrm{c}}$ . The telescope mirror has a radius $R$ . As long as $(\upsilon+\theta_{\mathrm{c}})$ is smaller than the half the camera field-of-view (FOV), the telescope will integrate light along the cord $D$ . The cord itself is a function of $\rho$ and the azimuth angle $\phi$ , defined with respect to a reference angle $\phi_{0}$ , normally chosen as the one at which the cord $D$ is largest.

The number of Cherenkov photons $N_{c}$ emitted in the photon energy range ( $\epsilon_{1},\epsilon_{2}$ ) per unit path length $x$ by a muon in air can then be calculated from the Frank-Tamm formula, provided that the muon energy is above the Cherenkov emission threshold (see, e.g., chapter 34.7.1. of Patrignani, 2017):

[TABLE]

where $\alpha\approx 1/137$ is the fine-structure constant, $h$ the Planck constant, $\beta_{\mu}$ the velocity of the muon, in units of the speed of light in vacuum $c$ , and $n$ the refractive index of the surrounding air. The velocity of the muon, in turn, is related to the Cherenkov angle via

[TABLE]

Requiring $(\cos\theta_{\mathrm{c}}<1)$ leads to a threshold energy $E_{t}$ for the muon to emit Cherenkov light, namely,

[TABLE]

where $m_{\mu}\approx 105.7$ MeV is the muon rest mass and $E_{t}\approx 5$ GeV for a reference altitude of 2200 m a.s.l.

Using the notation $n=1+\varepsilon$ and assuming $\varepsilon\ll 1$ , always valid for air, and $E_{\mu}\gg m_{\mu}c^{2}$ , the Cherenkov angle can be approximated to good accuracy as

[TABLE]

The approximation is better than one in thousand for all cases considered in this paper.

Because muon Cherenkov light observable by CTA telescopes is always emitted from heights of less than 1.5 km above ground, the refractive index can be assumed constant along $x$ to first order, depending only on the observatory altitude. We can then simplify Eq. 1 to the more convenient form (Patrignani, 2017, chap. 35):

[TABLE]

where we have introduced the photon optical bandwidth $\mathrm{d}\epsilon$ . Note that the photon spectrum Eq. 6 shows no explicit dependency on the energy of the Cherenkov photon, in contrast to the conventional parameterization in terms of the wavelength (Vacanti et al., 1994). The spectrum depends only indirectly on the photon energy through its slight dependency of the refractive index, and hence $\sin\theta_{\mathrm{c}}$ .

Integrating over the path length $L=D(\vec{\rho},\phi)/\tan\theta_{c}$ , visible for the camera, we obtain a prediction for the number of observed photoelectrons $N_{\mathrm{pe}}$ , detected by the light sensors of the camera along the ring:

[TABLE]

In the case that the telescope mirror may be approximated by a circle of radius $R$ , the following analytical form for $D(\rho,\phi)$ may be derived:

[TABLE]

Eq. 7 and 10 neglect the influence of the inclination angle on the observable muon track length but are accurate to $O(\tan\theta_{\mathrm{c}}\cdot\tan\upsilon)$ for small inclination angles (see Appendix B). Here we have introduced the possibility of obscuration of Cherenkov light by shadowing elements (such as the camera) in the parameter $\Theta(\vec{\rho},\vec{\upsilon},\phi)$ . In the case of roundish cameras with radius $r$ and small inclination angles, or central holes in the mirror dish, $\Theta(\vec{\rho},\vec{\upsilon},\phi)$ can be approximated as

[TABLE]

Finally, the total optical bandwidth444Note that $B_{\mu}$ is essentially the same as $I$ in Vacanti et al. (1994), although their definition is given in units of nm-1 and the integration limits of 300 nm and 600 nm have been approximated to a normal range of sensitivity of IACT photomultipliers.. of the traversed atmosphere and the detector is

[TABLE]

where $t_{\mathrm{\mu}}(\epsilon)$ denotes the atmospheric transmission, namely, the probability of a Cherenkov photon of energy $\epsilon$ , emitted by the muon, not to get absorbed or scattered out of the reach of the detector. The transmission term $t_{\mathrm{\mu}}(\epsilon)$ has a slight dependency on the impact distance and the telescope pointing zenith angle $\vartheta_{\mathrm{tel}}$ . We will show, however, that the dependency is negligible. Finally, $\xi_{\mathrm{det}}(\epsilon)$ is the probability of detecting (i.e. converting to a photoelectron) a Cherenkov photon that has hit the telescope mirror. Modern IACTs, with a sensitive range from $\sim$ 2 eV ( $\sim$ 600 nm) to $\sim$ 4 eV ( $\sim$ 300 nm), show a bandwidth of about $B_{\mu}\sim 0.65$ eV.

The detector efficiency $\xi_{\mathrm{det}}(\epsilon)$ can be written itself as a combination of mirror reflectance $\xi_{\mathrm{mirr}}$ , camera window transparency $\xi_{\mathrm{window}}$ , light-concentrator efficiency $\xi_{\mathrm{concentrator}}$ and detector photon detection efficiency (PDE) $\xi_{\mathrm{pde}}$ :

[TABLE]

where we have omitted, for simplicity, secondary dependencies of each contribution to $\xi_{\mathrm{det}}(\epsilon)$ on $\rho,\upsilon$ and $\phi$ , which will be treated in later sections of this paper. Figure 2 shows a typical example for several of the mentioned contributions.

Muon calibration then means measuring $B_{\mu}$ , with the best possible resolution and accuracy, from estimators of $N_{\mathrm{pe}}$ , $\theta_{c}$ , $\rho$ and $\phi_{0}$ . We cannot use muon calibration to retrieve any information about $\xi_{\mathrm{det}}(\epsilon)$ , beyond its (weighted) integral. It is, however, possible to determine dependencies of $\xi_{\mathrm{det}}$ on the impact point of a Cherenkov photon of the mirror (through its dependency on $\rho$ ), and the camera pixel (through its reconstruction of the inclination angle), if sufficient statistics have been accumulated.

Many of the gamma-ray parameters are, however, retrieved in an IACT analysis by employing the full shape of $\xi_{\mathrm{det}}(\epsilon)$ as, for instance, in Monte Carlo (MC) simulations of the detector response to Cherenkov light, $B_{\gamma}$ , emitted from gamma-ray-induced air shower particles:

[TABLE]

Here $t_{\gamma}(h,d,\vartheta_{\mathrm{tel}},\epsilon)$ is the atmospheric transmission for Cherenkov light from gamma-ray showers, the properties of which depend on the emission height of the shower particles $h$ , impact distance $d$ and observation zenith angle $\vartheta_{\mathrm{tel}}$ . The detector efficiency $\xi_{\mathrm{det}}(\epsilon)$ is normally measured at or before the construction of an IACT, and particularly its degradation with time requires calibration.

We can formally relate $B_{\mu}$ and $B_{\gamma}$ :

[TABLE]

and need to ensure that $C_{\mu-\gamma}(h,d,\vartheta_{\mathrm{tel}})$ cannot degrade with time beyond a predefined limit, by design of the optical elements of the telescope and for all relevant combinations of $h,d$ and $\vartheta_{\mathrm{tel}}$ 555This part of the muon analysis has often been greatly overlooked in the past, sometimes yielding strongly biased results..

Atmospheric transmission of muon light $t_{\mu}(\epsilon)$ is typically higher than $t_{\gamma}(h,d,\vartheta_{\mathrm{tel}},\epsilon)$ and extends to larger energies (smaller wavelengths), because both the strong extinction of UV light by Rayleigh scattering and absorption by ozone are greatly reduced during the short path traversed by the muon Cherenkov light.

We will hence be obliged to interfere already in the design of the CTA telescopes and cameras to ensure that their possible (chromatic) degradation cannot exceed a predefined limit, as will be shown in Section 6.

Whereas the (effective) mirror radius is determined by the chosen hardware, the impact distance has to be retrieved by the modulation of light intensity around the ring. For full rings, the light intensity reaches a maximum at $\phi=\phi_{0}$ , and is proportional to $(R+\rho)$ , while the minimum is reached at $\phi=\phi_{0}+180^{\circ}$ , where the amount of light scales with $(R-\rho)$ . If the impact distance is zero, the modulation along the ring becomes flat.

It follows immediately from Eq. 10 that for $\rho_{R}<1$ , a full ring is visible in the camera (i.e. all values of $\phi$ are allowed), while otherwise only a part of the ring is imaged within an azimuth angle range, fulfilling:

[TABLE]

Furthermore, a maximum emission height $h_{\mathrm{max}}$ can be derived from which a full muon ring is visible in the camera:

[TABLE]

where the range applies from muons with zero impact distance to those hitting the edge of the mirror.

One can also derive the amount of charge (in photoelectrons) received by one single camera pixel, characterized by its pixel FOV ( $\omega$ ), which covers an azimuth angle range of $\Delta\phi=\omega/\theta_{\mathrm{c}}$ :

[TABLE]

where the last line has been obtained by using a Cherenkov angle of 1.23∘.

Using the pixel FOV, one can roughly estimate the number of pixels hit by the muon light (assuming that the imaged ring width is slimmer than the size of one pixel and that the hit pixels are aligned on the ring with a constant separation angle) and in the following a maximum impact distance to trigger the readout:

[TABLE]

The total number of observed photoelectrons can be obtained by integrating Eq. 7:

[TABLE]

where we have (following the notation of Fegan and Vassiliev, 2007) split the total muon image size into the total muon size detected at zero impact distance ( $U_{0}$ ), the reconstructed Cherenkov angle $\theta_{\mathrm{c}}$ , and a Legendre elliptic integral of the second kind that can be evaluated numerically. In the case of full rings, and neglecting $\Theta(\rho,\upsilon)$ , $E_{0}(\rho_{R})$ becomes $2/\pi$ times the complete elliptic integral of the second kind $E(\rho_{R})$ 666see, e.g., http://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html..

Figure 3 shows the total number of photoelectrons (using $B_{\mu}=0.65$ eV) impinging on the telescope for each of the three telescope types of the CTA, as a function of the muon impact distance, according to Eq. 27. The difference between full and dashed lines shows the effect of the central hole in the primary mirrors, subtracted from the full mirror estimate using again Eq. 27.

For the propagation of uncertainties of the reconstructed impact distance, it is useful to know the derivative of $E(\rho_{R})$ , namely (Whittaker and Watson, 1990, p. 521):

[TABLE]

where $K$ is the complete elliptic integral of the first kind. The function $F$ ranges from $F=2$ at $\rho_{R}=0$ to $F\approx 1.34$ at $\rho_{R}=0.5$ and $F=0$ at $\rho_{R}=1$ .

The total uncertainty of $B_{\mu}$ relates then to the uncertainties, with which the ingredients of Eq. 26 can be reconstructed, as

[TABLE]

where we have not yet considered the uncertainty introduced by shadows.

3 Muon ring reconstruction

Geometrical reconstruction of the muon ring in order to recover the Cherenkov opening angle $\theta_{\mathrm{c}}$ , and the inclination angle $\vec{\upsilon}$ of the muon typically involves careful signal extraction, pixel selection, and event cuts to achieve robust and unbiased results. In this section, we present the main ideas, with the associated limitations and accuracy for the CTA telescopes presented later.

The schemes described below are viable for all CTA telescopes, using regular data taken close to or contemporary with normal science observations, but require improving some of the techniques previously applied. At least for the MSTs and SSTs, an efficient single telescope trigger and muon flagging scheme is presumed to identify candidate local muons.

3.1 FADC pulse integration

IACT cameras typically digitize photon detector signals with fast ADCs (Sitarek et al., 2013; Bartko et al., 2005; Aharonian et al., 2004; VERITAS Collaboration: Hays, E., 2007; Pühlhofer et al., 2014), which allow offline pulse integration and arrival time reconstruction (although readout systems without offline time reconstruction have also been built Barrau et al. (1998); Impiombato et al. (2013)) Algorithms that use sliding search windows or combine the determination of the readout window location with the image reconstruction (Shayduk et al., 2005; Shayduk, M. (2013), for the CTA Consortium) help to improve charge resolution in small $\gamma$ -ray images albeit at the price of a small bias for low pixel signals (Albert et al., 2008) due to night-sky background. As muon ring signals have a low intrinsic time spread of $O(200~{}\mathrm{ps})$ , a fixed pulse integration window used for all pixels around the average arrival time of the muon ring allows for the use of a fixed window pulse integrator without bias and eliminates the contribution of spurious night-sky background to the reconstructed number of photoelectrons $N_{\mathrm{tot}}$ . Such a procedure also reduces the muon image contamination by the parent shower.

3.2 Noise removal

Typically, Cherenkov images are cleaned from spurious signals due to the light of the night sky using a dual-threshold approach (the so-called “tail cut”), whereby all pixels with charge are retained provided that they are above the lower threshold. Additionally, time constraints may be applied, if available (Rissi, 2009; Aleksić et al., 2011; Lombardi et al., 2011; Bird, 2015). Whereas for $\gamma$ -ray images it is desirable to remove low, background-dominated signals, for muon calibration the entire signal needs to be integrated to avoid a biased measurement. The best results for muon calibration have previously been obtained with two passes: first cleaning the image with the traditional procedure to identify whether a ring is present, and then including all pixels with a secure distance from the reconstructed ring (Bolz, 2004; Mitchell, 2016), after fixing the ring parameters.

Only clean rings can then be fit in a robust manner (although promising attempts with wavelet filtering (Lessard et al., 2002) and machine-learning (Bird et al., 2018) methods have been made in the past that do not need to rely on a direct removal of night-sky background noise).

3.3 Preselection

For reasons of processing efficiency and to eliminate muon ring images obviously unsuitable for calibration purposes, candidate images are usually preselected, before fitting the ring (Bolz, 2004; Humensky et al., 2005). In order to remove frequently triggering small-amplitude images from showers or muons of very low energy, the aforementioned “tail cuts” (see section 3.2) are also used to clean the image. Preselection criteria may cut on the number of remaining pixels after the “tail cut” , the mean number of next neighbors $<\!\!\textit{NN}\!\!>$ (a well-focused muon image should contain pixels in a small ring, and hence $<\!\!\textit{NN}\!\!>\approx 2$ ), the number of pixels located at the camera edges (since those are typically affected by inefficiencies and aberration), the relative intensity variation of the remaining pixels (which is typically smaller for muons than for air shower images), the reconstructed Width parameter (Hillas, 1985) (smaller for elongated images and wider for ring-like ones), the Length parameter divided by the image Size (exploiting the fact that the amount of muon light per arc length is constant in zeroth order), or the standard deviation of the time spread in the arrival of light in the pixels, which is normally much smaller for muon rings (Mirzoyan et al., 2006), if compared with shower images.

3.4 Interpolation of non-active pixels

From this moment on, nonactive pixels may bias the ring reconstruction, depending on the algorithms used (Bolz, 2004). The MAGIC Collaboration interpolates the charge and time information of these pixels from their neighboring ones if at least three neighboring pixels are valid (Aleksić et al., 2012). For an optimized muon analysis, one might think of interpolations along the ring, instead of in a circle around the affected pixel.

Groups of broken pixels will more seriously impact the muon calibration results, and their effects must be studied with simulations (see, e.g., Bolz, 2004).

3.5 First-level Ring Fit

After preselection, typically less than 10% of all images remain to which a first muon ring reconstruction can be applied. Several algorithms have been used for this task (Chernov and Ososkov, 1984; Karimäki, 1991; Taubin, 1991; Chaudhuri and Kundu, 1993; Tyler, 2011), all of which calculate the distance of each preselected pixel to a tested ring center, using a candidate ring radius, and minimize its chi-square, possibly being weighted with pixel charges. It is also common to perform two fits in order to reduce biases due to hadronic shower contamination: a first fit to find a candidate ring radius and center, followed by a second fit that excludes pixels too far from the candidate ring. A good initializer for the ring center can be easily derived from the center of gravity of the image exploiting Eq. 36, and a typical Cherenkov angle (e.g. 1.2∘) for the radius. The thickness of the ring may be estimated by (Fegan and Vassiliev, 2007):

[TABLE]

Iterated fits can be made more robust after excluding pixels much farther away from the ring center than a minimum distance and closer to the ring center than a certain maximum distance. The total charge excluded may be used as a quality measure.

3.6 Quality of the Ring Fit

After the fit(s), it is important to reject rings reconstructed with poor quality. Common choices for such quality parameters include the $\chi^{2}/\textrm{NDF}$ of the ring fit, the ring radius (rejecting low-energy muons); and the ring containment.

Further choices may include the number of inactive pixels on the ring (Bolz, 2004), the ratio between image center-of-gravity distance to the camera center and the reconstructed ring radius (Fegan and Vassiliev, 2007) and the completeness of the ring (e.g. the number of nonempty azimuthal bins with a certain, loose threshold in pixel size). Generally, such cuts should be loose in order not to bias the final distribution of $B_{\mu}$ and its statistical properties.

Finally, both Bolz (2004) and Mitchell (2016) suggest the application of a cut on the amount of light found off the ring, which may identify muon rings most likely to be contaminated by background air shower light (see also Section 5.1.2), but do not implement it.

4 Determination of Impact Distance and Optical Bandwidth

Reconstruction of the impact distance $\rho_{R}$ , the azimuthal phase $\phi_{0}$ and the number of photoelectrons along the ring $\mathrm{d}N_{\mathrm{pe}}/\mathrm{d}\phi$ (see Eq. 7) allows one to retrieve the optical bandwidth $B_{\mu}$ on an event-by-event basis, typically on a selection of high-quality events.

4.1 One-dimensional Fits to the Light Distribution along the Ring

Fitting impact distance and optical bandwidth to the distribution of light along the ring, using Eq. 20, has been explored in Jiang et al. (1993); Rose (1995); Shayduk et al. (2003), and Goebel et al. (2005). All authors claim about 3–5% overall statistical uncertainties on their measurements of $B_{\mu}$ (which translates to slightly more than 100% per-event uncertainty, given the large sample sizes used) and around 3%–12% attributed to systematic uncertainties. Some obtain additional systematic biases of up to 30% when comparing MC simulated muons with those from real data. The biases might, however, be related to an incomplete correction of differences between the optical bandwidths of the telescopes to muon and gamma-ray shower light (Eq. 17), or from simulating pure muons instead of muons from hadronic showers. Only Shayduk et al. (2003) provide an estimate of the event-wise precision of the reconstructed impact parameter, namely, $\delta\rho_{R}\approx 0.07$ .

In general, the difficulty here is how to translate from the individual pixels to the distribution of light as a function of $\phi$ , and to make decisions about which pixels belong to the ring and which are considered background dominated instead.

Bolz (2004) and Iacoucci (1998) additionally smoothed the pixel-wise charges over an azimuthal range $\Delta\phi$ of $N_{\mathrm{smooth}}=4$ and $N_{\mathrm{smooth}}=3$ camera pixel widths, respectively (each of 0.17∘ in the case of H.E.S.S.-I):

[TABLE]

The obtained light intensity in an azimuthal bin was then renormalized again to the equivalent of one pixel, e.g. for $N_{\mathrm{smooth}}=4$ (Bolz, 2004):

[TABLE]

Data points $\overline{Q_{\mathrm{pe}}}$ were then constructed in distances of 0.5 $\overline{\omega}$ and fitted to a smoothed version of function Eq. 20.

After applying quality cuts on the results of that fit, a third iteration was performed including all pixels within $2\omega$ of the ring from the uncleaned image to avoid a bias due to the “tail cut.”777See section 3.2 For small Cherenkov angles, the allowed ring width was enlarged, according to the expected broadening of the ring. This procedure achieved an event-wise precision of $\delta\rho_{R}\approx 0.02$ . Bolz (2004) claims 5% per-event statistical and only 2% systematic uncertainty on the reconstructed $B_{\mu}$ 888Generally, the previous analyses used instead $\varepsilon=B_{\mu}/(hc\cdot\int_{\lambda_{1}=300~{}\mathrm{nm}}^{\lambda_{2}=600~{}\mathrm{nm}}1/\lambda^{2}\mathrm{d}\lambda)$ with the somewhat misleading name of “muon efficiency.”.

4.2 Center-of-gravity-based Calculations

An elegant way to estimate $B_{\mu}$ directly from Eq. 20, without the need to pass through (computationally extensive and hence slow) fits, has been outlined in Fegan and Vassiliev (2007).

The impact distance can be related to the distance of the center of gravity ( $\vec{cog}$ ) of the image from the camera center ( $\vec{\upsilon}$ ), and the reconstructed muon ring radius, according to Fegan and Vassiliev (2007):

[TABLE]

After performing the ring fit, the distance between the center of gravity of the image ( $\vec{cog}$ ) and the ring center ( $\vec{\upsilon}$ ) is calculated:

[TABLE]

With the resulting absolute value $|d|$ , Eq. 36 is employed, e.g. using tabulated solutions of $\rho_{R}/E_{0}(\rho_{R})$ , to obtain $\rho_{R}$ . The individual components of the impact distance ( $\rho_{x},\rho_{y}$ ) can then be calculated according to

[TABLE]

The method achieves a 9% statistical per-event uncertainty on the retrieved optical bandwidth. An estimate of the precision with which $\rho_{R}$ could be determined has not been provided.

4.3 Two-dimensional Fits Including Ring-broadening Effects

It was first realized by Rovero et al. (1996) that ring-broadening effects, implicitly integrated in Eq. 20, can be made explicit, and it is possible to fit the measured light in each pixel directly. The simplest way to do so leads to a Gaussian model for ring broadening, characterized by the ring width parameter $\sigma_{\theta}$ :

[TABLE]

Rovero et al. (1996) reached a statistical per-event uncertainty of $\delta B_{\mu}\sim$ 15% with a function of this type and claimed a systematic uncertainty of $\sim 10$ %. Chalme-Calvet et al. (2014), using the same function, achieved an overall uncertainty of 5% on $B_{\mu}$ , but did not specify how this number was obtained.

Rovero et al. (1996) did not assume that $\sigma_{\theta}$ was a constant, but rather a function of the distance $r$ of a given mirror to the center of the mirror dish and further ring-broadening effects, depending on $\theta_{\mathrm{c}}$ and hence $\sigma_{\theta}(r,\theta_{\mathrm{c}})$ . One could improve this approach by making $\sigma_{\theta}$ a function of $\rho$ and $\phi$ (along the lines started to be developed in Rovero et al. (1996)) and include here analytical models of ring-broadening effects, such as those treated in Section 7.

4.4 Two-dimensional Fits Including Photoelectron Statistics

Two algorithms have been introduced so far in the literature that include the intrinsic Poissonian fluctuations of observed charges around the expectation value $Q_{\mathrm{pe}}$ :

Leroy (2004) described a template-fitting procedure, based on Iacoucci (1998) and Le Bohec et al. (1998). The method constructs an expectation charge value for each pixel on the ring following Eq. 39 and constructs a $\chi^{2}$ function for the summed squared deviations of the observed charges $(Q^{\mathrm{obs}}_{i})$ from the expectation values $Q_{\mathrm{pe},i}(\theta_{\mathrm{c}},\vec{\upsilon},\rho_{R},\phi_{0},B_{\mu},\sigma_{\theta})$ , divided by the expected variances, assuming Poissonian statistics999In Leroy (2004) the variance terms were only partly divided between $Q_{\mathrm{obs}}$ and $Q_{\mathrm{pe}}$ , which is rectified here.:

[TABLE]

where $\sigma_{\mathrm{ped},i}$ represents the baseline fluctuations of pixel $i$ (mainly due to night-sky background) and $F_{i}$ the excess noise factor of the photon detector (pixel) $i$ , respectively (for the case of a photomultiplier tube, see Lewis et al. (1987); Bencheikh et al. (1992); Mirzoyan and Lorentz (1997)). The $\chi^{2}$ depends here on seven parameters and must be minimized numerically, although that not all parameters need to be left free in the fit on an event-by-event basis. For instance, some assumptions for the ring-broadening parameter based on the atmospheric conditions and measured PSF could be made. This line of development still awaits being explored in the future.

Leroy (2004) claimed that a fitting procedure minimizing the $\chi^{2}$ , Eq. 40, without taking into account ring-broadening effects, achieves a precision of $\delta\rho_{R}\approx 0.08$ (with a bias of 0.03), whereas Iacoucci (1998) obtained $\delta\rho_{R}\approx 0.13$ and $\sim$ 50% per-event statistical uncertainty on $B_{\mu}$ . The (systematic) uncertainty for $B_{\mu}$ has not been provided; however, Iacoucci (1998) finds an 80% discrepancy between measurement and expectations, attributed primarily to a loss of sensitivity of the telescope below 315 nm.

Mitchell et al. (2016a) instead used a log-likelihood fit of a smoothed image to Eq. 39, with $\rho_{R},\phi_{0},B_{\mu}$ and $\sigma_{\theta}$ as free parameters. The ring center $\vec{\upsilon}$ and radius $\theta_{\mathrm{c}}$ are kept as fixed with the values obtained from the previous ring fit. As previously, the observed distribution $Q_{\mathrm{obs},i}$ is first smoothed along $\phi$ (Eq. 35). The pixel likelihood is calculated from a Gaussian approximation to Poissonian statistics:

[TABLE]

where the summation goes over all pixels included in the (cleaned) image. The assumption of asymptotic Gaussian behavior is, however, only reasonable for pixel charges greater than $\sim$ 10 photoelectrons, which is often not the case, particularly for pixels that do not lie directly on the ring.

A version of the pixel likelihood that includes Poissonian statistics but assumes a Gaussian pedestal and single photoelectron PMT charge distribution can be written instead as follows:

[TABLE]

Table 2 compares the published results in terms of resolution and accuracy of the estimated parameters impact distance and optical bandwidth. One can immediately see that there is a large spread in the resolution obtained, and to a lesser extent in accuracy, for the estimated optical bandwidth ( $\hat{B}_{\mu}$ ), even within a similar method. This provides a hint that the fitting routines themselves seem to be less important for achieving precise and robust results than a careful ring reconstruction and event selection.

4.5 Quality cuts on the Impact Parameter Fit Results

Basically all analyses that use Eq. 20, or a more sophisticated version of it, cut on the $\chi^{2}/N_{\mathrm{pix}}$ of the fit and on an allowed range of the obtained impact distance. Low-impact distances are often excluded because of the shadowing of the camera (which is difficult to model) or holes in the center of the mirror dish, which are common in this field to host instrumentation for pointing (Acero et al., 2010) and camera calibration (Aharonian et al., 2004; Gaug, 2006; Hanna, 2008; Biland et al., 2014). Typical cut values are hence ( $0.2<\rho_{R}<1$ ), where the upper limit excludes the majority of incomplete rings.

Mitchell (2016) also cuts on the reconstructed image width: $0.04^{\circ}<\sigma_{\theta}<0.08^{\circ}$ .

4.6 Treatment of Camera Shadows

A correct treatment of shadows becomes important when the size of the shadow is non-negligible with respect to the size of the primary mirror. This is the case for several CTA cameras, particularly the MST, which hosts a 3 m broad camera to provide the required large FOV. Moreover, that camera is of quadratic shape, complicating further an analytical expression for its shadow (shown in Appendix C).

In case a muon passes either beside or through the camera, a part of the emitted Cherenkov light from above the camera plane will be shadowed. Instead, light emitted below the camera plane is not shadowed by the camera. Its fractional contribution to the total light can ascend to more than 5% for the case of an MST. Both contributions might be affected by additional shadows, like the camera support structure, and the hole in the center of the mirror dish. Figure 4 introduces the relevant contributions and their denomination, used in the following formulae to compute the camera shadow.

We will first derive the muon impact point at the altitude $F$ of the camera $\vec{\rho}^{\,\prime}$ , if the muon is inclined by an angle $\vec{\upsilon}$ (see Eq. 12):

[TABLE]

where the $f$ -number $f=F/2R$ has been used.

The chord of the shadow on the mirror can be approximated as

[TABLE]

using $A$ either as the radius of a spherical camera shape, used together with Eq. 10 for $D$ , or as the half side length of a quadratic camera, used together with Eq. C2 for $D$ .

Note that the projected camera shadow on the mirror will be centered on $-2f\cdot\tan\upsilon\cdot(\cos\psi,\sin\psi)$ and might not completely cover the central hole anymore. In that case, the residual contribution of the uncovered part of the hole must be added to the shadow.

Eq. 47 does not contain any correction for the Cherenkov light $F\cdot\tan\theta_{\mathrm{c}}$ , emitted below the camera, because its contribution to $D_{\mathrm{shadow}}$ cancels exactly with the additional shadow projected on the mirror dish from a height $F$ under the Cherenkov angle.

In the case that the central hole is completely covered by the shadow of the camera, its contribution can be estimated as

[TABLE]

4.7 Modifications for Non-radially Symmetric Mirror Geometries

The relation Eq. 10 for the chord $D(\rho,\phi-\phi_{0})$ assumes a circle to represent the shape of the mirror dish. Real mirrors, however, differ from exact circular geometries, some of them considerably, such as e.g. the H.E.S.S.-II telescope (Cornils et al., 2005). Mitchell (2016) has modified the algorithm such that the mirror dish edge is described by a series of points with an interpolation procedure used to find the intersection of the chord $D(\rho,\phi-\phi_{0})$ with the edge of the mirror. The effect of this change is shown in Figure 7.

4.8 Modifications for Dual-mirror Telescopes

The effect of a secondary mirror may be treated as a large shadow, as no light from the muon can pass through the mirror (see Figure 8).This approach, however, neglects the additional condition that the light must reflect off both mirrors in order to reach the camera. In predicting the amount of light arriving into a single pixel on the focal plane, the full path length at a given azimuthal angle due to both mirrors needs to be taken into account.

Preliminary studies on simulations of dual-mirror telescopes found that the sum of the chords ( $D_{p}+D_{s}$ ) provided a better description of the Cherenkov light distribution in the camera image from a muon, resulting in a more Gaussian distribution of reconstructed efficiencies and more linear behavior of the correction factor with optical degradation with respect to either neglecting the secondary mirror contribution to the chord or treating the secondary mirror as equivalent to a large hole (causing only loss of light; Mitchell (2016)).

In this case, however, to avoid overestimating the total number of photons, a coefficient $\beta$ is introduced to account for the reduction in light due to shadowing by the secondary mirror before the light reaches the primary mirror.

A precise analytical expression for $D(\rho,\phi)$ in the case of dual-mirror telescopes remains to be fully verified by means of a dedicated study. Investigations of the true light distribution around the muon ring using detailed ray-tracing would be necessary to verify both the dependence of the chord across the secondary mirror on $\phi$ and the correct normalization of the light yield.

We note also that the shape of the shadow caused by the secondary mirror may change with increasing off-axis angle; the significance of the impact this may have on the muon calibration remains for further study.

4.9 Monitoring of the Optical PSF

Fits to the equations of type Eq. 39 allow a precise estimator to be retrieved for the optical PSF size of the reflector, namely, the muon ring width $\sigma_{\theta}$ . Typically, such a parameter shows a dependency on the muon energy and hence the ring radius; however, this can be easily calibrated (see, e.g., Bretz, 2006). The distribution of such “calibrated ring widths” can then be compared with those from simulated muon images of a detector with different simulated PSFs (see Fig. 9 and Goebel et al. (2005); Garczarczyk (2006). Experience has shown that a low quantile (e.g. 0.2 is used by the MAGIC collaboration) of the distribution should be preferred over its median, in order to make the comparison for those rings which are not so strongly affected by other secondary ring broadening effects.

4.10 Monitoring of Individual Camera Pixels

The muon light intensity seen by one pixel has been approximated by Bolz (2004) as101010Note that Bolz (2004) introduced a (wrong) factor of two in $k_{c}$ and overlooked the additional normalization factor $1+\sigma_{\theta}/2\theta_{\mathrm{c}}$ .. :

[TABLE]

where $\Delta\theta=\theta-\theta_{\mathrm{c}}$ is the radial distance of the pixel from the ring. Equation 53 approximates the loss of angular coverage $\mathrm{d}\phi$ of a pixel as its position comes to lie farther away from the ring (and vice versa) by introducing the factor $k_{c}$ .

Just as for the overall image (Eq. 7), a pixel-wise optical bandwidth parameter $B_{\mathrm{pix}}$ can be derived that should correlate with the flat-fielding constants from the light flasher calibration (Aharonian et al., 2004; Gaug, 2006; Hanna et al., 2010; Gaug et al., 2014; Brown et al., 2015). Bolz (2004) finds a good correlation for the H.E.S.S.-I telescopes with a spread of 5% standard deviation.

This approach enables the identification of issues such as gradients across the camera (see Fig. 10) and differences in photodetector response between vertically incident light (such as used for calibration; Aharonian et al. (2004); Gaug (2006); Hanna (2008) and Biland et al. (2014)) and light arriving at the camera focal plane with a distribution of incident angles.

Although Eq. 53 has been used with some success by Bolz (2004), it is intrinically inconsistent: the infinity at $\Delta\theta=-\theta_{\mathrm{c}}$ leads to a nonconverging normalization integral. An improved model should include at least the non-negligible possibility for light scattered across the ring center into the opposite part of the circle and the finite pixel size at the ring center.

4.11 Monitoring of individual mirror facets

Plotting the optical bandwidth against the muon impact points $\rho_{x,y}$ enables possible influences of the individual mirrors and shadows on the focused reflectance to be identified (see Fig. 10). Whilst shadowing effects are already taken into account by the telescope model in simulations, effects such as the degradation of individual mirror facets are not usually included.

5 Secondary Effects Affecting the Reconstructed Optical Bandwidth

5.1 Instrumental Effects

5.1.1 Effects of Non-active Pixels

Pixels can become broken or be deactivated because of high illumination levels, e.g. due to bright stars (Aharonian et al., 2004; Aleksić et al., 2016). The latter typically lead to small clusters of 2–3 deactivated pixels111111The number of pixels deactivated due to bright stars mainly depends on the pixel FOV., whereas broken pixels tend to be more uniformly spread over the camera. Additionally, groups of pixels forming “clusters” or “drawers” may stop working simultaneously121212In the past, clusters of typically 7 or 16 pixels were affected, in a camera hosting about a thousand pixels. . The effect of both has been investigated by Bolz (2004) and Mitchell (2016) with dedicated MC simulations.

Both find stable results for the reconstructed optical bandwidth up to an average of 10% of broken pixels on the ring. Interestingly, the large number of deactivated pixels, averaged over neighboring signal amplitudes, leads to an overestimation of the optical bandwidth in Bolz (2004), because the reconstructed intensities along the ring become flatter (if the missing pixel intensities are replaced by the average of their neighbors’ intensities, as explained in Section 3.4) and hence the impact parameter is underestimated. Mitchell (2016), who employed the likelihood fit given by Eq. 41, found robust behavior even up to 30% of deactivated pixels. Both studies come to the conclusion that clustering of broken pixels does not alter the behavior of this systematic effect much compared to randomly removed pixels. Mitchell (2016) showed that the effect of broken pixels on the quality cuts applied later is minor, except for the criterion of a required minimum number of pixels contained in the ring image. Nonactive pixels are unlikely to be a problem for muon calibration with CTA, given that all cameras are required to maintain 95% pixel availability during observations.

5.1.2 Effects of the Stereo Trigger

If a telescope array is triggered in stereo mode (Aharonian et al., 2006; Rajotte, 2014; Aleksić et al., 2016), single muons rarely trigger more than one telescope simultaneously. Consequently, the muon events obtained from local muons are often recorded together with the parent air shower, which produces an additional Cherenkov light contribution that may contaminate the muon ring images.

Studies with H.E.S.S. showed that this effect may lead to a positive bias of up to 10% on the reconstructed optical bandwidth from simulations (Mitchell, 2016), whilst only 3% difference was observed when comparing muon rings obtained from a single telescope trigger to those obtained using a stereo trigger (Bolz, 2004). Suggestions for dealing with this effect include simultaneous fits of a Hillas ellipse and a muon ring to the images in order to locate and remove the air shower contribution, although this remains to be demonstrated (Mitchell, 2016).

5.1.3 Effects of the Pixel Baselines

Drifts in the pixel baselines, which are known to vary with time and temperature (Aharonian et al., 2004; Sitarek et al., 2013; Biland et al., 2014), may lead to biasing the reconstructed number of photoelectrons in Cherenkov images. For CTA, the baseline of each pixel is required to be measured with an accuracy of better than 0.2 times the equivalent pulse height of one photoelectron and must therefore be regularly monitored. Within the muon analysis, monitoring of the Off-ring intensity 131313“Off-ring intensity” means here the summed pixel intensity for all pixels excluded from the muon ring reconstruction (Bolz, 2004), which should peak around zero, may help to identify systematic drifts or issues in the hardware, such as a few slightly defocused mirrors, or to detect possible biases in the muon analysis itself.

5.1.4 Effects of the Readout Window

The Cherenkov light pulse received by a camera pixel from a local muon is much shorter $O$ (100–200 ps), with a standard deviation arrival time spread across the camera of $<0.7$ ns, than typical signals from air showers, which may spread for several nanoseconds across the camera and show a typical light pulse width of 0.8–3 ns for $\gamma$ -ray shower images (Mirzoyan et al., 2006). The choice of telescope optics and differences in the transit time spread of the photodetectors (Albert et al., 2008) may further smear the arrival time distribution. The MST telescope structure has been modified to reduce the $O$ (5 ns) time spread of a traditional Davies-Cotton reflector (Oakes et al., 2015), whilst for dual-mirror telescopes the arrival time spread exhibits an angular dependency, ranging from isochronous for on-axis rays to $O$ (1 ns) at $\sim 5^{\circ}$ off-axis (Vassiliev et al., 2007; Rulten et al., 2016). Shifts of up to 3 ns in mean intensity maxima between mono and stereo muons have been observed (Bolz, 2004). It is important to ensure that the readout window fully contains the muon signal. The readout window may, however, be further optimized to reduce contamination by an associated air shower signal.Additionally, the gamma-ray analysis should integrate the full gamma-ray pulse such that an unbiased conversion from the optical bandwidth obtained from muons to those of gamma-rays is obtained. Otherwise, correction factors need to be applied that take into account the difference in pulse coverage.

5.1.5 Optical Aberrations

Aberrations affect the estimated muon parameters, apart from the global broadening of the ring width. In particular, coma and distortion (Schliesser and Mirzoyan, 2005; Vassiliev et al., 2007) can move the image centroid away from its ideal position and hence bias the reconstructed ring width, if not properly taken into account. The apparent radius and center of the ring image on the focal plane are hence not completely independent of the impact point and inclination angle of the muon. For example, rays from an on-axis muon with $\rho_{R}=0$ preferentially impact the mirrors at smaller angles than is the case for a muon with large impact distance producing a different pattern of aberrations for each case.

The tangential image centroid coordinate $<\!\xi\!>$ can be approximated for parallel rays, uniformly distributed over the reflector of diameter $\mathcal{D}$ and focal length $F$ , with an angle $\upphi$ with respect to on-axis incidence, as (Vassiliev et al., 2007; Fegan and Vassiliev, 2018):

[TABLE]

where $\kappa=1$ corresponds to the ideal classic Davies-Cotton design (Davies and Cotton, 1957), $\kappa=1/2$ to the case of a parabolic reflector, and $\kappa=0.83$ to the “modified Davies-Cotton” reflector proposed for the MST (Farnell et al., 2015). A similar approximation can be used for the Schwarzschild-Couder telescope (Vassiliev et al., 2007).

Eq. 54 translates into a $\sim$ 4% shift with respect to the $\upphi$ -scaling for the parabolic LST141414even $\sim$ 7% for the MAGIC telescopes (Garczarczyk, 2006). and a $\lesssim$ 2% deviation for the MST/SST-1M. Corrections could be implemented straightforwardly through an adequate redefinition of the plate scale. Such a “plate-scale calibration” fits the observed image centroid (e.g. of astrometric standard stars) in camera pixel coordinates to the tangential and sagittal angular distances $(\xi,\eta)$ . The distribution of light becomes, however, strongly asymmetric along the tangential coordinate, at least for coma aberration, as one moves to high incidence angles. In such a case, a Gaussian ring width model, as used in Eq. 39, may result in biased ring parameters, as the muon inclination angle increases. Such biases need to be investigated with carefully simulated muon ring images. Cuts on the reconstructed impact parameter and/or muon inclination angle may be necessary to limit this effect.

Moreover, Eq. 54 has been obtained by integrating parallel light rays uniformly covering the reflector, an assumption that does not hold for the conical wavefront of Cherenkov light from local muons, which hit the reflector preferentially toward the muon impact point. The plate-scale correction will hence be different in our case and probably depend on muon impact distance and incidence angle.

We will conservatively assume that the error made for the determination of $\theta_{\mathrm{c}}$ by the (gamma-ray-optimized) plate-scale calibration is less than half of the calibration correction itself, until a dedicated study on this issue using ray-tracing has been carried out. Such a study will be part of our second paper on this topic.

5.1.6 Finite Camera Focuses

Large IACTs have rather limited depths of field and are unable to focus all parts of shower images at the same time (Hofmann, 2001). This applies particularly to the LST and the MST, which will focus their cameras at the mean altitude at which air shower maxima are observed, normally chosen as 10 km distance above the telescope.

For an ideal telescope of focal length $F$ , focused at infinity, i.e. with the camera placed at $z_{f}=F$ , the Cherenkov angle of the emitted rays equals the radius of the imaged ring. For a telescope focused at a distance $x_{f}$ , i.e. with the camera placed at $z_{f}$ with $1/z_{f}=1/F-1/{\color[rgb]{0,0,0}x}_{f}$ , the rays from an on-axis muon radiating at distance $x$ will form a ring of radius

[TABLE]

The values of $x=0$ and ${\color[rgb]{0,0,0}x_{f}}=\infty$ give $\theta_{\textit{xf}}=\theta_{\mathrm{c}}$ as expected. For $x=x_{f}$ we get $\theta_{\textit{xf}}=0$ ; this is, of course, the definition of what it means for the telescope to be focused at $x_{f}$ . At its highest point, the Cherenkov light from a local muon is imaged into a smaller ring with radius smaller than $\theta_{\mathrm{c}}$ , by the time the muon impacts the mirror $(x=0)$ the angle is the nominal value. So this results in a systematic error if the ring radius is directly interpreted as the mean Cherenkov angle. The bias introduced on the reconstructed Cherenkov angle amounts to (see Appendix F)

[TABLE]

For the MST with ${\color[rgb]{0,0,0}x_{f}}=10$ km, the bias ranges from $\sim-2\%$ to $\sim-4$ %, depending on the muon energy, whereas the LST can suffer biases of as high as $-8\%$ .

5.1.7 Shadowing Effects

For the MST, the large 3 m camera and masts cause an average shadowing in excess of 15% of the muon light. The effect of the central shadow (or of a hole in the mirror dish) may be approximated by subtracting from Eq. 7 a chord of the same shape, but using the radius of the hole (or the camera), instead of the mirror. Fig. 3 has been obtained in this way; however, the square shape of the MST camera makes this approach inaccurate (see Section 5.3.2). For dual-mirror SSTs, shadowing in excess of 35% is anticipated mainly owing to the secondary mirror, but also telescope masts (Rulten et al., 2016), although the exact value and its variation (between muons and air showers, and with $\rho$ and $\psi$ ) remain to be studied with dedicated simulations.

5.1.8 Nonuniformities of the Camera Acceptance and Mirror Reflectance

Nonuniformities of the camera acceptance can be controlled by a pixel-wise reconstruction of the optical bandwidth (see Section 4.10). However, acceptance gradients and particularly changes of the dependency of the pixel response with the incidence angle of light may bias the obtained optical throughput correction for gamma-ray shower images. Gradients of up to 3% have been found by Bolz (2004) for the H.E.S.S.-I telescopes, after plotting the mean $B_{\mu}$ against the inclination azimuth angle.

5.1.9 Nontrivial Dependencies of the Optical Bandwidth on Incidence Angle

Whereas $\gamma$ -ray-induced showers may be assumed to illuminate the reflector homogeneously, muons illuminate the reflector preferentially close to the impact point. Together with the variation in efficiency of the photodetectors (PMTs or SiPMs) with incidence angle, this can lead to a nontrivial dependence of the reconstructed $B_{\mu}$ on $\rho_{R}$ .

Photo-detection devices, such as PMTs (Okumura et al., 2017) or SiPMs (Otte et al., 2017), used for IACT cameras exhibit quite pronounced dependencies of their photon detection efficiencies on incidence angle. IACT camera pixels are, moreover, often preceded by light concentrators (or funnels, often approximating as much as possible the form of a compound parabolic concentrator (CPC); Winston et al. (2005) or Bézier curves Okumura (2012)) to reduce dead spaces and enhance overall optical efficiency (Hénault et al., 2013; Aguilar et al., 2015; Hahn et al., 2017; Okumura et al., 2017). Such light concentrators achieve only approximately uniform photon detection efficiencies up to a carefully chosen maximum incidence angle. Residual dependencies of several percent per degree remain and need to be taken into account. Moreover, light concentrators in front of the photodetector plate may degrade with time (unless they are kept in the camera housings and are hence not exposed to the outer environment Okumura et al. (2017)), particularly the dependency of their acceptance with respect to light incidence angle can change.

The dependence of the photon yield as a function of impact point and nonuniformity has been modeled by an analytic approximation (see Appendix E), where a parameter $\beta$ describes the enhancement $(\beta>0)$ or suppression $(\beta<0)$ of the response for photons reflected from the edge of the reflector. This is shown in Figure 11: for a muon with low impact parameter ( $\rho_{R}\approx 0$ ), the difference in total muon light received (and hence the systematic error) scales roughly as $-\beta/3$ , until $\rho_{R}\approx 0.3$ , whereupon the error reduces to zero at $\rho_{R}=1/\sqrt{2}\approx 0.7$ (by construction) and finally reaches $\approx\beta/9$ at $\rho=1$ . Hence, the muon light yield varies by about the same amount as $\beta$ across the possible ranges of $\rho_{R}$ , whereas by construction the telescope’s response to gamma-ray shower light remains the same, and independent of $\beta$ .

Incidence-angle-dependent Photon Detection Efficiencies

Typical values of $\beta$ are $\sim$ 0.07 for the PMT-based cameras, mainly introduced by the light funnels in front of each PMT151515However, a quite pronounced wavelength dependency of this value is expected (see Fig. 10 of Hayashida et al., 2017). , and $\sim$ -0.06 for the solutions involving SiPMs.

Incidence-angle-dependent Transmission of the Protecting Camera Entrance Window

The camera entrance window will introduce further incidence-angle-dependent effects: the amount of light lost to reflections at the two window surfaces changes with incidence angle such as the optical path length through the window. The latter will introduce a decrease in transmissivity with increasing incidence angle and hence a shift of the “cutoff” wavelength of the window transmission spectrum. The conversion of the optical bandwidth to muon Cherenkov light to that of gamma-ray-induced Cherenkov light (see Section 6) becomes then slightly dependent on incidence angle. Simulations show that these effects can cause an additional contribution to $\beta$ of up to $-0.1$ and will be treated in more detail in our second paper.

A flat camera entrance window shifts the light ray in tangential direction, as a function of the the relative distance of the impinging light ray from the camera center $r_{R}$ , the window thickness $d$ and its refractive index $n$ as

[TABLE]

Inserting realistic numbers for $d\lesssim 6$ mm and $n\approx 1.5$ for plexiglas, one can see that the angular shift $\Delta\upphi$ is always smaller than 0.005∘ and hence negligible.

The situation changes slightly with a conic camera protection window, as foreseen for the LST (Ambrosi et al., 2015) and one proposed MST camera solution (Glicenstein, 2016). In that case, the incidence angle of the light ray on a conic window becomes

[TABLE]

where $R_{\mathrm{window}}$ is the window’s curvature radius and $h_{\mathrm{offset}}$ the distance of the window’s center from the focal plane. Inserting realistic values of $R_{\mathrm{window}}\approx 3.4$ m and $h_{\mathrm{offset}}\approx 0.4$ m yields now a maximum shift of up to 0.01∘. Figure 12 shows an example. On average, such an offset will translate, however, into an error of the reconstructed muon Cherenkov angle of considerably less than 1%.

Additionally, deflection of the rays occurs owing to the conic window shape, and the incidence and exit direction of the ray to and from the window are now different. The net effect seems to be negligible, though for all proposed solutions so far.

5.1.10 Cherenkov Light Emitted by the Muon between the Primary and Secondary Mirror

Additional considerations for a Schwarzschild-Couder telescope include the amount of light emitted between the secondary and primary mirrors if the muon travels through the secondary mirror. This contribution amounts to a maximum of about 3% for the dual-mirror SSTs, considering a distance of about 3 m between primary and secondary mirror (Canestrari et al., 2013). However, the small opening angle of the Cherenkov light cone may result in a large proportion of this light being lost through the hole in the primary mirror for muons passing near the center of the dish. This contribution is hence applicable only to muons with medium impact distance, which avoid the hole in the primary mirror but still hit the secondary. A careful analysis might detect that contribution and correct for it; otherwise, a suitable impact distance cut will remove it.

5.1.11 Trigger and Selection Biases

If the mirror area is small, the reflected muon light intensity may not be strong enough to ensure a stable trigger efficiency above a certain muon energy threshold. In that case, only Poissonian upward fluctuations will launch a trigger, and the distribution of reconstructed optical bandwidths $B_{\mu}$ will become asymmetric and its mean value biased. Mirror degradations of reflectance will not scale linearly with $<B_{\mu}>$ anymore, and a bias correction needs to be applied. The corresponding correction factor needs to be retrieved from simulations; however, CTA will probably not be able to continuously monitor reflectance of each single mirror, and simulations will therefore tend to level out real fluctuations among the different mirrors. The correction factors may hence suffer from a considerable systematic uncertainty.

We have tried to simulate the effect using a Gaussian distribution of a certain mean muon image size and its square root as Gaussian width. Then, the lower parts of the distribution were removed until a given trigger efficiency was reached (see Fig. 13 left panel). The difference between the statistical mean of the new distribution and the original image size has then been defined as a toy-model trigger bias (Fig. 13 right panel). This extremely simplified model may give an order of magnitude of the expected real trigger biases.

However, if the distribution of $B_{\mu}$ shows an identifiable cut toward lower values, one can try to fit only the peak of that distribution, instead of calculating its mean or median. If such an approach can help, it needs further investigation. In any case, all CTA telescopes are required to ensure that their triggers are sensitive enough to local muon light, such that they can ensure that the trigger bias becomes negligible, even if the telescopes’ optical throughput degrades within the permitted ranges.

5.1.12 Mirrors

The mirror reflectance degrades with time, most probably in a wavelength-dependent way. The overall average degradation is a free parameter of the calibration procedure and is part of the obtained optical throughput. A general worsening of the PSF of the mirror should be retrieved by a careful analysis of the reconstructed muon ring width.

More difficult is the estimation of the effect of a few very misaligned mirrors, possibly on very short time scales, due to malfunctioning of the Active Mirror Control for that particular set of mirrors. In that case, the muon light gets reflected into the margins of or even outside the ring by those particular mirrors and possibly misinterpreted as off-ring contribution by the muon analysis. In this case, the disentanglement of the worsened PSF from the effects of the mirror reflectance becomes less accurate. For the case of the LST, one misaligned mirror out of $\sim$ 200 (Ambrosi et al., 2014) may create a systematic bias of $\lesssim$ 0.5%, while for the MST one mirror out of 86 (Schlenstedt, 2014) creates a bias of $\lesssim$ 1%.

5.2 Atmospheric Effects

5.2.1 Bending of the Muon Trajectory in the Geomagnetic Field

Bending of the muon trajectory in the geomagnetic field has not been considered important so far for local muons. It is not even mentioned as a source of ring broadening by Vacanti et al. (1994). A muon of momentum $p$ propagating a distance $D$ perpendicular to the field of strength $B$ will be deflected by an angle $\Delta\theta$ (see, e.g., chap. 35.12 of Patrignani, 2017):

[TABLE]

where $\vartheta,\varphi$ are the local muon direction angles and $B_{\perp}$ the magnetic field component perpendicular to its velocity. The maximum achieved at La Palma at a zenith angle of $\sim 40^{\circ}$ is of $B\approx 40~{}\mu$ T (Commichau et al., 2008), whereas at Paranal maximum vertical $B$ -fields of $\sim 25~{}\mu$ T are achieved (Hassan et al., 2017). The angle under which the muon Cherenkov light is imaged in the camera gets then slightly shifted by (see Appendix G)

[TABLE]

The shift, averaged over the ring, has an expectation value $E[\Delta\theta_{\mathrm{c}}]$ of

[TABLE]

The bias for the reconstructed Cherenkov angle $B(\theta_{\mathrm{c}})=(E[\Delta\theta_{\mathrm{c}}])/\theta_{\mathrm{c}}$ can reach up to 2% for the worst case of a full muon ring imaged into an LST camera, under the maximum impact distance and the azimuthal component of the $B_{\perp}$ parallel azimuth component of the impact parameter.

Geomagnetic field effects can be studied straightforwardly by comparing data from pointings toward the south (where the perpendicular magnetic field is strongest) and the north (where it is weakest, almost reaching zero; Commichau et al. (see 2008)) in the case of La Palma.

5.2.2 Atmospheric Transmission

A photon of energy $\epsilon$ emitted at a certain distance $r$ above the telescope suffers molecular and aerosol extinction before reaching the mirror. One can write the atmospheric transmission for that photon very generically as

[TABLE]

where $\alpha_{\mathrm{mol}}$ and $\alpha_{\mathrm{aer}}$ are the volume extinction coefficients from molecular and aerosol extinction at an altitude $h$ , respectively, and the telescope itself points to the sky under a zenith angle $\vartheta$ .

Assuming that the muon light emitted along the track $L$ from $r_{\mathrm{max}}=D(\vec{\rho},\phi)/\tan\theta_{\mathrm{c}}$ to the telescope mirror dish, we can write

[TABLE]

The transmission expressed in Eq. 64 is independent of the azimuthal coordinate of the impact distance $\phi_{0}$ , because of the radial symmetry of the system. It must then be folded with the telescope efficiencies to obtain the optical bandwidth of the full system:

[TABLE]

Molecular Transmission

The volume scattering cross section for molecular extinction in the atmosphere can be written as (McCartney, 1976; Bucholtz, 1995; Tomasi et al., 2005; Gaug, 2014):

[TABLE]

where $P(h)$ and $T(h)$ are the atmospheric pressure and temperature at altitude $h$ , and $P_{s}$ and $T_{s}$ are the reference values of the U.S. standard atmosphere (NASA, 1976).

Assuming a typical winter atmosphere at the ORM (Gaug, 2016), the term $\left(P(h)/P_{s}\right)\cdot\left(T_{s}/T(h)\right)$ can be described by a scale height of $H_{0}\approx 9.7$ km; hence,

[TABLE]

The result using Eq. 68 for $\alpha_{\mathrm{mol}}(\epsilon,h)$ to calculate $t_{\mu,\mathrm{mol}}(\epsilon,\rho_{R};\vartheta)$ (Eq. 65) is displayed in Fig. 14 for three CTA telescope types. Since the LST is most affected by molecular extinction, the effect of different impact distances is also shown. We studied also different telescope observation zenith angles; however, their effect was even smaller, i.e. well below 0.1%, and is hence not displayed. Averaging $B_{\mu}(\rho_{R})$ over the different impact distances yields $<B_{\mu}>=0.6348\pm 0.0006$ for the LST, hence a systematic effect of less than 0.1%.

One can see that only the LST shows non-negligible molecular extinction of muon light; however, its contribution to $B_{\mu}$ can be safely estimated to lie always below 3%, independent of the impact distance of the muon. In order to assess the range of variation of $t_{\mu,\mathrm{mol}}$ due to changes in the atmospheric profile for the LST, we used the data of one full year from the “Global Data Assimilation System” (GDAS) “Final Analysis” database161616https://rda.ucar.edu/datasets/ds083.3/ for the closest grid point to La Palma, whose temperature and pressure predictions show excellent correlation with the measurements of the MAGIC weather station (Gaug et al., 2017), apart from a constant 2∘C bias, which can be corrected. The full height-resolved temperature and pressure profiles were used to derive $B_{\mu}$ for fixed telescope efficiencies $\xi_{\mathrm{det}}(\epsilon)$ for the LST. The resulting distribution of $B_{\mu}$ shows a peak-to-peak variation of less than 0.2%. We can conclude that the uncertainty of the molecular profile on $t_{\mathrm{\mu,mol}}$ is absolutely negligible, at least for CTA-N, once the correct average atmospheric profile has been selected.

Aerosol Transmission

Astronomical sites, like Armazones and La Palma, are characterized by extremely clean environments and small aerosol content close to the ground. Patat et al. (2011) find a median aerosol extinction of 0.045 mag airmass*-1*, and the semi-interquartile range is 0.009 mag airmass*-1* at 400 nm wavelength for the VLT site at Paranal (about 500 m higher than the CTA-S site). The height profile of the aerosol extinction is not provided, however. We can nevertheless assume that stratospheric aerosol accounts for about 0.005 mag, and the rest forms part of the nocturnal boundary layer close to ground. Hence, about $(4\pm 1)$ % of the Cherenkov light from gamma-ray showers is expected to be scattered out of the FOV in the first kilometer on the ground, probably slightly more on the CTA-S site, due to the lower altitude.

At La Palma, García-Gil et al. (2010) measure a median extinction of $0.129\pm 0.002$ mag airmass*-1* in the $V$ band, throughout the whole year, except for the summer months, at 2369 m a.s.l. (about 100–150 m higher than CTA-N). Subtracting molecular extinction of 0.079 mag airmass*-1* and stratospheric ozone absorption of 0.023, an aerosol extinction of about $\approx$ 0.027 mag airmass*-1* is obtained.

The wavelength dependency of the aerosol extinction is typically expressed by the Ångström coefficient $A$ , where

[TABLE]

with $A\approx 1.4$ for Paranal (Patat et al., 2011), and ranging from $A\approx 0.8$ to 2.0 for La Palma (Whittet et al., 1987). Similar values have been obtained for Tenerife by Maring et al. (2000), ranging from $A\approx 1.5$ to 2 during nondusty nights, and Andrews et al. (2011), which tend toward $A\approx 1.2$ for small scattering coefficients.

The MAGIC Collaboration has found that the altitude profile of the aerosol extinction coefficient for clear nights is exponential with a scale height of $H_{\mathrm{aer}}\approx$ (500–700) m (Gaug et al., 2017), with an aerosol optical depth (AOD) on the ground of about (AOD $\approx 0.02$ ) on average, however ranging from practically zero to almost 0.1 for normal clear nights.

Stronger aerosol extinction is possible only during Saharan dust intrusions (called “calima” in the Canary Islands). Such nights show severely enhanced optical depths (up to close to or greater than 1.0), Ångström coefficients close to zero, and constant aerosol extinction from the ground to approximately 5 km altitude (Lombardi et al., 2011).

We have calculated the aerosol transmission for muons $t_{\mu,\mathrm{aer}}$ for an average case (AOD ${}_{\mathrm{532~{}nm}}=0.03$ , $H_{\mathrm{aer}}=600$ m, $A=1.2$ ) and two extreme cases of clear nights with the exponential profile, created to yield strong extinction of muon Cherenkov light (AOD ${}_{\mathrm{532~{}nm}}=0.08$ , $H_{\mathrm{aer}}=500$ m, $A=2.0$ ) and tiny extinction (AOD ${}_{\mathrm{532~{}nm}}=0.01$ , $H_{\mathrm{aer}}=700$ m, $A=1.0$ ), just to show the range within which $t_{\mathrm{\mu,aer}}$ can vary for clean nights. The case of Saharan dust intrusions (“calima”) is treated as another extreme example: AOD $=0.5$ (supposing that this is the absolute limit for observation), $A=0.$ and constant aerosol extinction up to 3 km above the ground. The results are shown in Figure 15 and Table 4.

One can see that only in the case of strong aerosol densities, very close to the ground and with an unusually high contribution of the accumulation mode particles, does a correction of less than 3% need to be applied for the MSTs and LSTs. This problem can be either circumvented either by adequate data selection (e.g., using LIDAR data) or by applying a bias correction, which can be calculated once the aerosol profile is assessed. Such a correction is expected to take place only rarely, when the observatory operates under nonoptimal conditions. Variations due to different muon impact distances are smaller than $\pm 0.2$ % for the average aerosol case, and smaller than $\pm 0.5$ % for the extreme high case.

5.3 Effects Related to Analysis and Reconstruction of the Muon Rings

5.3.1 Biased Signal Estimators

Many signal estimators show a small bias toward low signal amplitudes (Albert et al., 2008), which will distort the amount of reconstructed muon light. It is recommended to switch to unbiased estimators, at the latest in the last analysis step prior to determination of the optical bandwidth and impact parameter. The error introduced by keeping biased estimators (typically around 0.5–1 p.e. at signal expectations below 2 p.e.; Albert et al. (2008)) can be estimated as follows: assuming that on average one pixel on either side of the ring is affected by this bias at each azimuthal angle, the total bias is roughly 1.5 p.e. per number of pixels on the ring ( $N_{\mathrm{pix}}$ , Eq. 23), hence $1.5\cdot 2\pi\cdot\theta_{\mathrm{c}}/\omega$ , where $\omega$ is the pixel FOV. For $\theta_{\mathrm{c}}\approx 1.23^{\circ}$ and $\omega=[0.1^{\circ},0.24^{\circ}]$ , systematic biases in the range from 8% to 15% of the total signal are obtained.

5.3.2 Non-negligible Inclination Angles

Muon inclination angles were so far limited to $|\upsilon|\lesssim{\color[rgb]{0,0,0}1^{\circ}}$ for $\lesssim 5^{\circ}$ FOV cameras, with associated effects scaling as $\approx\tan\theta_{\mathrm{c}}\tan\upsilon+\tan^{2}\upsilon/2{\color[rgb]{0,0,0}\lesssim 0.05\%}$ for $|\upsilon|<{\color[rgb]{0,0,0}\lesssim 1^{\circ}}$ (see Appendix B) and hence negligible. However, wide-field cameras for CTA with up to 10∘ FOV can regularly obtain muon images with ring centers up to $3.5^{\circ}$ from the camera center, resulting in a maximum effect on the modulated signal along the ring of 0.3%.

The situation becomes more complicated when central shadowing objects are considered, like the camera. Traditionally, relatively small cameras (with small focal ratios and hence small dimensions compared to the reflector) have been employed by IACTs (Barrau et al., 1998; Cogan, 2005; Kildea et al., 2007; Anderhub et al., 2009; Aleksić et al., 2016; Giavitto et al., 2017). In these cases, the shadowing function $\Theta(\vec{\rho},\vec{\upsilon})$ can be approximated by Eq. 13. For telescopes with a focal ratio of $f=F/D$ , the shadow produced by the camera moves with inclination angle approximately as

[TABLE]

slightly more than the intrinsic resolution with which $\rho_{R}$ can be obtained (see Section 4.4). The apparent shift of the camera shadow can, however, be included in the reconstruction of the impact parameters and optical bandwidth using Eq. 43.

Furthermore, the square shape of some camera designs (like the MST and the LST) introduces dependencies of the reconstructed optical bandwidth on the inclination angle; see Figs. 16 and 17. Using the procedure described in Appendix C, the effect of such a quadratic shape can be included in the analytical model, as well as the additional shadow of the central hole (Eq. 50).

We considered the MST as the worst-case scenario for the effects of camera shadows. The MST will host a camera of $\sim 8^{\circ}$ FOV in an $f/D\sim 1.3$ configuration (Glicenstein et al., 2013; Oakes et al., 2015; Puehlhofer et al., 2015), resulting in a 3 m diameter square shadow on a 12 m diameter optical mirror structure. Figures 16 and 17 (top) show the number of photons for a muon incident on the mirror at different impact distances with $0^{\circ}$ and $2^{\circ}$ inclination angle ( $\upsilon$ ), the latter under various azimuthal projection angles ( $\psi$ ). The predictions have been produced by ray-tracing simulations of the entire telescope, including all shadows and the Cherenkov light of the muon emitted below the camera. One can clearly see that the amount of expected photoelectrons depends not only on $\upsilon$ , but also sensitively on $\psi$ , throughout a range of impact distances from $\rho_{R}\approx 0.1$ to $\rho_{R}\approx$ 0.7. If only the solution by Vacanti et al. (1994) (Eq. 26) was used, both for the chord on the mirror and the chord of the shadow (Eq. 13), peak-to-peak differences of up to 25% are observed for the simulated amount of light, only varying the $\psi$ angle. Accordingly, an error of the same magnitude can be obtained in the reconstruction (middle panels). If instead the displaced square camera shadow is used for the model, the maximum error is reduced to below 5% (bottom panels). The latter includes now all the additional shadows from ropes, the camera support structure, etc., which have not been taken into account in the analysis. One can also see in both figures that an average error of less than 2% can be expected in that case.

6 Systematic Effects due to the Different Optical Bandwidth of the Cherenkov Light Produced by Muons and Gamma-ray Showers

Given the different energy spectra of the Cherenkov light from local muons and (distant) gamma-ray showers after traversing the atmosphere, unrecognized chromatic changes in the telescope optics and camera may lead to unrecognized changes in the conversion factor $C_{\mu-\gamma}$ (defined in Eq. 17). This factor, defined as the ratio between $B_{\gamma}$ and $B_{\mu}$ , depends on the detector efficiency $\xi_{det}(\epsilon)$ and on the atmospheric transparency to Cherenkov light from muons $t_{\mu}$ and gamma-ray showers $t_{\gamma}$ , respectively.

Figure 18 shows that at wavelengths lower than $\sim$ 290 nm (photon energies greater than $\sim 4.3$ eV) the Cherenkov light from gamma-ray showers is almost completely absorbed by the atmosphere, at variance with that from muons. On the other hand, PMTs with super-bialkali photocathodes, protected by borosilicate glass, can feature considerable quantum efficiency even below 250 nm. This is the range in which muon calibration was most strongly affected by systematic errors in the past (Leroy, 2004; Fegan and Vassiliev, 2007).

Because $B_{\gamma}$ and $B_{\mu}$ contain only an integrated detector spectral efficiency $\xi_{\mathrm{det}}(\epsilon)$ , a recalibration of the telescope optical throughput would lead to wrong results for gamma-ray analyses, if the telescope loses efficiency at less than 300 nm wavelength degrading $B_{\mu}$ . Since there is practically no Cherenkov light from gamma-rays in this energy band, corrections based only on $B_{\mu}$ will overestimate $B_{\gamma}$ .

Figure 19 compares the combined detector and atmospheric efficiency for Cherenkov photons from local muons and gamma-ray showers, for two different telescope types. The loss of efficiency for gamma-ray shower light is clearly visible for the LST, which uses PMTs, whereas the effect is much lower for the SST, where the camera sensors are SiPMs. Also visible is the carefully chosen cutoff around 290 nm for the LST, produced by the camera protection window (Ambrosi et al., 2015).

The CTA Consortium has made a considerable effort to eliminate exposure of the PMT cameras to light with wavelengths below 290 nm, through the choice of an adequate protection window. That window,171717Shinkolite #000 from Mitsubishi Rayon, see https://www.m-chemical.co.jp/shinkolite/. made of polymethyl methacrylate (PMMA), filters out such short wavelengths, while maintaining a high transmittance to photon wavelengths above 300 nm (see Figure 2). We treat in the following the residual systematic effect due to chromatic degradation of individual elements, after taking into account the spectral cut produced by the camera protection window.

Since optical elements tend to degrade in a chromatic manner, losing efficiency rather toward the blue and UV part of the spectrum, we consider here the extreme case of a complete sensitivity loss (“blindness”) below a given wavelength while maintaining full efficiency throughout the rest of the spectrum. Such a complete unrecognized loss of sensitivity is very improbable, compared to partial losses. In that sense, this procedure provides an upper limit on the possible error that can be made by neglecting the effect of chromatic degradation of optical elements, neglecting additional possible, though unexpected and less harmful, chromatic effects at longer wavelengths. Note that an additional degradation at longer wavelengths normally reduces that error because the sensitivity loss will become less chromatic. Also, the difference between muon- and gamma-ray-generated Cherenkov light spectra becomes negligible at long wavelengths.

In that hypothetical extreme scenario, the measured optical bandwidth degrades to

[TABLE]

An analyzer would (wrongly) conclude that the telescope’s optical throughput will need to be corrected by the factor

[TABLE]

where $\epsilon_{\mathrm{min,orig}}$ and $\epsilon_{\mathrm{max,orig}}$ are the originally assumed detector sensitivity limits (or atmospheric cutoff, if applicable). The error made by the analyzer is then

[TABLE]

The results of Eq. 73 are shown in Figure 20, for cameras both with and without the protection window.

One first notes that even in the high-energy limit, where supposedly no degradation has occurred, the error does not converge exactly to zero. This has to do with the chosen integration limits $\epsilon_{\mathrm{min,orig}}$ and $\epsilon_{\mathrm{max,orig}}$ for the reference atmosphere. One can calibrate the offset with MC simulations.

An unrecognized chromatic degradation of the detector elements, expressed by the rise of $\Delta B_{\gamma}/B_{\gamma}$ when moving to smaller energies (longer wavelengths), produces an intrinsic error of about 10% if the detector becomes blind below wavelengths of 400 nm for the PMT-based telescopes, and below 3% for the SiPM-based ones. If the camera protection window is left out, errors double.

Several possible causes can contribute to the chromatic changes in optical elements of the telescopes. Their effects are summarized in the following:

•

Chromatic degradation of the light detectors: photomultiplier tubes (PMTs) may have considerable quantum efficiency below 290 nm, if UV-transparent glass (fused silica or borosilicate) is used for the photocathode (Kildea et al., 2007; Otte et al., 2011). This is the range where Cherenkov light from air showers is almost completely absorbed, but the one from muons still gets through to the camera. Experience with the Whipple Hamamatsu R 1398 photomultipliers has shown that the photon detection efficiency degraded from 5% to 20% (depending on the individual PMT tube), in the wavelength range from 290 to 450 nm (Daniel, 2015). VERITAS observed a 10% sensitivity increase below 300 nm and a 10%–30% drop above 550 nm with the Hamamatsu R10560-100-20 and claims a serious aging of spectral sensitivity of the previous Photonis XP2970 (Gazda et al., 2016). If such a behavior occurs with the CTA PMTs, the muon calibration will overcorrect the loss by 1%–2% for the PMT-equipped cameras.

•

Chromatic degradation of the PMMA: the Shinkolite #000 PMMA, currently favored for use in both LST and MST cameras, can permanently lose up to 45% transmittance in the range from 290 and 450 nm, after strong exposure to UV light, whereas the rest of the light spectrum is unaffected (after some recovery time; Förster (2015); Motohashi (see also Figs. 4.5–4.10 of 2016)). In this case, the muon calibration will overcorrect the loss by $\sim 3-4$ % for an affected MST or LST, both for the worst-case scenario of an extreme exposure to sunlight. Direct exposure of the material to the Sun during 10 days on the CTA-N observatory has not, however, been able to reproduce such a chromatic transmittance loss, but instead a constant drop by 2–3% (Ambrosi et al., 2015).

•

Chromatic degradation of the focused mirror reflectance: experience from the H.E.S.S. mirrors has shown that reflectance losses affect the wavelength range around 300 nm 5–10% more strongly than at wavelengths around 500 nm (Förster, 2015). Such a scenario would translate into an error in the muon calibration of the order of 1% for MSTs and LSTs and $<$ 1% for the SSTs (unless new SiPMs are developed with much higher PDE around 300 nm).

Chromatic effects can be characterized from time to time by external devices, as foreseen for CTA (Segreto et al., 2016). If the telescope’s optical throughput is characterized with respect to its relative response to two reference (laser) wavelengths, above and below 400 nm (e.g. the doubled and tripled emission line of an NdYAG laser at 532 and 355 nm, respectively), with a precision of the order of 10% for the response ratio between the two wavelengths, the unrecognized chromaticity effects will add maximally 1% to the systematic uncertainty of the muon calibration method.

7 Secondary Effects Broadening the Ring

Broadening of the ring influences the detectability of muons (by both the telescope trigger and the offline analysis) and the reconstructed optical PSF of the telescope.

Normally, muon rings become detectable if their radii become considerably larger than the ring width, i.e. $\Delta\theta_{\mathrm{c}}/\theta_{\mathrm{c}}\ll 1$ . This is the case for energetic muons with large Cherenkov angles. Most effects that cause a broadening of the ring width have been described by Vacanti et al. (1994), but they are often overestimated. We review them here and estimate their magnitudes for the CTA telescopes. Additionally, we add further, hitherto overlooked, effects.

7.1 Instrumental Effects

7.1.1 Finite Camera Focuses

Finite camera focuses were not treated by Vacanti et al. (1994). We calculate the variance of the ring displacement produced by a camera focused at $u_{f}$ (Eq. F9) and obtain

[TABLE]

For telescopes focused at 10 km, this effect broadens the ring of an LST by $<5$ %, that of an MST by $<3$ %, and $<1$ % for the case of an SST.

7.1.2 Optical Aberrations

The contribution of the spherical mirror aberration to the broadening of the ring has been described as constant by Vacanti et al. (1994). Additionally, coma and astigmatism must be taken into account. According to Fegan and Vassiliev (2018) (see also Schliesser and Mirzoyan (2005) for the pure parabolic case and Vassiliev et al. (2007) for the ideal Davies-Cotton (DC) telescope), aberrations for a single-mirror telescope can be approximated from third-order aberration theory as

[TABLE]

where $\Delta\xi$ denotes the tangential and $\Delta\eta$ the sagittal standard deviation of the PSF. The angle $\upphi$ denotes inclination with respect to on-axis. We have introduced here also a common spread of $\Delta_{0}$ subsuming the intrinsic optical quality of the individual mirror facets (mainly due to form deviations from ideal mirror shapes and to surface roughness; Tayabaly et al. (2015)), the aberrations caused by the fact that most facets are inclined relative to the optical axis and the precision with which the individual mirror facets are aligned (Cornils et al., 2003). Spherical aberrations of each (spherical) mirror facet scale with $r^{3}/4R^{2}$ , where $r$ is the radius of the individual mirrors and $R$ its radius of curvature (Lewis, 1990; Cornils et al., 2003; Vassiliev et al., 2007); $N_{\mathrm{mirrors}}$ denotes the number of mirror facets (see Table 1).

Optical aberration will become the dominant effect for the ring broadening of images from the highest-energy muons for all CTA telescopes, although requirements on their optical PSF are rather strict (see Table 5).

Moreover, the PSF of a typical IACT varies with pointing zenith angle $\theta$ , because a support structure that does not deform under the influence of gravity would be too costly. Cornils et al. (2003) found that the size of the PSF of the H.E.S.S.-I telescopes can be parameterized as

[TABLE]

where $\Delta(\xi,\eta)_{0}$ is the PSF value with the telescope pointing to zenith and $d\approx 1$ mrad (in the case of H.E.S.S.-I) is a parameter describing the influence of the deformation of the mirror support under inclination. The LST will, however, make use of an active mirror control system (Garczarczyk, 2006; Biland et al., 2008; Aleksić et al., 2016) that effectively removes the zenith angle dependency $d$ .

Additionally, wind may degrade the quality of the PSF, particularly the larger telescopes. Careful simulation studies have shown, however, that the magnitude of related effects is of the same order as or smaller than $\Delta_{0}$ (Ambrosi et al., 2015; Farnell et al., 2015).

Figure 21 (top) shows an example for the expectation of tangential and sagittal rms point spread for each telescope type, together with the requirements. A muon ring of radius $\theta_{\mathrm{c}}$ is affected by a wide range of these parameters. The average radial distribution of the ring depends then on the inclination angle, the impact distance, and particularly the azimuthal component of the inclination parameter with respect to the azimuthal component of the impact parameter.

We will quantify this effect now by weighting each ring segment with the chord $D$ and the projected component of the tangential $(\sin\delta)$ and sagittal $(\cos\delta)$ aberration onto the ring, evaluated at the distance $d$ of each point on the ring to the camera center (see Appendix D):

[TABLE]

Figure 21 (bottom) shows the resulting average ring width, for an impact distance of $\rho_{R}=0.9$ . As expected, the ring width depends not only on the muon inclination angle but also on the relative azimuthal components of the impact parameter and inclination angle. The effect of the latter can even create a peak-to-peak spread larger than the average ring width at a given inclination angle, particularly for the LST.

These results show that in order to improve precision, both inclination and impact parameter must be taken into account in a ring width analysis, leading to an estimate of the PSF of a CTA telescope.

7.1.3 Discrete Pixel Width

The discrete pixel is mentioned here for completeness but introduces a negligible uncertainty on the reconstruction of the muon ring in our case (Vacanti et al., 1994):

[TABLE]

where the number of hit pixels, Eq. 23, has been used. Given the fine pixelation of the CTA cameras, the relative ring broadening introduced by this effect is always below 5%.

7.2 Atmospheric Effects

7.2.1 Multiple Scattering

Multiple Coulomb scattering of particles in the air along their trajectory causes the emitted Cherenkov photons to appear scattered around their mean positions in the focal plane, according to an approximately Gaussian distribution. The standard deviation of that distribution can be described as (Lynch and Dahl (1991); Patrignani (sect. 34.3 of 2017))

[TABLE]

where $p$ is the momentum of the muon and $X_{0}$ its radiation length ( $\approx$ 367 kg/m2). The density of air is approximated by $\rho_{0}\cdot\exp(-h/H_{0})$ , where $\rho_{0}\approx$ 1.225 kg m*-3* and $R_{0}=X_{0}/(\rho_{0}\cdot\exp(-H_{\mathrm{obs}}/H_{0}))\approx 378$ m for an observatory altitude $H_{\mathrm{obs}}$ of 2200 m a.s.l. and $H_{0}\approx 9.7$ km. Since the ring probes different path lengths in the atmosphere, according to its azimuth value, we calculate an expectation value for the average ring width (see Appendix H):

[TABLE]

The new Equation 98 yields only about half the contribution predicted by Vacanti et al. (1994), who assumed that the full muon path is seen by the telescope through the entire ring, and a constant atmospheric density. Figure 22 shows an example for its expectation for each telescope type. One can see that the peak-to-peak variation, due to different impact distances, amounts to a bit less than half the mean ring width produced by this effect.

7.2.2 Variation of the Refractive Index

The refractive index of air varies with altitude, temperature, atmospheric pressure, relative humidity, and wavelength. Dependencies of the refractive index of air on these parameters have been discussed in great detail in Ciddor (1996, 2002) and later in Tomasi et al. (2005). Atmospheric conditions, especially temperature and water vapor content, vary with time and altitude, but also, to a lesser extent, with the concentration of CO2. In general, a good description is obtained by assuming a variation of $\varepsilon=(n-1)$ with the density of air, which in turn can be approximated by an exponential decrease with a scale height of $H_{0}=9.7$ km 181818Note here that (Vacanti et al., 1994; Leroy, 2004) assumed erroneously a scale height of 8.4 km which is valid only for the atmospheric pressure. at La Palma (Gaug et al., 2017). Assuming an exponential scaling law $\epsilon$ , and $\Delta\theta_{c}\simeq\Delta x\cdot\partial\theta_{\mathrm{c}}/\partial h$ , we obtain for the standard deviation of the angular deviation (see Appendix I):

[TABLE]

Eq. 99 differs from the one proposed by Vacanti et al. (1994) by the square root factor, which ranges from $\sim$ 0.3 for $\rho_{R}=0$ to $\sim$ 0.5 for $\rho_{R}=1$ .

Additionally, since $\theta_{\infty}^{2}\approx 2\varepsilon$ , dependencies of the refractive index and hence the Cherenkov angle on the atmospheric conditions can be expected. Drifts of $\theta_{\infty}$ calculated using formulae in Tomasi et al. (2005) are shown in Figure 23. The temperature dependency is the largest effect, affecting $\theta_{\infty}$ by up to 7%, whilst only negligible dependencies of the refractive index on humidity and CO2 concentration are observed.

Since the Cherenkov angle is reconstructed directly from ring radius and the accumulated signal in the camera, $Q_{\mathrm{tot}}$ , scales linearly with the Cherenkov angle, the residual bias on the reconstructed optical bandwidth is expected to be negligible. Nevertheless, Bolz (2004) found measurable effects with H.E.S.S. data of the same size as those suggested in Figure 23, implying that temperature- and pressure-dependent effects should be monitored, at least during the commissioning phase.

7.2.3 Energy Loss of the Muon

Muons lose energy through ionization, which reduces the Cherenkov angles with path length. This effect has been neglected in Eqs. 7 and 26, but introduces a small bias into the results.

It is convenient to write the average rate of muon energy loss as (Barrett et al., 1952)

[TABLE]

Here $a(E)$ is the ionization energy loss and $b(E)$ the sum of $e\!^{+}e\!^{-}$ -pair production, bremsstrahlung, and photonuclear interactions. The second term $b(E)\cdot E$ dominates energy loss only above 1.1 TeV in air (Groom et al., 2001), where the muon flux is already very small, compared to energies in the range from 10 to 100 GeV. The ionization energy is given by the Bethe-Bloch formula (see, e.g., sect. 34.6 of Patrignani, 2017) and amounts to about $-\mathrm{d}E/\mathrm{d}x=2.15$ MeV g*-1cm2* for a muon producing a Cherenkov angle of $\theta_{\mathrm{c}}=0.4^{\circ}$ and $-\mathrm{d}E/\mathrm{d}x=2.43$ MeV g*-1cm2* for $\theta_{\mathrm{c}}=1.2^{\circ}$ .

Assuming $\Delta\theta_{c}\simeq\Delta E_{\mu}\cdot\partial\theta_{\mathrm{c}}/\partial E_{\mu}$ , one obtains the variation of $E_{\mu}$ while the particle traverses the path length $L$ losing $\mathrm{d}E/\mathrm{d}x$ . Calculating expectation values for the changes on the ring width (see Appendix J), one obtains for the average ring width

[TABLE]

where the muon rest mass $m_{\mu}\simeq 105.7$ MeV has been used and the relationship $\theta_{\mathrm{c}}^{2}\simeq\theta^{2}_{\infty}-(m_{\mu}/E_{\mu})^{2}$ . As before, the average ring width results are smaller than that estimated by Vacanti et al. (1994) by a factor of 0.3 to 0.5.

Figure 22 shows examples for its expectation for each telescope type. One can see that this contribution to the ring broadening is smallest at high muon energies.

7.2.4 Bending of the Muon Trajectory in the Geomagnetic Field

We calculate the variance of the angular shift (Eq. 61), introduced by the bending of the muon path in the Earth’s magnetic field (see Appendix G):

[TABLE]

For a muon of $p=10$ GeV in a magnetic field perpendicular to its velocity, the average fractional ring broadening amounts to up to 0.3 for an LST, for the worst case of a telescope pointing perpendicular to the Earth’s magnetic field and detecting a small ring (see Figure 22).

7.2.5 Combined Effects

Natural ring broadening is dominated by multiple scattering of the muon, with broader rings at small impact distances and large observation zenith angles. On the other hand, variation of the refractive index, energy loss, and the magnetic field create broader rings for larger impact distances. The net effect yields slightly broader rings for smaller impact distances and large observation zenith angles.

Figure 24 shows the combined effects of the ring broadening for all telescope types, together with the range of values expected from aberrations. As expected, natural broadening dominates small rings, whereas telescope aberrations dominate the largest rings. There is, however, a broad transition region where either one or the other may dominate, according to the muon impact distance and inclination angle.

8 Muon Spectra and Expected Event Rates

The muon spectrum at ground level in the range 10 GeV – 1 TeV was published by Hebbeker and Timmermans (2002) and modeled with a parameterization based on the theoretical formula proposed by Bugaev et al. (1998):

[TABLE]

where $\Phi$ is the muon flux, $p_{\mu}$ is the muon momentum, and $H_{1}=0.133\pm 0.002$ , $H_{2}=-2.521\pm 0.002$ , $H_{3}=-5.78\pm 0.004$ , $S_{2}=-2.11\pm 0.03$ , and $C=0.86\pm 0.06$ . The parameters $H_{1}$ , $H_{2}$ and $H_{3}$ denote the logarithm of the differential flux at 10, 100, and 1000 GeV, $S_{2}$ represents the exponent of the differential flux at 100 GeV, while $C$ is the absolute normalization of the muon spectrum.

More recent data are given by Schmelling et al. (2013) for the range 100 GeV – 2.5 TeV and by the IceCube Collaboration (Aartsen et al., 2016) at higher energies ( $E_{\mu}>15$ TeV ). The latter is reproduced by the following formula:

[TABLE]

The muon flux dependence with altitude has been measured by balloon experiments (e.g. De Pascale et al., 1993; Bellotti et al., 1996; Haino et al., 2004). The flux for muon momenta above 10 GeV can be parameterized by the following formula at altitudes less than 1000 m (Hebbeker and Timmermans, 2002):

[TABLE]

where $h$ is the altitude in meters.

The dependency of the muon flux on zenith angles computed using air shower simulations can be described as

[TABLE]

where the parameter $a(p_{\mu})$ is $-1.5$ at $p_{\mu}=10$ GeV and $-1.28$ at $p_{\mu}=30$ GeV (Hebbeker and Timmermans, 2002). The behavior for energies $E_{\mu}\cdot\cos\vartheta>100$ GeV and zenith angles $<70^{\circ}$ is better described by the following formula (Sect. 29 of Patrignani, 2017):

[TABLE]

We estimate, in the following, per-telescope mono trigger muon event rates, required by CTA to get triggered in an unbiased way and flagged, unless another neighboring telescope triggers at the same time and activates a stereo trigger. Such a requirement ensures that the stereo trigger bias is eliminated and muon rates are kept at a useful level for calibration. Only the LSTs will be allowed to detect muon images only under a stereo trigger.

To compute an estimate for the maximally achievable muon rates, we assumed a reference observatory altitude of 2200 m a.s.l., a minimum muon energy of $E_{\mathrm{min}}=8$ GeV (i.e. $\theta_{c}\gtrsim 1.0^{\circ}$ ), an impact distance value contained within $0.5<\rho_{R}<1$ , and two reference observation zenith angles: $\vartheta=0^{\circ}$ and $\vartheta=60^{\circ}$ . The energy spectra relative to the two considered angles are plotted in Fig. 25 as a function of the muon energy. Using Eqs. 103 and 106 for energies lower than 100 GeV and Eq. 107 above this energy, we calculate a maximally achievable rate $R_{\mathrm{max}}$ of useful muon images:

[TABLE]

where $\rho_{\mathrm{max}}$ is the maximum impact parameter and $\upsilon_{\mathrm{max}}$ the maximum muon angle to achieve a fully contained ring in the camera. The latter has been defined as $\textit{FOV}/2-\theta_{\infty}-0.2^{\circ}$ to account for the broadening of the ring.

Under these assumptions, the estimated maximally available rate of useful muons ranges from $\sim$ 0.4 Hz (for the SSTs) to $\lesssim$ 10 Hz (for the LSTs) (see Table 6). These values are sufficiently high to yield acceptable statistics for nightly calibration.

Requiring a statistical precision of $\leq$ 1% for the derived optical bandwidth with muons $B_{\mu}$ , and using 5% for the events-wise resolution of $B_{\mu}$ for LST and MST (see, e.g., Table 2) and 20% for the SSTs, we estimate a required minimum number of $\sim$ 25 and $\sim$ 400 usable muon images for the larger telescopes and for the SST, respectively. With these estimates, we derive minimum time scale $T_{\mathrm{monitor}}$ within which optical bandwidth estimates can be updated for monitoring purposes. These time scales lie always below the typical run time of an observation, albeit one order of magnitude larger for the SSTs, as expected. All results are shown in Table 6.

We calculate now the rates of muon images that are required to contain a given pixel in the camera, such that we can use it for an alternative camera flat-fielding approach (Bolz, 2004). To do so, we use Eqs. 23 and 26, and obtain

[TABLE]

where $\omega$ and $N_{\mathrm{pixels}}$ can be read off from Table 1. In order to estimate the monitoring time scales, we realize that the resolution is now dominated by the intrinsic Poissonian fluctuations of $Q_{\mathrm{pix}}$ , which will be of the order of $\pm\sqrt{60}~{}\mathrm{p.e.}$ for the LST and MST and $\pm\sqrt{10}~{}\mathrm{p.e.}$ for the SST, assuming a ring width smaller than one pixel. This entails again $E_{\mathrm{min}}\approx 8$ GeV and an event-wise resolution of $\sim$ 0.13 (LST and MST) and $\sim$ 0.32 for the SST. Here we ask for 3% resolution of $Q_{\mathrm{pix}}$ (a typical value for the accuracy of the pixel-wise gain adjustments (Aharonian et al., 2004; Aleksić et al., 2016)) and end up again with at least 25 events for the larger telescopes and 150 events for the SST. As shown in Table 6, considerably longer monitoring time scales are necessary for a pixel-wise calibration update, albeit well within one typical observation night.

In order to assess rates of muon images that can be used to determine the average time delay of a given pixel in the camera, with respect to the rest of the camera, we use Eq. 109, albeit with a slightly relaxed minimum ring radius of $\theta_{\mathrm{c}}>0.7^{\circ}$ , hence $E_{\mathrm{min}}\approx 6$ GeV, and no lower impact distance ( $\rho_{\mathrm{min}}=0$ ). To determine the number of events needed to be accumulated, we realize that one event is already sufficient to obtain a time resolution of $<$ 1 ns for an LST or MST pixel, given the isochronous parabolic reflector (or the modified Davies-Cotton design of the MST, optimized to minimize timing dispersion; Farnell et al. (2015); Fegan and Vassiliev (2018)), and the transit time spread of the PMTs used (Toyama et al., 2015). For the SST, we assume a time spread $<$ 1.5 ns across the camera (della Volpe et al., 2016; Fegan and Vassiliev, 2018), and the excellent timing characteristics of the SiPM (Piemonte et al., 2012) and arrive at the same conclusions. This is sufficient to detect a single pixel’s time behavior change with respect to its camera average (however, not for gradual drifts across the camera, for which more statistics are required).

Finally, Table 6 also lists maximum rates with which the PSF of a telescope can be monitored. For this, we use the results of the previous section and require muon rings with $\theta_{\mathrm{c}}>1.15^{\circ}$ ( $E_{\mu}>14.3$ GeV for the LST), $\theta_{\mathrm{c}}>1.1^{\circ}$ ( $E_{\mu}>11.2$ GeV for the MST), and $\theta_{\mathrm{c}}>0.9^{\circ}$ ( $E_{\mu}>7.2$ GeV for the SST). On the other hand, impact distance cuts do not need to be so strict for this purpose, allowing $\rho_{R}>0.2$ . Here we assume 30% event-wise resolution for all telescopes (see, e.g., Fig. 9)191919We multiply the standard deviation by $\sqrt{\pi/2}$ to account for the variance of the 0.2 quantile instead of the mean, as outlined in Section 4.9. , but only 5% desired overall resolution, given the systematics outlined in the previous section. Hence, about 40 events needed to be accumulated within $T_{\mathrm{monitor}}$ for all telescopes.

The estimates of Table 6 have to be understood as upper limits since neither trigger losses nor analysis cut efficiencies have been included. However, realistic estimates for these effects cannot be too far from unity, if an unbiased trigger and analysis scheme needs to be achieved. For this reason, we believe that the rates and monitor times of Table 6 predict at least a correct order of magnitude.

9 Discussion and available precision

We summarize the expected systematic uncertainties in Table 7 and obtain an overall accuracy of the optical throughput calibration using muon rings, after adding each contribution quadratically. The obtained numbers require careful analysis and corrections, going well beyond the algorithms presented in the past. Nevertheless, the resulting accuracy is considerably worse than that estimated in the past.

The situation is, however, not as dramatic if we assume that telescopes arrive very well calibrated and undegraded on site, and a careful first correction can be made by comparing MC simulated muon images with those obtained during telescope commissioning. The first eight items in Table 7 can then be assessed and corrected. For the subsequent monitoring of the degradation of the telescope, only the latter parts are important (provided that the hardware is not modified or replaced).

The last row of Table 7 shows these uncertainties applicable to the monitoring of the optical bandwidth. One can see that the numbers are already quite acceptable for the LST and the MST.

In order to achieve this accuracy, several hardware requirements have to be met by the telescope designers:

•

An unbiased event flagging and forced readout has been required for muon images triggering only one single MST or SST telescope. This measure has been taken to ensure that sufficient high-quality ring images are obtained in order for the muon calibration to provide at least a monitoring value per observation night and to avoid trigger biases (particularly the stereo trigger bias).

•

PMT-based cameras have selected a camera protection window that is sufficiently opaque to light below 290 nm wavelength. This choice ensures that the overcorrection produced by muon calibration due to chromatic degradation of the optical elements of a telescope remains below the values provided in item 19 of Table 7.

•

Additionally, CTA will be equipped with a dedicated device to monitor the wavelength dependency of the telescopes’ optical throughput (Segreto et al., 2016).

•

Generally, requirements on the quality of the optical system, required to avoid significant biases in the reconstructed ring parameters, are much stricter than that achieved with the current generation of IACTs.

Moreover, the muon analysis must apply several additional calibration steps, which were not among the standard procedures in the past:

•

A plate-scale calibration must be carried out to a precision better than 1%, in order to eliminate biases due to coma and astigmatism of the primary mirrors.

•

The effect of the finite camera focus must be corrected, according to Eq. 55.

•

A strong impact distance cut of $\rho_{R}>0.6$ must be applied to reduce the effect of the shadows from the large-FOV cameras.

•

The atmospheric conditions, particularly aerosol extinction, must be regularly monitored, and data with large boundary layer extinction excluded. Atmospheric transmission must be included (at least on average) in the derivation of the optical bandwidth and its transformation to the bandwidth seen by Cherenkov light emitted from gamma-ray showers.

•

The readout hardware and pulse integration software must ensure that the full muon light pulse gets integrated in an unbiased way. This requires modifications of the typical low-level analysis carried out for gamma-ray analysis.

A selection of the analyzed muon images can also be used for monitoring purposes: the overall quality of the telescope PSF can be monitored for muon rings with $\theta_{\mathrm{c}}>1.15^{\circ}$ (for the LST), $\theta_{\mathrm{c}}>1.1^{\circ}$ (MST) and $\theta_{\mathrm{c}}>0.9^{\circ}$ (SST). We predict event rates that allow such monitoring to be carried out with sufficient precision of $<$ 1% on time scales of minutes to the LST and the MST and less than 5 minutes for each SST. Preferably, events are selected that were recorded at low observation zenith angle and have large impact distances and small inclination angles.

With very large statistics of selected muon events, a pixel-wise monitoring of relative photon detection efficiencies should be envisaged to obtain an independent and complementary cross-check for the standard flat-fielding procedures. We predict that such a data set will be available after less than 10 minutes (LST and MST), and within one average observation night (SST) with a resolution of 3%.

10 Conclusions

Muon calibration has been used as both an absolute calibrator and a relative monitor of the optical throughput of all IACTs in the past. Using the Vacanti et al. (1994)’s algorithms and a series of suitable estimators for the ring parameters and the impact distance of the muon, several authors have claimed a theoretical accuracy of as good as $\approx 2$ % in the past. We have carefully reviewed all documented algorithms and systematic effects and found several that have not been taken into account so far, or which appear non-negligible for the large mirrors or large square-shaped cameras of CTA, or the dual-mirror telescope designs. Most of these effects can be corrected by modifications of the standard analysis at several critical steps, like the pulse integration, pixel selection, carefully selecting the remaining ring images, and using new models for the effect of camera shadows and central holes in the reflector surface. Nevertheless, a few of the systematic effects will require dedicated studies or even adaptation in the design of the hardware employed for the CTA telescopes to achieve the desired accuracy. These effects are mostly related to the differences in the observed spectrum of Cherenkov light received from local muons and distant gamma-ray showers. The conversion of measured optical bandwidth with the muon method to its equivalent for the telescopes’ response to Cherenkov light from gamma-ray showers, to account for wavelength-dependent degradation, requires hardware that efficiently cuts the camera response to photon wavelengths below about 290 nm. Finally, close comparison with MC simulated telescope response to Cherenkov light emitted by local muons is necessary at the beginning of operation, to correct for residual shadows and spaces between individual mirrors.

We find that after such a carefully prepared calibration scheme, the optical bandwidth of the full telescope can be determined with a resolution of $<$ 1% (LST and MST) and $<$ 5% (SST) on time scales of better than a minute and an accuracy of better than 4% (LST and MST) and 5% (SST). This is considerably better than other proposed methods so far (Gaug, 2014; Brown, 2018; Stefanik et al., 2019) and only slightly worse than the cross-calibration scheme proposed by Mitchell et al. (2016b), which allows, however, only for a relative calibration between telescopes or telescope types. Other direct methods (Segreto et al., 2016) perform even better and show more flexibility but require a dedicated setup and cannot be carried out at the same time as the telescopes perform science observations. The muon calibration method is therefore the first option for online and offline monitoring of the optical bandwidth of all CTA telescopes (Gaug et al., 2014).

Muon calibration allows one to monitor additionally, and without much effort the flat-fielding of the camera with a resolution of 3% on time scales of at least one typical 20-minute data run (for the LST and MST) and within one typical observation night for the SSTs. For those telescopes that record the event arrival times, the muon images will provide the time resolution of each pixel and can serve to determine the relative time offsets with a precision of better than 1 ns within time scales of less than 2.5 minutes. Similarly, differences in reflectance among the mirrors can be detected by plotting the retrieved optical bandwidth as a function of reconstructed impact point.

Following the approach of the MAGIC analysis, the 0.2 quantile of the ring width distribution shall be used to monitor changes in the telescopes’ optical PSF. Contrary to stars, the light from muon rings gets registered within sub-nanosecond time windows, which allows one to greatly reduce residual backgrounds from the night sky and unresolved stars. The stronger concentration of Cherenkov light from local muons close to their impact points may even enable point-spread monitoring resolved to individual mirror facets or groups of mirrors. We have shown that the larger telescope mirrors, together with very strict requirements on the optical quality of the reflector, make such a monitoring scheme more challenging for CTA than in previous IACTs. Particularly, the LSTs need to select only the largest rings with a radius $\theta_{\mathrm{c}}>1.0^{\circ}$ to be able to disentangle degradations of the optical point spread from the natural ring-broadening effects. Further selection of impact distances and the orientation of the telescope with respect to the Earth’s magnetic field may be necessary. Even under these conditions, monitoring points can be obtained at least per minute with a precision of 5%.

The authors thank the internal CTA referees for attentive reading and comments on the paper. This work has been funded by the grant FPA2015-69210-C6-6-R of the Spanish MINECO/FEDER, EU. The CTA Consortium gratefully acknowledges financial support from the agencies and organisations listed at https://www.cta-observatory.org/consortium_acknowledgments. This paper has gone through internal review by the CTA Consortium.

Appendix A Useful Relations

Several useful relations for calculations of muon Cherenkov light have been used in the following sections of this appendix. All of them assume that in air $(n-1)=\epsilon\ll 1$ ,

[TABLE]

Appendix B Derivation of the Expected Photon Yield for Inclined Muons

Inclined muons are observed over a path length $L$ different from $D(\vec{\rho},\vec{\upsilon})/\tan\theta_{\mathrm{c}}$ (see Fig. 26), namely,

[TABLE]

where $\upsilon_{1}$ is the projection of the angle $\vec{\upsilon}$ onto a vertical plane rotated by $\psi$ with respect to the $\phi$ .

Since the cosine of the azimuthal projection of the incidence angle minus the azimuth part of the chord varies in the range $[-1,+1]$ , the received signal gets modified by maximally

[TABLE]

where a Cherenkov angle of $1.23^{\circ}$ has been inserted in the last line.

Appendix C Derivation of the Expected Photon Shadow for a Quadratic Camera

C.1 Impact Point within the Square

The chord $D(\rho,\phi,\phi_{0})$ of a square-shaped camera with half-diameter $A$ can be expressed as the one of a roundish camera of equivalent radius $A$ (see e.q. Eq. 13), plus the contribution $b$ between circle and square (see Fig. 27). Whenever the length $L$ achieves $A$ , the triangles are swapped and $b$ is calculated until the upper/lower edge of the square or the right/left edge, respectively. The swap condition $S$ can be obtained by requiring that the height $H$ is found in the range $[-1,1]$ :

[TABLE]

from which the new full chord $D(\rho,\phi,\phi_{0})$ can be obtained,

[TABLE]

C.2 Impact Point outside the Square

For an impact point outside the square ( $\rho_{A}>1$ ), we determine first the four angles $\phi_{\textit{XY}}$ under which the corners of the square are seen from the impact point:

[TABLE]

where $\Theta(x)$ is the Heaviside function. In order to ensure that all four angles are correctly sorted, we require, additionally, that if any of the four angles are negative and any of the four angles are larger than $\pi$ , then

[TABLE]

The visible angular range is then found between $[\min(\phi_{\textit{XY}}),\max(\phi_{\textit{XY}})]$ . The sorted four angles enclose three angular ranges: within the central range, the square is fully crossed by the chord, whereas the other two ranges lead to a chord that only cuts through an edge of the square. We need to treat the three cases separately. As above, a swap condition will be established to distinguish whether the chord crosses the square from the bottom to the top, or from the left to the right:

[TABLE]

Then, the length of the chord fully crossing the square is

[TABLE]

and the length of the chord partly crossing the square is (i.e. between the first and second sorted corner angles and between the third and fourth)

[TABLE]

The condition if the chord exits crossing an upper/lower or a left/right side of the square was determined with the following condition for $\phi_{0}$ and the corner angles $\min(\phi_{XY})$ (for the chord crossing between first and second corner angles) and $\max(\phi_{XY})$ (for the chord crossing between third and fourth corner angles), respectively,

[TABLE]

This condition determines whether the impact point is found beside the square or above or below. If not (i.e. in one of the four quarters in diagonal direction from the square), the orientation of the first and last corner point angles determines the crossing mode.

Figure 29 compares the square shadow with a roundish one for different azimuthal angles.

Appendix D Derivation of the Tangential and Sagittal Weight along the Muon Ring

We consider that a muon analysis projects the ring’s radial distribution along the azimuthal angle $\phi$ . The distance $d$ of a point on the ring to the camera center and the angle $\delta$ of the radial projection with the line connecting the ring with the camera center can be expressed as

[TABLE]

The weight given to the tangential component, evaluated at the corresponding distance of the ring from the camera center, is then $\sin(\delta)$ , whereas the sagittal component gets weighted by $\cos(\delta)$ .

Appendix E Model of a Nonuniform Pixel Response to Light under Different Incidence Angles

The following model has been adopted from Fegan and Vassiliev (2007) and slightly modified. We can model the amount of light in a muon image under the reasonable assumption that the response of the instrument, $\Omega(r)$ , approximates to a quadratic function of the distance from the center of the reflector:

[TABLE]

where the distance from the mirror center, in units of mirror radius $r_{R}=r/R$ and $\beta$ describes the enhancement $(\beta>0)$ or suppression $(\beta<0)$ of the response for photons reflected from the edge of the reflector. The function is normalized such that the response to a distant source that illuminates the full reflector uniformly is independent of $\beta$ .

[TABLE]

The response function is shown in Figure 11 (left). With $\beta<0$ the response to photons reflected from the center of the telescope is enhanced with respect to those from the edge. If $\beta>0$ , the efficiency from the edge of the reflector is enhanced.

The length of the chord $D$ (Eq. 10), convolved with the response of the instrument, is then

[TABLE]

where

[TABLE]

Solving the integral Eq. E3 yields

[TABLE]

Then, integration over $\phi$ yields

[TABLE]

and $U_{0}$ and $\Phi$ as defined in Eq. 27. For $\rho_{R}<1$ , the integral $E_{2}$ can be translated into a combination of Legendre elliptic integrals of first $F(\rho_{R})$ and second $E(\rho_{R})$ kind, namely,

[TABLE]

Appendix F Derivation of Bias and Variance of Ring Width from Finite Camera Focuses

From Eq. 55, we derive the net change in observed Cherenkov angle, if the telescope is focused to a distance $x_{f}$ , instead of infinity. Then,

[TABLE]

where $x$ denotes the distance of the photon emission point to the telescope.

We calculate the expectation values for the net shift:

[TABLE]

To derive the variance, we calculate the expectation value of $(\Delta\theta_{\mathrm{c}}/\theta_{\mathrm{c}})^{2}$ :

[TABLE]

Appendix G Derivation of Bias and Variance of Ring Width from Muon Bending in the Earth’s Magnetic Field

The angle under which the muon Cherenkov light is imaged in the camera is shifted by (see, e.g., Eq. B2, using $\tan\theta_{\mathrm{c}}\simeq D/L$ )

[TABLE]

where a small bending effect (i.e. $\mathrm{d}\theta/\mathrm{d}x\cdot L\ll\theta_{\mathrm{c}}$ ) has been assumed in the last line. Using realistic values for the Earth’s magnetic field at the CTA sites, the approximation is then better than 5% for all telescopes and useful muon Cherenkov angles.. Moreover, assuming constant bending (i.e. $\mathrm{d}\theta/\mathrm{d}x=\mathrm{const}$ ), the net shift, averaged over the ring, has an expectation value of

[TABLE]

To derive the variance, we calculate the expectation value of $(\Delta\theta_{\mathrm{c}}/\theta_{\mathrm{c}})^{2}$ :

[TABLE]

Appendix H Derivation of the Variance of Ring Width from Multiple Scattering of the Muons

We denote with $h$ the altitude of the emitted Cherenkov light above the telescope and with $\vartheta$ the telescope’s pointing zenith angle. The distance of the photon emission point to the telescope is $x$ . Multiple scattering of the muon produces an azimuth-dependent average widening of the ring according to

[TABLE]

where the density of air is approximated by $\rho_{0}\cdot\exp(-h/H_{0})$ , and $\rho_{0}\approx$ 1.225 kg m*-3* and $R_{0}=X_{0}/(\rho_{0}\cdot\exp(-H_{\mathrm{obs}}/H_{0}))\approx 378$ m for an observatory altitude $H_{\mathrm{obs}}$ of 2200 m a.s.l. We calculate expectation values for $\Delta\theta_{\mathrm{c}}/\theta_{\mathrm{c}}$ :

[TABLE]

Appendix I Derivation of the Variance of Ring Width from Variations of the Refractive Index with Altitude

A change in the emission altitude produces a relative change of the Cherenkov angle, according to

[TABLE]

where the density of air is approximated by $\rho_{0}\cdot\exp(-h/H_{0})$ , and $\rho_{0}\approx$ 1.225 kg m*-3*. We calculate expectation values for $\Delta\theta_{\mathrm{c}}/\theta_{\mathrm{c}}$ and for $(\Delta\theta_{\mathrm{c}}/\theta_{\mathrm{c}})^{2}$ :

[TABLE]

Appendix J Derivation of the Variance of Ring Width from Muon Energy Loss

[TABLE]

Then, analog to the previous section, we can calculate the variance of $\frac{\displaystyle\Delta\theta_{\mathrm{c}}(x)}{\displaystyle\theta_{\mathrm{c}}}$ .

[TABLE]

Bibliography153

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Aartsen et al. (2016) Aartsen, M. G., Abraham, K., Ackermann, M., et al., May 2016. Characterization of the atmospheric muon flux in icecube. Astropart Phys. 78, 1–27.
2Acero et al. (2010) Acero, F., Aharonian, F., Akhperjanian, A. G., et al., 2010. Localizing the VHE γ 𝛾 \gamma -ray source at the Galactic Centre. MNRAS 402, 1877–1882, ar Xiv:0911.1912.
3Acharya et al. (2013) Acharya, B. S., et al., 2013. Introducing the CTA concept. Astroparticle Physics 43, 3–18.
4Actis et al. (2011) Actis, M., et al., 2011. Design concepts for the Cherenkov Telescope Array CTA: an advanced facility for ground-based high-energy gamma-ray astronomy. Experimental Astronomy 32, 193–316, ar Xiv:1008.3703.
5Aguilar et al. (2015) Aguilar, J. A., Basili, A., Boccone, V., et al., Jan. 2015. Design, optimization and characterization of the light concentrators of the single-mirror small size telescopes of the Cherenkov Telescope Array. Astroparticle Physics 60, 32–40.
6Aguilar et al. (2014) Aguilar, J. A., et al., 2014. The Cherenkov Telescope Array Single-mirror Small Size Telescope project: status and prospect. In: Proc. SPIE Astronomical Telescopes + Instrumentation. Vol. 9145 of Ground-based and Airborne Telescopes V. p. 914508.
7Aharonian et al. (2004) Aharonian, F., et al., 2004. Calibration of cameras of the H.E.S.S. detector. Astropart.Phys. 22, 109–125, ar Xiv:astro-ph/0406658.
8Aharonian et al. (2006) Aharonian, F., et al., 2006. Observations of the Crab nebula with HESS. Astronomy & Astrophysics 457, 899–915.