# Fluorescence Lifetime Hong-Ou-Mandel Sensing

Ashley Lyons<EMAIL_ADDRESS>School of Physics and Astronomy,
University of Glasgow, Glasgow, G12 8QQ, UK  Vytautas Zickus  School of
Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK  Department
of Laser Technologies, Center for Physical Sciences and Technology, LT-10257,
Vilnius, Lithuania  Raúl Álvarez-Mendoza  School of Physics and Astronomy,
University of Glasgow, Glasgow, G12 8QQ, UK  Danilo Triggiani  School of
Mathematics and Physics, University of Portsmouth, Portsmouth, PO1 3QL, UK
Vincenzo Tamma  School of Mathematics and Physics, University of Portsmouth,
Portsmouth, PO1 3QL, UK  Institute of Cosmology and Gravitation, University of
Portsmouth, Portsmouth PO1 3FX, UK  Niclas Westerberg  School of Physics and
Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK  Manlio Tassieri  James
Watt School of Engineering, University of Glasgow, Glasgow, G12 8QQ, UK
Daniele Faccio<EMAIL_ADDRESS>School of Physics and Astronomy,
University of Glasgow, Glasgow, G12 8QQ, UK

###### Abstract

Fluorescence Lifetime Imaging Microscopy in the time domain is typically
performed by recording the arrival time of photons either by using electronic
time tagging or a gated detector. As such the temporal resolution is limited
by the performance of the electronics to $100$’s of picoseconds. Here, we
demonstrate a fluorescence lifetime measurement technique based on photon-
bunching statistics with a resolution that is only dependent on the duration
of the reference photon or laser pulse, which can readily reach the $1$-$0.1$
picosecond timescale. A range of fluorescent dyes having lifetimes spanning
from $1.6$ to $7$ picoseconds have been here measured with only $\sim 1$
second measurement duration. We corroborate the effectiveness of the technique
by measuring the Newtonian viscosity of glycerol/water mixtures by means of a
molecular rotor having over an order of magnitude variability in lifetime,
thus introducing a new method for contact-free nanorheology. Accessing
fluorescence lifetime information at such high temporal resolution opens a
doorway for a wide range of fluorescent markers to be adopted for studying yet
unexplored fast biological processes, as well as fundamental interactions such
as lifetime shortening in resonant plasmonic devices.

## Introduction

Fluorescence Lifetime Imaging Microscopy (FLIM) measures the exponential decay
time of fluorophores excited by an ultrafast source. FLIM is used across the
bio-imaging community to provide key information about local biological
environments as it can be dependent on pH, temperature, viscosity, and
chemical concentrations. The lifetimes of commonly used fluorophores are
within the range of $100$’s of picoseconds to nanoseconds. However, lifetimes
can shorten down to $10$’s of femtoseconds depending on the nature of the
fluorophore [1]. As such, access to higher resolutions enables the use of
femtosecond scale fluorophores [2, 3], coupling to plasmonic devices which can
increase quantum yield [4, 5], and exploration of vibrational and rotational
modes that act on the femtosecond scale [6].
In the time domain, the lifetime is measured by a combination of a single
photon detector and an electronic timing mechanism such as either gating the
detector or by assigning each individual photon an arrival time through Time
Correlated Single Photon Counting (TCSPC). The resolution is then typically
limited in the $50-100$ picosecond scale or longer [7]. This is due, for
example, to the Instrument Response Function (IRF) of the detector, which is
limited by effects such as electron diffusion [8].
Nonlinear optical gating methods where the fluorescent signal is combined with
a laser pulse in a nonlinear frequency-conversion crystal, have demonstrated
femtosecond scale resolution for fluorescence lifetime measurements [9, 10,
11, 12, 13]. However, the overall low efficiency of the nonlinear effects
requires measurements to be performed with high power lasers and to have a
limited wavelength flexibility for a given system, due to phase matching
considerations and choice of nonlinear crystals. For instance, Kerr-gating is
an alternative nonlinear optical method that can achieve femtosecond-scale
resolutions, but with efficiencies reaching up to 50% at best [14, 15, 16].
Alternatively, streak cameras can achieve sub-picosecond resolutions albeit at
the cost of significant hardware cost and limited dynamic range, specifically
when operated at high temporal resolutions [17]. Similarly, Single Photon
Superconducting Nanowire Detectors (SNSPDs) have recently been reported with
few-picosecond temporal resolutions [18, 19]
Hong-Ou-Mandel (HOM) is a two-photon interference effect whereby two
indistinguishable photons that meet at a beamsplitter will leave via the same
output port [20]. This photon bunching leads to a characteristic “HOM dip” in
the photon coincidence counts measured by the second order correlation
function $g^{(2)}(\tau)$, as a function of the relative delay time, $\tau$.
Originally proposed for measuring sub-picosecond time intervals, the use of
HOM interferometry for photon time-of-flight sensing has since been adopted
for measuring optical delays due to birefringence [21], and path-delay
measurements [22, 23, 24] with attoseconds time-scale resolution [25, 26].
Recently, HOM interferometry has been applied to molecular spectroscopy,
demonstrating an ability to measure dephasing times on the $100$ fs time-scale
from an organic dye in solution [27]. Notably, two-photon HOM interference can
be observed between completely independent light sources as first demonstrated
by Rarity et al. [28] and can even be extended to incoherent light sources.
This is well illustrated by the work performed by Deng et al. [29] showing HOM
interference between photons generated by a quantum dot and photons collected
from the Sun.

Figure 1: Schematic layout of the experimental setup consisting of a Mach-
Zehnder interferometer within one arm of which we place the fluorescent
sample. Not shown in the figure are spectral filters to remove the 1040 nm
pump before the microscope objective and various lenses for re-sizing the
beam. Red lines indicate the reference and emitted fluorescence beam-paths at
660 nm, orange shows the 1040 nm pump beam for the SHG, and green shows the
resulting 520 nm excitation beam. An illustration of the measured second order
correlation function, $g^{(2)}(\tau)$, is shown in the top right inset from
which the lifetime is estimated.

Fluorescence lifetime can also be measured by monitoring the photon anti-
bunching when pumped by two femtosecond pulses with a relative temporal delay
[30], this approach is however restricted to isolated single molecules. Photon
bunching has also been used for higher densities of molecules [3, 31] with
schemes that produce an auto-correlation of the fluorescence signal.
Recovering the actual lifetime from the autocorrelation requires the solution
of an ill-posed inverse problem [32].
In this work, we propose HOM interference between fluorescence photons and
secondary reference photons and demonstrate experimentally the application in
the classical case where the reference photon is replaced by an attenuated
laser pulse. Our approach, which we call Fluorescence Lifetime Hong-Ou-Mandel
(FL-HOM), offers the benefit of measuring directly the fluorescence lifetime
signal by using a reference photon or laser pulse that is much shorter than
the fluorescence lifetime. FL-HOM is able to operate at low photon numbers and
does not rely on any optical nonlinear components. We demonstrate the
potential of this approach by characterising a range of fluorescence lifetimes
in the picosecond range, including those from so-called molecular rotors that
are sensitive to the environment viscosity. Previous measurements with these
rotors were limited to high viscous systems (i.e., long lifetimes). With FL-
HOM, their lifetime can be resolved in the ps regime, thus enabling
‘contactless nanorheology’ measurements of very low viscosity fluids.

## Results

Experimental Layout. The setup is made of a Mach-Zehnder interferometer using
two outputs from a single femtosecond laser source running at an $80$ MHz
repetition rate (Coherent Chameleon Discovery NX). The first output, which we
label “Excitation Pulse” (EP), has a fixed wavelength at $1040$ nm and $140$
fs pulse duration. This is frequency-doubled in a BBO crystal to $520$ nm and
focused on to the fluorescent sample using a Nikon $60\times$, $0.7$ NA
microscope air objective. The fluorescence is collected along the same beam
path in a confocal geometry, and spatially filtered by coupling into a single
mode polarisation maintaining fiber. An interference filter is placed before
the fiber so as to select a spectral band around $660$ nm with two different
bandwidth options (details below). The second laser output, which we label
“Reference Pulse” (RP), synchronised to the EP, and is tuned to the
fluorescent $660$ nm wavelength and coupled into a single mode fiber with the
same spectral filtering as the fluorescence arm. The fiber coupler is mounted
on a translation stage so as to control the relative path delay between the
two interferometer arms. Both interferometer arms are then directed to a
fiber-coupled $50:50$ beamsplitter and photon correlations are measured
between the two beamsplitter outputs by a pair of Single Photon Avalanche
Diode (SPAD) detectors and a commercial Time Correlated Single Photon Counting
(TCSPC) unit (ID Quantique ID801).

Lifetime Sensing. The intensity profile of the fluorescent emission is given
by the following convolution integral:

$I_{\text{FP}}(t)=I_{\text{ex}}(t)*F_{\text{fl}}(t),$ (1)

where $I_{\text{ex}}(t)$ is temporal intensity distribution of the excitation
pulse, $F_{\text{fl}}(t)=Ae^{-t/\mu}$ is the fluorescent single exponential
decay response of the sample, with amplitude $A$ and decay constant (i.e.
fluorescence lifetime) $\mu$. Here, we focus on the limit where the excitation
pulse is of short duration with respect to $\mu$ such that
$I_{\text{FP}}(t)\simeq F_{\text{fl}}(t)$. The general form of the normalised
second order intensity correlation function is:

$g^{(2)}(t^{\prime})=\frac{\left<I_{1}(t)I_{2}(t+t^{\prime})\right>}{\left<I_{1}(t)\right>\left<I_{2}(t+t^{\prime})\right>},$
(2)

where $I_{1}(t)$ is the intensity measured at one detector at time $t$ and
$I_{2}(t)$ is the intensity measured at the second detector at a time
$t+t^{\prime}$. From Eq. (2), it is possible to calculate the normalised
number of coincident photon pair events (which we define as those occurring
within some time window of the same detection time i.e. $t^{\prime}=0$) as a
function of the optical delay, $\tau$, between the two arms of the
interferometer controlled by a translation stage (see Fig. 1):

$C(\tau)\sim 1-4C_{0}(\sigma)e^{-\tau/\mu},$ (3)

where $\sigma$ is the temporal duration of the reference photon or laser
pulse, and $C_{0}$ is a constant that encapsulates the visibility of the
interference depending on the amplitude of the two pulses, their durations,
and the reflection and transmission constants of the beamsplitter. Here we
have used that the reference pulse is much shorter than the lifetime (i.e.
$\sigma\ll\mu$), and that $\tau>0$. A full derivation of $C(\tau)$ can be
found in the SI.

From Eq. (3) it is possible to observe that the lifetime can be measured
directly from the number of coincident events. We also note that the
resolution can be controlled by changing the pulse duration of the reference
and/or excitation pulses, which can be achieved by using spectral bandpass
filters before the detectors.

Figure 2: (a) Fluorescence lifetime curve for the 4-DASPI sample (the curve
has been inverted and normalised so as to appear more familiar compared to
standard FLIM curves). Blue shows experimental data (2 s acquisition per delay
point) and red depicts a fit from our model described in the Supplementary
Material using Equation S11. (b) The same measurements with a 3.5 s total
acquisition time, retrieving the same approximate lifetime. based on fitting
to the points shown in green. (c) Shows the statistical distribution of
retrieved lifetimes for N=260 measurements (3.5 s acquisition time each),
indicating a mean of 6.8 ps and a standard deviation of 0.7 ps. (d) The
standard deviation of the lifetime estimation as a function of the total
measurement acquisition time. Blue points are the experimental data points and
the red line illustrates the expected reduction in uncertainty proportional to
the square root of the measurement time.

Measurements of picosecond fluorescence decay. Figure 2(a) shows an example of
a normalised FL-HOM trace for a sample consisting of $4$-DASPI dye [33]. To
ensure a high level of indistinguishability between the two interferometer
arms, a narrowband $0.6$ nm FWHM spectral filter (Alluxa $660$-$0.5$-OD$6$)
was used in each arm. This has the overall effect of increasing the
interference visibility at the expense of reducing the achievable temporal
resolution. Indeed, a wider bandwidth filter of $10$ nm showed, under the same
conditions, a maximum visibility on the order of ten times smaller dependent
on the lifetime (see below). Typical photon count rates were of the order
$10^{6}$ photons/second with a photon coincidence count rate of $\sim 14\times
10^{3}$ photons/second outside of the HOM dip. The data shown in Fig. 2(a) was
obtained with a long acquisition time of $2$ s per delay point, for a total
measurement time of approximately $40$ minutes, and the number of coincident
events was normalised by the total number of photon counts to account for
fluctuations in laser power and fiber coupling. The step size of the
translation stage equates to an optical delay of $16.7$ fs. Figure 2(a) also
shows a fit to our model without any approximations about the length of the
reference pulse (i.e. before Eq. (3)) and the measured IRF, from which we
recover a lifetime of $7.22\pm 0.04$ ps. Previous work in [33] measured a
lifetime of 11 ps for 4-DASPI in water. Whilst this is relatively close to our
measurement, we attribute our shorter lifetime to potential differences in the
temperature, pH, and emission wavelength chosen, all of which can affect the
lifetime. To help validate our approach, we also conducted a measurement of
4-DASPI in ethanol which is known to produce a longer lifetime. From this we
observed a lifetime of 65 $\pm$1.2 ps which is in close agreement with
previous measurements of 62 ps [34] (see Supplementary Material). In order to
investigate and ensure the capability of our approach, we significantly over-
sampled the data and purposely measured for long times, thus resulting in a
much lower noise level than a typical time domain FLIM measurement. As a means
of comparison, Fig. 2(b) shows a similar measurement, but using a total
acquisition time of only $3.5$ s. By using a Markov Chain Monte Carlo approach
for the lifetime retrieval, our model still returns a reliable lifetime of
$6.8$ ps. A full distribution of lifetimes derived under these conditions
(i.e., $3.5$ s total acquisition time, 260 individual measurements) is shown
in Fig. 2(c), which has a standard deviation of $\sigma_{\mu}=0.7$ ps. In Fig.
2(d) we report $\sigma_{\mu}$ versus total measurement time. This data were
obtained by taking multiple short exposures of $50$ ms for each optical delay
point and summing them to obtain longer acquisition times. This was then
repeated over multiple iterations to statistically evaluate the uncertainty of
the measured lifetime (typical measurements and fits for each acquisition time
can be found in the SI). In Fig. 2(d) the red line indicates a gradient of
-0.5 i.e. a line proportional to the square root of the total measurement time
which we expect to observe given that we have a set of independent and
statistically equivalent measurements. It is also useful to determine how the
number of coincident events themselves scale as a function of the measurement
time. We evaluate this (see SM) finding that the Signal to Noise Ratio
improves proportionally to $\nicefrac{{N}}{{\sqrt{1+2N}}}$. For $N>>1$ the
$\textrm{SNR}\sim\sqrt{N}/\sqrt{2}$, indicating that coincidence counting
brings a penalty of $\sqrt{2}$ in the SNR compared to direct intensity
measurements.

Figure 3: Lifetime estimation of different samples. The natural log of the
(inverted) normalised coincidence events is shown for comparison with typical
FLIM data. Crosses indicate measured data points and solid lines show a linear
fit to the exponential tail. Figure 4: (a) The effect of spectral filtering
for two different bandwidths is shown for Allura Red fluorescence. The wider
bandwidth results in a higher temporal resolution measurement at the expense
of reducing the interferometric visibility. The signal from the 10 nm
broadband filter has been scaled by a factor of 10x for comparison. The black
dashed line shows the normalised single photon counts at one of the detectors
for the 0.6 nm filter measurement. (b) Model of the visibility as a function
of the ratio between the pulse duration and fluorescence lifetime.

To further demonstrate the FL-HOM approach we repeated the measurements for
two alternative fluorescent dyes i.e., Allura Red and Pinacyanol Iodide.
Figure 3 shows the related data, where we plot the logarithm of the normalised
coincident events. From the analysis of these measurements we retrieve
lifetimes of $6.23$ (4-DASPI), $1.60$ (Allura Red) and $3.66$ (Pinacyanol
Iodide) picoseconds. As the lifetimes measured here are close to the effective
IRF (see below), we validated our shortest lifetime measurement using the
Allura Red against our MCMC retrieval. From this, we found a lifetime of 1.51
ps indicating that simple tail fits provide a reasonable lifetime estimation
at these short timescales. We also observed a variation in the lifetime by
changing the temperature using the $4$-DASPI sample, showing a lifetime range
of $5$-$8.25$ picoseconds (see SI).
Interference visibility. The maximum achievable interference visibility is
dependent on the overlap integral between the fluorescence and reference
pulses. Figure 4(a) shows the normalised number of coincident events for the
Allura Red sample with both a $0.6$ nm and a $10$ nm bandpass filter, keeping
all other conditions fixed. We measure the IRFs with these two filters to be
2.08 ps and 0.17 ps FWHM respectively (See Methods). In the case of the
broadband filter, the temporal overlap between the two pulse always remains
low due to the relatively large difference between the $140$ fs reference
pulse and the much longer $1.6$ ps duration of the lifetime. With the
narrowband filter on the other hand, the maximum temporal overlap becomes much
larger as a 0.6 nm bandwidth corresponds to a transform limited pulse of order
1 ps duration. We also show (dashed black line) the normalised single photon
counts for one of the detectors in our setup under conditions where we
achieved maximum HOM interference visibility. No interference fringes can be
observed confirming that we are operating in a regime where there is no first
order interference. In Fig. 4(b) we further model how the HOM interference
visibility varies with the reference photon or laser pulse duration by
evaluating Eq. (3). We observe that there is a trade-off between the temporal
resolution and the FL-HOM visibility with optimal visibility being achieved
when the pulse duration is $0.7\times$ the lifetime. For longer/shorter pulse
durations the pulses become less indistiguishable in the temporal domain
thereby reducing the visibility.
Molecular Rotor Lifetime as a Probe for Viscosity. As an example of FL-HOM
capability, we show that lifetime measurements of rotor molecules with
picosecond resolution enables contact-free nanorheology over a broad range of
viscosities. Rheology is the study of flow of matter, and its applications
have revealed valuable information on biological specimens [35] at different
time and length scales, down to intracellular dynamics at microscales [36].
Recent work has also revealed the need for contact-free nanoscale rheology to
better understand the impact of viscosity on carbon capture and microbe
populations in the oceans [37]. Nanorheology has been enabled in cells using
atomic force microscopy, where a metal tip is placed in ‘contact’ with the
sample [38]. However, such an approach is not viable for intracellular or
remote measurements.
Rotor molecules interact with the surrounding medium through their rotational
degrees of freedom. The formation of twisted intramolecular charge transfer
(TICT) states upon photoexcitation leads to competing de-excitation pathways:
fluorescence emission and non-radiative de-excitation from the TICT state.
Since TICT formation is viscosity-dependent, the emission intensity and
lifetime of molecular rotors depends on the solvent’s viscosity [39, 40].
However, intensity-based measurements are influenced by (i) the optical
properties of the surrounding environment, (ii) the dye concentration and
(iii) solvent–dye interactions, thus requiring a calibration of the
measurement system to a specific solvent. Lifetime measurements instead,
provide the opportunity of a calibration-free approach [41]. However, previous
studies have been limited to the measurement of very high viscous fluids
(i.e., nearly pure glycerol, which is circa thousand times more viscous than
water) as only in this regime the rotor molecules (such as 4-DASPI) have
sufficiently long lifetimes for these to be measured and/or resolved with
$\sim 100$ ps of standard FLIM approaches. For instance, 4-DASPI has been used
for viscosity measurements by using both intensity [42] and sub-nanosecond
lifetime measurements [43] for which the shortest lifetime measurements were
of order of $700$ ps.
We used 4-DASPI dye in water and glycerol/water mixtures to explore a range of
low viscosity values (i.e., of the order of mPa$\cdot$s), whilst maintaining a
constant pH and temperature (see Methods). The viscosity within cellular
environments can lie between 1 - 400 mPa$\cdot$s or even higher [44, 45, 46].
Here we focus on the lower 1-5% of this range to demonstrate the high
sensitivity of the technique. The viscosity of each mixture was measured by
means of conventional bulk-rheology measurements performed with a commercial
rheometer (see Methods). Figure 5 shows FL-HOM lifetime measurements as a
function of viscosity. We observe more than an order of magnitude change in
lifetime, with no appreciable change in the photon count levels thus revealing
that FL-HOM lifetime measurement is a much more sensitive metric than
intensity measurements at such low viscosity values. Interestingly, at lower
viscosity values, it is only the $2.8$ ps reference pulse duration after the
filter (our effective IRF) that is limiting the viscosity measurement to a
value as low as $0.0065$ mPa$\cdot$s, which is more than two orders of
magnitude lower than that of water (which is $\eta\sim 1$mPa$\cdot$s).
Moreover, FL-HOM returns a resolution of $<1\%$ and is thus competitive to
more established microrheology measurement techniques such as multiple
particle tracking or optical tweezers [47, 48]. At relatively high viscosity
values, the noise in the decay function measurement at long times is the
limiting factor affecting the measurements, yet the method still returns a
$\sim 2-4$% error in viscosity measurements. Notably, for longer lifetimes
than those explored here, one would simply switch to using the standard time-
correlated-single-photon techniques by analysing the data at just one
detector, with no changes required in the setup.

Figure 5: FL-HOM lifetime measurement of 4-DASPI as a viscosity probe.
Glycerol is added to the solution increase the viscosity, starting from that
of pure water. More than an order of magnitude dynamic range is observed in
the lifetime with no appreciable change in the photon count rate, indicating a
high level of sensitivity to the viscosity. Uncertainties in the lifetime are
shown although for most point, the error bars are smaller than the data point
and thus not visible.

## Discussion

In Fig. 2(b) we demonstrate that our FL-HOM system can retrieve a lifetime
measurement as low as $6.8$ ps within $3.5$ s measurement duration; with
potentials to be significantly improved. The SPADs were kept at a maximum
count rate of $10^{6}$ photons per second to avoid pileup effects with a
probability of detecting a photon of only around 1$\%$ for each of the $80$
MHz rate laser pulses and a $0.016$% probability of detecting a correlated
pair event. Bearing in mind that the probability of measuring a coincidence
scales with the product of the photon numbers at each detector, decreasing the
repetition rate of the laser whilst keeping the average photon number constant
can lead to an increase in the number of measured coincidences. For example,
if a repetition rate of $0.1$ MHz is used with an average (photon starved)
$10$% probability to measure a photon at each detector from each pulse, the
measurement time can be reduced by a factor of $\sim 10$x implying that the
measurements in Fig. 2(b) could be obtained in sub-second timescales.
The high temporal resolution of FL-HOM provides the possibility to study
faster fluorescent dynamics such as those from dyes with femtosecond-scale
lifetimes [2, 3]. This allows for a wider range of markers to be used for
FLIM, thereby also increasing the dynamic range. Individual molecular
processes can also be investigated with transitions relating to the
vibrational modes commonly falling in the femtosecond regime [6]. Although
longer lifetime fluorophores generally have a greater quantum yield than
shorter lifetime probes, interaction with plasmonic modes have also shown a
strong increase in the fluorescence intensity at the expense of shortening the
lifetime [4, 5]. We also note that as our method directly measures the
temporal distribution of the fluorescent light, it allows also for the
measurement of bi-exponential and more generalised fluorescent decays. It is
also not restricted to measurements of single molecules and in this way, is
complementary to techniques such as those shown in [30].
It is worth mentioning that our approach is related, but somewhat different,
to ‘fluorescence lifetime correlation spectroscopy’ (FLCS) [49]. In FLCS,
rapid fluctuations of the fluorescence intensity are recorded in time in a
Hanbury-Brown Twiss interferometer, i.e. a single beamsplitter is placed
_after_ the sample and correlations between the two outputs are recorded
without the presence of any additional reference photons. These correlations
can be resolved within the duration of the fluorescence decay but are still
limited in resolution by the IRF of the SPAD sensors, similarly to standard
FLIM approaches.
Finally, our experiments are implemented with a classical reference laser
pulse. One could use a single photon generated for example by a parametric
downconversion (PDC) or a quantum dot source. It has previously been shown
that the visibility of the HOM interference between a single photon source
and, for example, an attenuated laser [28] can approach unity where an
equivalent measurement with classical light would be limited to a $g^{(2)}$ of
$0.5$. Furthermore, a PDC source provides greater flexibility in terms of
tuning the reference wavelength to the fluorescence emission by operating at
non-degenerate photon pair wavelengths controlled through the phase-matching
of the generation crystal.

## Methods

Fluorescent dye sample preparation. 4-DASPI (Sigma-Aldrich 336408) and Allura
Red AC (Sigma-Aldrich 458848) were dissolved in purified water at
concentrations of $\approx$2.7 3mM and $\approx$5.04mM respectively, and the
Pinacyanol iodide (Alfa Aesar H31540) was dissolved in MeOH (Sigma-Aldrich
1.06002) at the concentration of $\approx$0.12 mM. The samples were measured
at room temperature, in 10 mm path length UV fused quartz cuvettes (Thorlabs
CV10Q35E). We note that in all cases, the probability of re-absoprtion in
these samples is low due to the relatively large Stokes shift. We also confirm
that there is a low probability of interaction between fluorophores whilst in
the excited state. A discussion of this can be found in the SM.

Glycerol/Water Mixtures. Samples for viscosity measurements were prepared by
varying the weight concentration of Glycerol (Sigma-Aldrich G5516) in purified
water and 4-DASPI (Sigma-Aldrich 336408) dye. 4-DASPI was kept at constant
concentration of $0.1$ wt% in all samples. The glycerol/water-4-DASPI
solutions were then titrated to a pH of $\approx 7.0$ while continuously
mixed.

Effective IRF measurement. To estimate the IRF, we perform a linear auto-
correlation of the 660 nm reference pulse. A 50:50 beamsplitter is placed
before the single mode fiber seen in Fig. 1(a) and the reflection is coupled
into another fiber. This then replaces the input from the fluorescence arm
into the same fiber coupled beamsplitter to ensure spatial mode matching. The
coincident photon counts at one of the detectors is measured as a function of
the delay stage position and a low pass filter is applied to remove any
residual first order coherence fringes and retrieve the envelope function.
From this, we measure a FWHM of the 2.08 ps and 0.17 ps for the 0.6 nm and 10
nm bandpass filters respectively.

Markov Chain Monte Carlo Lifetime Estimation for Short Acquisition Times. The
lmfit package in Python was first used to fit the experimental data to our
model by solving the non-linear Least Squares regression problem with the
Levenberg-Marquardt algorithm [50]. The parameter space around the best-fit
solution was explored by means of the Affine Invariant Markov Chain Monte
Carlo (MCMC) Ensemble sampler that is implemented in lmfit via the
Minimizer.emcee() method, which allowed us to refine the initial Least Squares
estimate. The sampling rate of the lifetime curves that were deconvolved using
this approach was one-tenth of that of the curves obtained using long
acquisition times, and simultaneously a smaller delay region was sampled where
the experimental noise is lowest to reduce total acquisition time. For the
lifetime curve shown in Fig. 2, this implies that the peak was sampled with
only 35 delay points, resulting in acquisition time of 3.5 s.

Conventional Rheology. Bulk rheology measurements were performed on each
explored sample to determine the ground-truth viscosity. These were performed
by using a stress controlled rheometer (Anton Paar Physica MCR $302$
Instruments) equipped with a cone-plate measuring system (CP60-1-SN42255). The
temperature was controlled by means of a Peltier system connected to a water
bath. The fluids’ viscosity was measured by performing a flow curve test at
shear rates varying from $50$ s-1 to $500$ s-1. For each sample the viscosity
had a constant value within the explored range of shear rates.

## Data Availability

All data used to produce the plots shown in this article is available at DOI
10.5525/gla.researchdata.1524.

## References

## References

  * [1] Mikhail Y. Berezin and Samuel Achilefu.  Fluorescence lifetime measurements and biological imaging.  Chemical Reviews, 110:2641–2684, 5 2010.
  * [2] Hideki Kandori, Yuko Katsuta, Masayoshi Ito, and Hiroyuki Sasabe.  Femtosecond fluorescence study of the rhodopsin chromophore in solution.  Journal of the American Chemical Society, 117:2669–2670, 3 1995\.
  * [3] Motoyoshi Baba, Ying Li, and Masahiro Matsuoka.  Intensity interference of ultrashort pulsed fluorescence.  Physical Review Letters, 76:4697–4700, 6 1996.
  * [4] O. L. Muskens, V. Giannini, J. A. Sánchez-Gil, and J. Gómez Rivas.  Strong enhancement of the radiative decay rate of emitters by single plasmonic nanoantennas.  Nano Letters, 7:2871–2875, 9 2007.
  * [5] Anika Kinkhabwala, Zongfu Yu, Shanhui Fan, Yuri Avlasevich, Klaus Müllen, and W. E. Moerner.  Large single-molecule fluorescence enhancements produced by a bowtie nanoantenna.  Nature Photonics, 3:654–657, 11 2009.
  * [6] Ahmed H. Zewail.  Femtochemistry: atomic-scale dynamics of the chemical bond.  The Journal of Physical Chemistry A, 104:5660–5694, 6 2000.
  * [7] Claudio Bruschini, Harald Homulle, Ivan Michel Antolovic, Samuel Burri, and Edoardo Charbon.  Single-photon avalanche diode imagers in biophotonics: review and outlook.  Light: Science & Applications, 8:87, 12 2019.
  * [8] G. Ripamonti and S. Cova.  Carrier diffusion effects in the time-response of a fast photodiode.  Solid-State Electronics, 28:925–931, 9 1985.
  * [9] Leslie A. Hallidy and Michael R. Topp.  Picosecond luminescence detection using type-ii phase-matched frequency conversion.  Chemical Physics Letters, 46:8–14, 2 1977.
  * [10] Jianhua Xu and Jay R. Knutson.  Chapter 8 ultrafast fluorescence spectroscopy via upconversion: Applications to biophysics.  In Fluorescence Spectroscopy, volume 450 of Methods in Enzymology, pages 159–183. Academic Press, 2008.
  * [11] Simin Cao, Haoyang Li, Zenan Zhao, Sanjun Zhang, Jinquan Chen, Jianhua Xu, Jay R. Knutson, and Ludwig Brand.  Ultrafast fluorescence spectroscopy via upconversion and its applications in biophysics.  Molecules, 26:211, 1 2021.
  * [12] M. MÜLLER, R. GHAUHARALI, K. VISSCHER, and G. BRAKENHOFF.  Double-pulse fluorescence lifetime imaging in confocal microscopy.  Journal of Microscopy, 177:171–179, 2 1995.
  * [13] A. H. Buist, M. Müller, E. J. Gijsbers, G. J. Brakenhoff, T. S. Sosnowski, T. B. Norris, and J. Squier.  Double‐pulse fluorescence lifetime measurements.  Journal of Microscopy, 186:212–220, 6 1997.
  * [14] B. Schmidt, S. Laimgruber, W. Zinth, and P. Gilch.  A broadband kerr shutter for femtosecond fluorescence spectroscopy.  Applied Physics B, 76:809–814, 7 2003.
  * [15] Jolie C. Blake, Jesus Nieto-Pescador, Zhengxin Li, and Lars Gundlach.  Ultraviolet femtosecond kerr-gated wide-field fluorescence microscopy.  Optics Letters, 41:2462, 6 2016.
  * [16] Zhenqiang Huang, Wenjiang Tan, Jinhai Si, Shijia Zeng, Zhen Kang, and Xun Hou.  High-transmittance femtosecond optical kerr gate with double gate pulses based on birefringence effect.  Photonics, 9:71, 1 2022.
  * [17] A. Campillo and S. Shapiro.  Picosecond streak camera fluorometry - a review.  IEEE Journal of Quantum Electronics, 19:585–603, 4 1983.
  * [18] Boris Korzh, Qing-Yuan Zhao, Jason P. Allmaras, Simone Frasca, Travis M. Autry, Eric A. Bersin, Andrew D. Beyer, Ryan M. Briggs, Bruce Bumble, Marco Colangelo, Garrison M. Crouch, Andrew E. Dane, Thomas Gerrits, Adriana E. Lita, Francesco Marsili, Galan Moody, Cristián Peña, Edward Ramirez, Jake D. Rezac, Neil Sinclair, Martin J. Stevens, Angel E. Velasco, Varun B. Verma, Emma E. Wollman, Si Xie, Di Zhu, Paul D. Hale, Maria Spiropulu, Kevin L. Silverman, Richard P. Mirin, Sae Woo Nam, Alexander G. Kozorezov, Matthew D. Shaw, and Karl K. Berggren.  Demonstration of sub-3 ps temporal resolution with a superconducting nanowire single-photon detector.  Nature Photonics, 14:250–255, 4 2020.
  * [19] Robert H. Hadfield.  Superfast photon counting.  Nature Photonics, 14:201–202, 4 2020.
  * [20] C. K. Hong, Z. Y. Ou, and L. Mandel.  Measurement of subpicosecond time intervals between two photons by interference.  Physical Review Letters, 59:2044–2046, 1987.
  * [21] D. Branning, A. L. Migdall, and A. V. Sergienko.  Simultaneous measurement of group and phase delay between two photons.  Physical Review A - Atomic, Molecular, and Optical Physics, 62:063808–063801, 2000.
  * [22] Daniel Giovannini, Jacquiline Romero, Václav Potoček, Gergely Ferenczi, Fiona Speirits, Stephen M Barnett, Daniele Faccio, and Miles J Padgett.  Spatially structured photons that travel in free space slower than the speed of light.  Science (New York, N.Y.), 347:857–60, 2015.
  * [23] Michał Jachura and Radosław Chrapkiewicz.  Shot-by-shot imaging of hong–ou–mandel interference with an intensified scmos camera.  Optics Letters, 40:1540, 4 2015.
  * [24] Ashley Lyons, Thomas Roger, Niclas Westerberg, Stefano Vezzoli, Calum Maitland, Jonathan Leach, Miles J. Padgett, and Daniele Faccio.  How fast is a twisted photon?  Optica, 5:682, 6 2018.
  * [25] Ashley Lyons, George C Knee, Eliot Bolduc, Thomas Roger, Jonathan Leach, Erik M Gauger, and Daniele Faccio.  Attosecond-resolution hong-ou-mandel interferometry.  Science Advances, 4:eaap9416, 5 2018.
  * [26] Danilo Triggiani, Giorgos Psaroudis, and Vincenzo Tamma.  Ultimate quantum sensitivity in the estimation of the delay between two interfering photons through frequency-resolving sampling.  Physical Review Applied, 19:044068, 4 2023.
  * [27] Audrey Eshun, Bing Gu, Oleg Varnavski, Shahaf Asban, Konstantin E. Dorfman, Shaul Mukamel, and Theodore Goodson.  Investigations of molecular optical properties using quantum light and hong–ou–mandel interferometry.  Journal of the American Chemical Society, 143:9070–9081, 6 2021\.
  * [28] J. G. Rarity, P. R. Tapster, and R. Loudon.  Non-classical interference between independent sources.  Journal of Optics B: Quantum and Semiclassical Optics, 7, 2005.
  * [29] Yu-Hao Deng, Hui Wang, Xing Ding, Z.-C. Duan, Jian Qin, M.-C. Chen, Yu-Ming He, Yu-Ming He, Jin-Peng Li, Yu-Huai Li, Li-Chao Peng, E. S. Matekole, Tim Byrnes, C. Schneider, M. Kamp, Da-Wei Wang, Jonathan P. Dowling, Sven Höfling, Chao-Yang Lu, Marlan O. Scully, and Jian-Wei Pan.  Quantum interference between light sources separated by 150 million kilometers.  Physical Review Letters, 123:080401, 8 2019.
  * [30] J. Schedlbauer, P. Wilhelm, L. Grabenhorst, M. E. Federl, B. Lalkens, F. Hinderer, U. Scherf, S. Höger, P. Tinnefeld, S. Bange, J. Vogelsang, and J. M. Lupton.  Ultrafast single-molecule fluorescence measured by femtosecond double-pulse excitation photon antibunching.  Nano Letters, 20:1074–1079, 2 2020.
  * [31] Ying Li, M. Baba, and M. Matsuoka.  Femtosecond measurement of fluorescence by two-photon interference.  Physical Review A, 55:3177–3183, 4 1997.
  * [32] A. J. Devaney.  The inverse problem for random sources.  Journal of Mathematical Physics, 20:1687–1691, 8 1979.
  * [33] Jiho Kim and Minyung Lee.  Excited-state photophysics and dynamics of a hemicyanine dye in aot reverse micelles.  The Journal of Physical Chemistry A, 103:3378–3382, 5 1999.
  * [34] W. Sibbett and J.R. Taylor.  Passive mode locking in the blue spectral region.  Optics Communications, 46:32–36, 6 1983.
  * [35] Leandro G Rizzi and Manlio Tassieri.  Microrheology of biological specimens.  Encyclopedia of Analytical Chemistry: Applications, Theory and Instrumentation, pages 1–24, 2006.
  * [36] Claude Verdier, Jocelyn Etienne, Alain Duperray, and Luigi Preziosi.  Review: Rheological properties of biological materials.  Comptes Rendus Physique, 10:790–811, 11 2009.
  * [37] Òscar Guadayol, Tania Mendonca, Mariona Segura-Noguera, Amanda J. Wright, Manlio Tassieri, and Stuart Humphries.  Microrheology reveals microscale viscosity gradients in planktonic systems.  Proceedings of the National Academy of Sciences, 118, 1 2021.
  * [38] Pablo D. Garcia, Carlos R. Guerrero, and Ricardo Garcia.  Nanorheology of living cells measured by afm-based force–distance curves.  Nanoscale, 12:9133–9143, 4 2020.
  * [39] Mark A Haidekker and Emmanuel A Theodorakis.  Environment-sensitive behavior of fluorescent molecular rotors.  Journal of Biological Engineering, 4:11, 2010.
  * [40] Klaus Suhling, James A. Levitt, Pei-Hua Chung, Marina. K. Kuimova, and Gokhan Yahioglu.  Fluorescence lifetime imaging of molecular rotors in living cells.  Journal of Visualized Experiments, 2 2012.
  * [41] James A. Levitt, Marina K. Kuimova, Gokhan Yahioglu, Pei-Hua Chung, Klaus Suhling, and David Phillips.  Membrane-bound molecular rotors measure viscosity in live cells via fluorescence lifetime imaging.  The Journal of Physical Chemistry C, 113:11634–11642, 7 2009.
  * [42] Alice Briole, Thomas Podgorski, and Bérengère Abou.  Molecular rotors as intracellular probes of red blood cell stiffness.  Soft Matter, 17:4525–4537, 5 2021.
  * [43] Federico Caporaletti, Marius R. Bittermann, Daniel Bonn, and Sander Woutersen.  Fluorescent molecular rotor probes nanosecond viscosity changes.  The Journal of Chemical Physics, 156:201101, 5 2022.
  * [44] Xiaojun Peng, Zhigang Yang, Jingyun Wang, Jiangli Fan, Yanxia He, Fengling Song, Bingshuai Wang, Shiguo Sun, Junle Qu, Jing Qi, and Meng Yan.  Fluorescence ratiometry and fluorescence lifetime imaging: Using a single molecular sensor for dual mode imaging of cellular viscosity.  Journal of the American Chemical Society, 133:6626–6635, 5 2011\.
  * [45] Marina K. Kuimova.  Mapping viscosity in cells using molecular rotors.  Physical Chemistry Chemical Physics, 14:12671, 10 2012.
  * [46] Tianyu Liu, Xiaogang Liu, David R. Spring, Xuhong Qian, Jingnan Cui, and Zhaochao Xu.  Quantitatively mapping cellular viscosity with detailed organelle information via a designed pet fluorescent probe.  Scientific Reports, 4:5418, 6 2014.
  * [47] Manlio Tassieri.  Microrheology with optical tweezers: peaks & troughs.  Current Opinion in Colloid & Interface Science, 43:39–51, 2019\.
  * [48] John A. McGlynn, Nan Wu, and Kelly M. Schultz.  Multiple particle tracking microrheological characterization: Fundamentals, emerging techniques and applications.  Journal of Applied Physics, 127:201101, 5 2020.
  * [49] Arindam Ghosh, Narain Karedla, Jan Christoph Thiele, Ingo Gregor, and Jörg Enderlein.  Fluorescence lifetime correlation spectroscopy: Basics and applications.  Methods, 140-141:32–39, 5 2018.
  * [50] Matthew Newville, Till Stensitzki, Daniel B. Allen, and Antonino Ingargiola.  LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python. https://doi.org/10.5281/zenodo.11813. Accessed 28/09/23.

## Acknowledgements

The authors acknowledge financial support from the Royal Academy of
Engineering under the Chairs in Emerging Technology and Research Fellowships
schemes, and the U.K. Engineering and Physical Sciences Research Council
(Grants No. EP/T002123/1, EP/T00097X/1). DF, AL and MT acknowledge financial
support from the U.K. Engineering and Physical Sciences Research Council
(Grants No. EP/X035905/1). NW wishes to acknowledge support from the Royal
Commission for the Exhibition of 1851. VZ received funding from European
Social Fund (project No 09.3.3-LMT-K-712-23-0132) under a grant agreement with
the Research Council of Lithuania. VT acknowledges support from the Air Force
Office of Scientific Research under award number FA8655-23-1-7046.

## Author Contributions Statement

The concept for the work was developed by AL, VZ, and DF. Experimental work
was undertaken by AL, VZ, RAM whilst DT, VT, and NW developed the theoretical
framework. MT provided valuable insights for the rheology sections of the work
as well as ground truth viscosity measurements. The project was supervised by
AL and DF.

## Competing Interests Statement

The authors declare no competing interests.


# Calibration of the Air Shower Energy Scale of the Water and Air Cherenkov
Techniques in the LHAASO experiment

F. Aharonian26,27, Q. An4,5, Axikegu20, L.X. Bai21, Y.X. Bai1,3, Y.W. Bao15,
D. Bastieri10, X.J. Bi1,2,3,
Y.J. Bi1,3, H. Cai23, J.T. Cai10, Z. Cao1,2,3,⋆ Z. Cao4,5, J. Chang16, J.F.
Chang1,3,4, X.C. Chang1,3,
B.M. Chen13, J. Chen21, L. Chen1,2,3, L. Chen18, L. Chen20, M.J. Chen1,3, M.L.
Chen1,3,4, Q.H. Chen20,
S.H. Chen1,2,3, S.Z. Chen1,3, T.L. Chen22, X.L. Chen1,2,3, Y. Chen15, N.
Cheng1,3, Y.D. Cheng1,3, S.W. Cui13,
X.H. Cui7, Y.D. Cui11, B.Z. Dai24, H.L. Dai1,3,4, Z.G. Dai15, Danzengluobu22,
D. della Volpe31, B. D’Ettorre Piazzoli28, X.J. Dong1,3, J.H. Fan10, Y.Z.
Fan16, Z.X. Fan1,3, J. Fang24, K. Fang1,3, C.F. Feng17, L. Feng16, S.H.
Feng1,3, Y.L. Feng16, B. Gao1,3, C.D. Gao17, Q. Gao22, W. Gao17, M.M. Ge24,
L.S. Geng1,3, G.H. Gong6, Q.B. Gou1,3, M.H. Gu1,3,4, J.G. Guo1,2,3, X.L.
Guo20, Y.Q. Guo1,3, Y.Y. Guo1,2,3,16, Y.A. Han14, H.H. He1,2,3, H.N. He16,
J.C. He1,2,3, S.L. He10, X.B. He11, Y. He20, M. Heller31, Y.K. Hor11, C.
Hou1,3, X. Hou25, H.B. Hu1,2,3, S. Hu21, S.C. Hu1,2,3, X.J. Hu6, D.H. Huang20,
Q.L. Huang1,3, W.H. Huang17, X.T. Huang17, Z.C. Huang20, F. Ji1,3, X.L.
Ji1,3,4, H.Y. Jia20, K. Jiang4,5, Z.J. Jiang24, C. Jin1,2,3, D. Kuleshov29, K.
Levochkin29, B.B. Li13, C. Li1,3, C. Li4,5, F. Li1,3,4, H.B. Li1,3, H.C.
Li1,3, H.Y. Li5,16, J. Li1,3,4, K. Li1,3, W.L. Li17, X. Li4,5, X. Li20, X.R.
Li1,3, Y. Li21, Y.Z. Li1,2,3, Z. Li1,3, Z. Li9, E.W. Liang12, Y.F. Liang12,
S.J. Lin11, B. Liu5, C. Liu1,3, D. Liu17, H. Liu20, H.D. Liu14, J. Liu1,3,
J.L. Liu19, J.S. Liu11, J.Y. Liu1,3, M.Y. Liu22, R.Y. Liu15, S.M. Liu16, W.
Liu1,3, Y.N. Liu6, Z.X. Liu21, W.J. Long20, R. Lu24, H.K. Lv1,3, B.Q. Ma9,
L.L. Ma1,3⋆, X.H. Ma1,3, J.R. Mao25, A. Masood20, W. Mitthumsiri32, T.
Montaruli31, Y.C. Nan17⋆, B.Y. Pang20, P. Pattarakijwanich32, Z.Y. Pei10, M.Y.
Qi1,3, D. Ruffolo32, V. Rulev29, A. Sáiz32, L. Shao13, O. Shchegolev29,30,
X.D. Sheng1,3, J.R. Shi1,3, H.C. Song9, Yu.V. Stenkin29,30, V. Stepanov29,
Q.N. Sun20, X.N. Sun12, Z.B. Sun8, P.H.T. Tam11, Z.B. Tang4,5, W.W. Tian2,7,
B.D. Wang1,3, C. Wang8, H. Wang20, H.G. Wang10, J.C. Wang25, J.S. Wang19, L.P.
Wang17, L.Y. Wang1,3, R.N. Wang20, W. Wang11, W. Wang23, X.G. Wang12, X.J.
Wang1,3, X.Y. Wang15, Y.D. Wang1,3, Y.J. Wang26⋆, Y.P. Wang1,2,3, Z.
Wang1,3,4, Z. Wang19, Z.H. Wang21, Z.X. Wang24, D.M. Wei16, J.J. Wei16, Y.J.
Wei1,2,3, T. Wen24, C.Y. Wu1,3, H.R. Wu1,3, S. Wu1,3, W.X. Wu20, X.F. Wu16,
S.Q. Xi20, J. Xia5,16, J.J. Xia20, G.M. Xiang2,18, G. Xiao1,3, H.B. Xiao10,
G.G. Xin23, Y.L. Xin20, Y. Xing18, D.L. Xu19, R.X. Xu9, L. Xue17, D.H. Yan25,
C.W. Yang21, F.F. Yang1,3,4, J.Y. Yang11, L.L. Yang11, M.J. Yang1,3, R.Z.
Yang5, S.B. Yang24, Y.H. Yao21, Z.G. Yao1,3, Y.M. Ye6, L.Q. Yin1,3, N. Yin17,
X.H. You1,3, Z.Y. You1,2,3, Y.H. Yu17, Q. Yuan16, H.D. Zeng16, T.X. Zeng1,3,4,
W. Zeng24, Z.K. Zeng1,2,3⋆, M. Zha1,3, X.X. Zhai1,3, B.B. Zhang15, H.M.
Zhang15, H.Y. Zhang17, J.L. Zhang7, J.W. Zhang21, L. Zhang13, L. Zhang24, L.X.
Zhang10, P.F. Zhang24, P.P. Zhang13, R. Zhang5,16, S.R. Zhang13, S.S.
Zhang1,3, X. Zhang15, X.P. Zhang1,3, Y. Zhang1,3, Y. Zhang1,16, Y.F. Zhang20,
Y.L. Zhang1,3, B. Zhao20, J. Zhao1,3, L. Zhao4,5, L.Z. Zhao13, S.P. Zhao16,17,
F. Zheng8, Y. Zheng20, B. Zhou1,3, H. Zhou19, J.N. Zhou18, P. Zhou15, R.
Zhou21, X.X. Zhou20, C.G. Zhu17, F.R. Zhu20, H. Zhu7, K.J. Zhu1,2,3,4, X.
Zuo1,3,
(The LHAASO Collaboration)  1 Key Laboratory of Particle Astrophyics &
Experimental Physics Division & Computing Center, Institute of High Energy
Physics, Chinese Academy of Sciences, 100049 Beijing, China
2University of Chinese Academy of Sciences, 100049 Beijing, China
3TIANFU Cosmic Ray Research Center, Chengdu, Sichuan, China
4State Key Laboratory of Particle Detection and Electronics, China
5University of Science and Technology of China, 230026 Hefei, Anhui, China
6Department of Engineering Physics, Tsinghua University, 100084 Beijing, China
7National Astronomical Observatories, Chinese Academy of Sciences, 100101
Beijing, China
8National Space Science Center, Chinese Academy of Sciences, 100190 Beijing,
China
9School of Physics, Peking University, 100871 Beijing, China
10Center for Astrophysics, Guangzhou University, 510006 Guangzhou, Guangdong,
China
11School of Physics and Astronomy & School of Physics (Guangzhou), Sun Yat-sen
University, 519082 Zhuhai, Guangdong, China
12School of Physical Science and Technology, Guangxi University, 530004
Nanning, Guangxi, China
13Hebei Normal University, 050024 Shijiazhuang, Hebei, China
14School of Physics and Microelectronics, Zhengzhou University, 450001
Zhengzhou, Henan, China
15School of Astronomy and Space Science, Nanjing University, 210023 Nanjing,
Jiangsu, China
16Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain
Observatory, Chinese Academy of Sciences, 210023 Nanjing, Jiangsu, China
17Institute of Frontier and Interdisciplinary Science, Shandong University,
266237 Qingdao, Shandong, China
18Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical
Observatory, Chinese Academy of Sciences, 200030 Shanghai, China
19Tsung-Dao Lee Institute & School of Physics and Astronomy, Shanghai Jiao
Tong University, 200240 Shanghai, China
20School of Physical Science and Technology & School of Information Science
and Technology, Southwest Jiaotong University, 610031 Chengdu, Sichuan, China
21College of Physics, Sichuan University, 610065 Chengdu, Sichuan, China
22Key Laboratory of Cosmic Rays (Tibet University), Ministry of Education,
850000 Lhasa, Tibet, China
23School of Physics and Technology, Wuhan University, 430072 Wuhan, Hubei,
China
24School of Physics and Astronomy, Yunnan University, 650091 Kunming, Yunnan,
China
25Yunnan Observatories, Chinese Academy of Sciences, 650216 Kunming, Yunnan,
China
26College of Sciences, Northeastern University, 110819 Shenyang, Liaoning,
China
27Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, 2 Dublin,
Ireland
28Max-Planck-Institut for Nuclear Physics, P.O. Box 103980, 69029 Heidelberg,
Germany
29 Dipartimento di Fisica dell’Università di Napoli “Federico II”, Complesso
Universitario di Monte Sant’Angelo, via Cinthia, 80126 Napoli, Italy.
30Institute for Nuclear Research of Russian Academy of Sciences, 117312
Moscow, Russia
31Moscow Institute of Physics and Technology, 141700 Moscow, Russia
32Département de Physique Nucléaire et Corpusculaire, Faculté de Sciences,
Université de Genève, 24 Quai Ernest Ansermet, 1211 Geneva, Switzerland
33Department of Physics, Faculty of Science, Mahidol University, 10400
Bangkok, Thailand
caozh@ihep.ac.cn;zengzk<EMAIL_ADDRESS>
wangyanjin<EMAIL_ADDRESS>
###### Abstract

The Wide Field-of-View Cherenkov Telescope Array (WFCTA) and the Water
Cherenkov Detector Arrays (WCDA) of LHAASO are designed to work in combination
for measuring the energy spectra of various cosmic ray species over a very
wide energy range from a few TeV to 10 PeV. The energy calibration of WCDA can
be achieved with a proven technique of measuring the westward shift of the
Moon shadow of galactic cosmic rays due to the geomagnetic field. This
deflection angle $\Delta$ is inversely proportional to the energy of the
cosmic rays. The precise measurements of the shifts by WCDA allows us to
calibrate its energy scale for energies as high as 35 TeV. The energy scale
measured by WCDA can be used to cross calibrate the energy reconstructed by
WFCTA, which spans the whole energy range up to 10 PeV. In this work, we will
demonstrate the feasibility of the method using the data collected from April
2019 to January 2020 by the WFCTA array and WCDA-1 detector, the first of the
three water Cherenkov ponds, already commissioned at LHAASO site.

LHAASO, WCDA, WFCTA, Moon Shadow, Energy Scale, composition

††preprint: AIP/123-QED

##  I Introduction

Cosmic ray experiments based on the extensive air shower (EAS) technique
usually feature different types of ground based detectors, such as a
scintillation counter array, water Cherenkov detectors, or imaging air
Cherenkov telescopes, each measuring shower properties in different ways.
Therefore, it is mandatory to establish a way to calibrate the shower energy
measurement between the different detectors, a non-trivial task given that it
has to be done directly using cosmic ray data. Many experiments, such as ARGO-
YBJ ARGO-YBJ , have successfully their energy-scale calibration by using the
deficit of cosmic rays blocked by the Moon, as it moves inside the field of
view (FoV) of the detectors, usually referred to as the Moon shadow of cosmic
rays. The geo-magnetic-field (GMF), deflecting the charged cosmic rays, shift
the Moon shadow on the ground with respect to the Moon real position. The
displacement of the shadow is clearly dependent on the cosmic ray rigidity and
becomes negligible at high energies. However, at energies below 40 TeV the
shadow shift westward is clearly observable.

The Water Cherenkov Detector Array (WCDA) of the Large High Altitude Air
Shower Observatory (LHAASO) has a detection threshold of about 1 TeV for
cosmic rays. WCDA has measured the Moon shadow shifts as a function of
$N_{fit}$. Here, $N_{fit}$ is the number of triggered units with trigger time
within 30 ns from the conical front. Units are the 5 m $\times$ 5 m cells in
which the pond is subdivided. A data set are collected from 2019/05/01 to
2020/01/31 with 731.2 hours of the Moon observation with zenith angle smaller
than $45^{\circ}$. The reconstruction of arriving directions and shower cores
can be seen in the reference WCDA-on-Crab . During the calculation of the
significance, events with cores located both inside and outside the pond are
used. The significance of the Moon shadow for each $N_{fit}$ group is greater
than $10\sigma$, enabling a good measurement of both the location and the
width of the distribution. The deflection angle of the shadow as a function of
$N_{fit}$ is plotted in Fig. 1. The measured deflection can be approximated by
the simple form $\Delta=k\cdot N_{fit}^{\alpha}$, where $k=-39.2\pm 2.9$ and
$\alpha=-0.939\pm 0.028$.

Figure 1: The displacement of the Moon shadow from the real direction of the
Moon as a function $N_{fit}$.

At the same time, the shadow shift can be evaluated using a ray tracing
algorithm with a detailed GMF model GMF-model  and for different cosmic ray
compositions, i.e., a mixture with different ratios of protons and alpha
particles as a function of the their energies. The deflection angle measured
by the Moon shadow shift provides a means to correlate the energy estimator
and the energy of primary particles, thus providing a simple method for
measuring the energy of the showers. In the energy range below 10 TeV,
$N_{fit}$ is selected as the shower energy estimateWCDA-on-Crab . Each group
of events in $N_{fit}$ bins is assigned for the energies corresponding to the
shifts just like what have done in ARGO-YBJ experimentARGO-YBJ .

The Wide FoV Cherenkov Telescope Array (WFCTA) of LHAASO measures shower
characteristics using the Cherenkov light generated by secondary particles
produced in the air during the whole shower development. The faint Cherenkov
light is overwhelmed by moonlight, and can not be detected near the direction
of the Moon, making it impossible to calibrate the energy scale by measuring
the Moon shadow shift. Nevertheless, we can use events that trigger both WCDA
and WFCTA, to bridge the energy scale of WCDA to the WFCTA measurements. Since
the trigger threshold of WFCTA is higher than for WCDA, only events with
energies above 15 TeV can be used for the cross-calibration. At such high
energies, the energy estimator of $N_{fit}$ for WCDA-1 is found saturated
because of the natural limit of the 900 units of the detector. A more adequate
energy estimator than $N_{fit}$ has to be defined first. The real challenge
lies in the fact that the shift of the Moon shadow becomes small at such high
energies and, at the same time, the number of events for this measurement
decreases as energy increases. Therefore, at this early stage of the LHAASO
experiment, the limited statistics is the dominant uncertainty in this
calibration procedure. Clearly, accumulating data over time, the statistical
error will be reduced until it becomes less important than the systematic
uncertainties due to the unknown composition of the cosmic rays and the
hadronic interaction models used to simulate showers. The overall uncertainty
of the energy scale should become lower than 10% after four years of
observation.

In this work, firstly we briefly describe the LHAASO Observatory detectors in
section II. In section III, we discuss how the energy scale for WCDA can be
established using the Moon shadow shift measurement for energies above 6 TeV.
The uncertainties of the energy scale are also discussed in the section.

How the WCDA energy scale is propagated to WFCTA for calibrating the
reconstruction of shower energies is described in section IV and, finally, in
section V we draw some conclusions.

##  II The LHAASO WCDA and WFCTA arrays

The LHAASO Observatory is based on the so-called ’hybrid’ approach for the
measurement of shower characteristics, consisting in the simultaneous
detection of atmospheric showers with different types of detectors. The
observatory is built around the three ponds of water Cherenkov Detector Array
(WCDA), featuring 3120 gapless detecting units to instrument an area of 78,000
m2. Near WCDA 18 wide field of view Cherenkov telescopes (WFCTA) are
installed. They survey the sky above the whole array with a coverage of 4608
square degrees. This core of the array is surrounded by 5195 scintillation
counters (ED) and 1188 muon detectors (MD), which constitute an array covering
an area of 1 km2 (KM2A). The construction of LHAASO has nearly been completed.

Figure 2: The LHAASO layout. The three ponds of WCDA are represented by the
cyan rectangles at the center of the site, each divided into square cells of 5
$\times$ 5 m2 (as can be seen in the zoomed view at the bottom right). The 18
WFCTA telescopes, positioned near the WCDA, can be seen in the zoomed view at
the top right. The remaining KM2A extends over an area of about 1 km2,
instrumented with the electromagnetic detector (ED) array of scintillation
counters (small red dots) and the muon detector array (big blue dots).

As shown in Fig. 2, WCDA is composed of two ponds with an area of 150 m
$\times$ 150 m each and a third larger one with an area of 300 m $\times$ 110
m. The smaller pond in the South-West direction, named WCDA-1, has started
science operations since April 2019. It has 900 units, or cells, of 25 m2,
each equipped with a large (8”) PMT used also for timing and a small (1.5”)
PMT at the center of each unit at 4 m of depth from the water surface. The use
of two PMTs watching upwards allows us to cover a wide dynamic range spanning
from 10 to 30,000 photo-electrons, which are produced in water by the shower
charged secondary particles. To suppress the light cross-talking effect and
improve the timing resolution, black plastic curtains delimit the units. The
front-end electronics (FEE) of the large PMTs is designed to achieve a time
resolution of 0.5 ns, which together with the large dynamic range, enables the
measurement of the particle density distribution in the shower cores without
any saturation even for energetic showers up to 10 PeV. This allows to achieve
determination of the core location with a precision better than 3 m over a
wide energy range. WCDA-1 can measure shower directions with a resolution
better than 0.2∘ above 10 TeV and 1.0∘ above 600 GeV WCDA-on-Crab . The water
transmission can be monitored by the ”muon” peak observed by each unit WCDA-
on-Crab .

Six WFCTA Cherenkov telescopes in the southwest corner of WCDA-1, and two more
in the southeast corner were put into science operation during winter 2020.
The centers of the six telescopes are located 100 m south and 38.5 m west of
the wcda-1 center. As shown in Fig. 2, the final layout will include 18
telescopes. The telescopes have mirrors with an area of about 5 m2, composed
of aluminized spherical mirror facets. Each telescope can be tilted up and
down in elevation from 0∘ to 90∘. A camera with 32 $\times$ 32 square pixels,
realizing a FoV of 16${}^{\circ}\times$16∘, is in the focal plane of the
telescope, at 2870 mm away from the mirror center. Each pixel is realised by a
silicon photomultiplier (SiPM) coupled to a square-surface light-funnel of 1.5
cm $\times$ 1.5 cm, corresponding to a FoV of 0.5${}^{\circ}\times$ 0.5∘. Each
SiPM is composed by an array of 0.36 million avalanche photo diodes (APDs)
with a micro-cell size of 25 $\mu$m, covering a dynamic range from 10 to
40,000 photo-electrons with a measured non-linearity less than 5% SiPM . In
front of the photo-sensors and light-funnels, a window is installed. The
window is coated with a wide-band filter to suppress the light above 550 nm,
where the night sky background light (NSB) dominates with respect to the
Cherenkov signal. This is relevant since the photo-detection of SiPMs is
considerably larger than PMT one at these wavelengths. The performances of
SiPM are affected by the temperature and night sky background SiPM-performance
. A calibrated led is installed in front of the camera to calibrate the
performance of SiPMs in real time.

Each group of 16 pixels is clustered together in a FEE board connected with
the readout system. Each pixel has a high-gain and a low-gain channel, each
equipped with a 50 MHz 12-bit flash ADC to digitize the waveform. These two
channels allow to cover the whole dynamic range of the SiPMs. A trigger
signal, $T_{0}$, is generated by the high-gain channel of each pixel and
transmitted to the trigger logic that collects the signals from all pixels. A
pattern recognition algorithm is used to decide whether or not a shower has
been observed. It generates a signal $T_{1}$, which starts the readout of the
digital waveform from every pixel. The total charge and the average time of
the waveform are transmitted to the LHAASO data center with an absolute time
stamp. Further off-line triggers, in particular the inter-telescope trigger,
and noise pixel rejection are carried out on the CPUs of the data center. A
typical common trigger event with WCDA-1 and two telescopes is shown in Fig.
3.

Figure 3: A typical event that triggered both WCDA and WFCTA. The two tilted
panels on the left-hand side show the images recorded by two of the WFCTA
telescopes, both in vertical view, but with different orientations of
$15^{\circ}$ and $45^{\circ}$, respectively, from North towards West. The two
panels on the right, instead, show the maps of hits in WCDA units, where the
color scale indicates the number of recorded photo-electrons ($N_{pe}$). The
top right panel shows the map for the 8” PMTs, which exhibit a clear
saturation in the central region of the shower. The lower right panel shows
the map for the same event as recorded by the 1.5” PMTs. The corresponding
lateral distributions of $N_{pe}$ are shown in the lower panel in the middle
column. Here, the red dots represent the 8” PMT and the blue ones the 1.5”
PMTs, which are also scaled (green dots) and overlaid on the previous one to
reconstruct the profile over the full range. The upper panel shows the arrival
time (in ns) of shower secondary particles as recorded by the WCDA-1.

##  III Determination of Shower Energy Scale with WCDA-1

The energy region where WCDA-1 and WFCTA can observe the same showers, at the
maximum energy range for WCDA due to the higher energy threshold of WFCTA. In
the energy range, as an energy estimator $N_{fit}$, the saturation appears, so
a better energy estimator than $N_{fit}$ is needed for the energy scale
determination. The energy scale of WCDA-1 can be obtained as following setps:

1: Selecting total number of detected photo-electron by PMTs, $N_{pe}$, as the
energy estimator. The Moon shadow shift westward is measured as a function of
$N_{pe}$, i.e., $\Delta$ vs. $N_{pe}$.

2: Calculating the Moon shadow shifts by tracing the cosmic rays with certain
compositions through GMF as functions of CR promary energy, i.e. $\Delta$ vs.
energy.

3: Investigating the composition of cosmic ray events that trigger WCDA-1 in
the relevant ranges of $N_{pe}$ using the air shower and detector response
simulation. The composition is used as an input to establish the adequate
relation between the measured $\Delta$ and the primary energy.

4: Investigating the effect due to the shower energy resolution and the power-
law like CR spectrum to the match between $N_{pe}$ and the primary energy
according to the corresponding Moon shadow shift. Establishing the energy
scale, namely $N_{pe}$ versus $E_{median}$ in each group of $N_{pe}$.

###  III.1 $N_{pe}$ as the energy estimator for high energy showers

WCDA measures Cherenkov photons generated by the secondary charged particles
of showers, mainly $e^{\pm}$ and $\gamma$-rays, generated in air showers. The
summed number of photo-electrons recorded by the PMTs is dependent not only on
the number of secondaries, but also on their energies. Because the energy of
the secondary particle in the air shower is proportional to the number of
electrons in its induced shower in water, the total number of photo-electrons
measured by units with trigger time within 30 ns from conical plane, denoted
as $N_{pe}$, is a primary energy estimator. The simple form relating the
primary energy and the estimator $E=bN_{pe}^{\beta}$ can be used and the
parameters, $b$ and $\beta$, can be determined by fitting the experimentally
measured Moon shadow shift.

In Fig. 4 the distribution of $N_{pe}$ is shown for the events used in the
Moon shadow analysis, for different values of the minimal number of detector
units used in the reconstruction, $N_{hit}$, i.e., greater than 100, 200, or
300. Showers reconstructed in WCDA-1 using more than 200 units (hits) have a
rather good directional resolution $0.39^{\circ}$ WCDA-on-Crab . From Fig. 4
we can infer that above a certain energy, corresponding to $N_{pe}>50,000$,
the showers are detected by WCDA-1 with rather high efficiency. As a matter of
fact, the distribution of $N_{pe}$ has a clear power law behavior between
30,000 to 106 with a power index of -2.6, which is the same as the spectral
index of the proton spectrum measured around 13 TeV by DAMPE DAMPE-proton .
This indicates that the efficiency does not change in this energy range, and
$N_{pe}$ can be used as an energy estimator.

Figure 4: The distribution of the total number of photo-electrons, $N_{pe}$,
for shower events coming from a region around the Moon nominal position. The
histograms in red, green and blue correspond to events with number of hit
cells, $N_{hit}$, greater than 100, 200 and 300, respectively. The black solid
line indicates a power law $\propto$ $N_{pe}^{-2.6}$, which fits the
histograms for $N_{pe}>50,000$. The two vertical magenta lines indicate the
range used for the energy calibration using the Moon shadow.

###  III.2 Measurement of Moon shadow shifts for high energy showers

The data set in this analysis was collected from 2019/05/01 to 2020/01/31. The
total effective observation time to the Moon with zenith angles smaller than
$45^{\circ}$ is 731.2 hours. To calculate the deflection angle as a function
of the energy, events are grouped into five groups according to $N_{pe}$ as
shown in Table 1. During the significance analysis, only events with zenith
angles $<$ $45^{\circ}$ and $N_{hit}$ $>$ 200 are used.

For the data set in each group, the sky map in celestial coordinates near the
moon region is divided into a grid of $0.02^{\circ}\times 0.02^{\circ}$ bins
and filled with detected events according to their reconstructed arrival
directions. The number of events in each grid denoted as $n$. The number of
background events in each grid is estimated by the equal zenith method
Amenomori , denoted as $b$. The deficit significance in each grid is estimated
by Li & Ma formula lima . The maximum deficit significance around the Moon
region in each $N_{pe}$ group are also shown in Table 1, where the
corresponding Moon shadow westward shifts, uncertainties and are also listed.

For the $N_{pe}$ group with $N_{pe}>60,000$, a significance as high as 10.9
$\sigma$ can be achieved. In Fig. 5, the significance of the deficit is
plotted as a function of the arrival direction in right ascension (RA) and
declination (Dec). The significance map can be fitted by a bi-dimensional
Gaussian function to determine the location of the shadow center. Its shift
with respect to the nominal Moon position is $0.02^{\circ}\pm 0.03^{\circ}$ in
the Dec direction , while is quite small in the RA direction, i.e, $0.005\pm
0.03$. The statistical uncertainty is the dominant contribution given the
limited statistics. The standard deviation $\sigma$ of the Gaussian function
is found to be $0.33^{\circ}\pm 0.05^{\circ}$, which is the result of the
combination of the size of the Moon (angular extension 0.25∘) and the angular
resolution of the detector WCDA-on-Crab . The deficit number of events
contained in the angular region $\sigma_{DEC}^{2}+\sigma_{RA}^{2}<2\sigma^{2}$
is about 63% of the total deficit number, a value that is consistent with the
simulation. For lower energies, the shadow shifts towards the West and its
size increases as well, indicating a worsening of the angular resolution.

Figure 5: The significance map of the Moon shadow for shower events detected
by WCDA-1 with a $N_{pe}$ greater than 60,000. The deflection of the CRs by
the GMF, at such high energies, is seen to be quite small.

###  III.3 Calculation of Moon shadow shifts by raytracing in GMF

The expected Moon shadow shift westward has been calculated by using a ray-
tracing simulation which propagates protons and Helium nuclei coming from the
Moon direction through GMF. The Interplanetary Magnetic Field, due to the
solar wind, is by far less intense and can be neglected. The Tsyganenko-IGRF
model GMF-model  has been used to describe the magnetic field in the Earth-
Moon system. We find that the displacement obtained applying this model to the
propagation of protons and Helium nuclei can be represented above 5 TeV by the
simple dipole approximation

$\Delta=1.59^{\circ}/R(TV)$ (1)

where R(TV) is the particle rigidity E(TeV)/Z. Thus the expected shift for
Helium nuclei is a factor of 2 greater than the shift of protons of the same
energy. However, for a given energy, the shower size $N_{e}$ of Helium nuclei
is less than the size of proton-induced shower. According to the simulation,
the energy of an Helium nuclei generating at the LHAASO altitude a shower size
equal to the size of a proton of energy $E_{P}$ is $E_{He}\approx 1.9E_{P}$.
The rigidity of the Helium nuclei with the same shower size $N_{e}$ of proton
is then $R\approx 1.9E_{P}/2$ , very close to the proton rigidity. This result
has been already obtained in the Moon shadow analysis performed by ARGO-YBJ
ARGO-YBJ . Thus, at a given selected interval of $N_{pe}$, corresponding to a
shower size $N_{e}$ interval, the displacement of Helium-induced showers is
practically equal to that of showers generated by proton primaries. Moreover,
fixing a $N_{pe}$ interval we select a primary cosmic ray beam with a fraction
Helium nuclei : protons is 2. after triggering assuming the same flux for both
components as measured by the CREAM experiment in this energy range CREAM-
proton-helium . The mean energy of this beam is $E_{u}\approx 1.3E_{P}$. Since
the westward shift of these particles is equal to the displacement of protons
of energy $E_{P}$, that is $\Delta$ = 1.59/$E_{P}$ = k/$E_{u}$, we find the
relation $\Delta\approx 2.1/E_{u}$. This result is fully confirmed by a
detailed simulation which provides the amount of the shift $\Delta=(2.1\pm
0.3)/E$ where E is the median energy of the cosmic beam selected by fixing
$N_{pe}$. This relation is shown in Fig. 6 where the deflection angles of
proton and Helium nuclei pure beams are also reported for comparison. Using
the results of Fig.6 we can attribute the cosmic ray median energy to each
$N_{pe}$ selected interval reported in Table 1. Therefore, the energy scale
can be fixed in the $N_{pe}$ range 6,000- 60,000 as shown in Fig. 7 . In the
energy range from 6.6 TeV to 35 TeV $N_{pe}$ can be used as energy proxy
according to the relation $E[GeV]=bN_{pe}^{\beta}$, $\beta=0.95\pm 0.17$ and
$b=1.33^{+5.26}_{-1.06}$.

.

Figure 6: The deflection angle of westward Moon shadow shift is graphed versus the inverse of the CR particle energy for pure protons (red line) and pure Helium nuclei (black line) according to the ray-tracing simulations. The blue squares, which are the measured mean shift of the Moon shadow with its statistical error in the corresponding $N_{pe}$ bins, lie along the pink line representing the mixture of protons and Helium nuclei with ratio of 2.0:1. The uncertainty in the composition is represented by the shaded area around the pink line, which is so small that it is barely visible. Table 1: Moon shadow shifts in RA, the significance of the Moon shadow, and the reconstructed median energy with its error, grouped bins of $N_{pe}$. Range of Npe | Shift of the | Significance | Median E
---|---|---|---
|  Moon shadow (∘) |  $(\sigma$) | (TeV)
6,000-10,000 |  -0.32$\pm 0.04$ | 18.2 |  6.6$\pm 0.8$
10,000-15,000 |  -0.25$\pm 0.04$ | 14.0 |  8.4$\pm 1.3$
15,000-20,000 |  -0.15$\pm 0.04$ | 11.6 |  14.0$\pm 3.7$
20,000-30,000 |  -0.11 $\pm 0.03$ | 11.9 |  19.1$\pm 5.2$
30,000-60,000 |  -0.06 $\pm 0.03$ | 10.8 |  35$\pm 17.5$
$>$60,000 |  -0.005$\pm 0.03$ | 10.9 |  $>$50.0
Figure 7: The average shower energy measured using the Moon shadow shift
versus $N_{pe}$, the total number of photo-electrons detected by WCDA-1
detector.

###  III.4 Uncertainties of Energy Scale

The uncertainty of the energy scale is due mainly to statistics, which are
caused by the uncertainty of the measurements of Moon shadow shifts. The
uncertainty of the Moon shadow shift can be transferred to the energy scale by
using the fact that $d\Delta/\Delta=dE/E$. From Table 1 it can be seen that
the average position of the shadow has an error of 0.04∘ in the lower energy
bins, to become 0.03∘ at higher energies. These errors result in a rather
large uncertainty in the energy scale reflected by the big error bars in Fig.
7, i.e., from 12% at 6.6 TeV to 50% at 35 TeV. These uncertainties are
expected to drop to 3% and 12%, respectively, after four years of operations
of the full WCDA detector.

Two systematic uncertainties may affect this analysis, one of them is the
uncertainty of the ratio of the proton and Helium, whcich is found fluctating
with an aplitude of 10% over the energy range from 5 to 50 TeV according to
the simulation. The other one is caused by the energy and angular resolution
of the detector about 4%.

##  IV Energy Scale of Shower Measurements for WFCTA

In the previous section before the energy scale for WCDA-1 in the energy range
from 6 TeV to 40 TeV has been determined. The correspondence between energy
and its estimator $N_{pe}$ is modelled as $E[GeV]$ = $1.33~{}N_{pe}^{0.95}$.

The methods presented here can serve as a calibration of the energy estimation
only in the specific energy range used in the data analysis, while at higher
energies, the energy scale has to rely solely on simulations of the air
showers and detector response.

Figure 8: (Upper) The distribution of zenith angles of all events used in the
WCDA-1 Moon shadow analysis. The shaded areas represent the zenith range not
covered by WFCTA telescopes. (Lower) The significance map of the Moon shadows,
for events with $N_{pe}$ between 20,000 and 60,0000 and in the zenith angle
range covered by the WFCTA telescopes, i.e., 22${}^{\circ}\leq\Theta_{Z}\leq$
38∘

To detect the faint Cherenkov light, WFCTA telescopes have to be operated on
dark nights or at most with partial moonlight, thus avoiding having the Moon
in the FoV and making it impossible to measure Moon shadow shifts for
determining the energy scale. At the same time, the energy scale of WCDA-1
obtained with the Moon shadow shifts can be used to provide a reference for
the energy calibration of WFCTA, using events triggered by both. This cannot
be done on an event by event basis, but instead with a set of commonly
triggered events, which need therefore to have similar features with the set
used for WCDA-1 analysis.

The WFCTA telescopes are pointed at $30^{\circ}$ in zenith, which results in a
zenith angle coverage from $22^{\circ}$ to $38^{\circ}$, taking into account
the 16∘ FoV of the telescopes. While the events used to determine the WCDA-1
energy scale come from a region defined by the zenith angle range, i.e.,
mainly between $10^{\circ}$ to $45^{\circ}$ in section 4. Fig. 8 shows the
zenith range covered by WCDA-1 and in white the zenith range covered by WFCTA.
Fig. 9 (upper) shows the distribution of $N_{pe}$ for the events used in the
Moon shadow analysis of WCDA (in blue) and the commonly triggered events (in
red) obtained discarding the events falling in the shaded area as of Fig. 8.
Looking at the distribution in Fig. 9, it can be seen that WCFTA triggers for
100% of WCDA-1 events only for $N_{pe}$ above 60,000, and that the efficiency
drops below 75% for $N_{pe}$ below 20,000. This is due to the fact that the
trigger efficiency of WCFTA telescopes decreases with increasing distance from
the shower core and also for decreasing energy. On the contrary, the Moon
shadow shift is better reconstructed for lower energies ($N_{pe}<$60,000) and
when the core is well inside the ponds, i.e., further away from the
telescopes.

Therefore, the common data set has been selected choosing the best compromise
using the following criteria: $N_{hit}>200$, $20,000<N_{pe}<60,000$ and
$22^{\circ}<\Theta_{z}<38^{\circ}$. Fig 9 (lower) shows the distribution of
zenith angles of commonly triggered events with the center of gravity if the
images within 5∘ from the camera center.

Using these selection criteria leads to a new significance map, shown in Fig.
8, with the same analysis method as used for WCDA. The value obtained for the
shift in this case is $\Delta=(0.1\pm 0.03)^{\circ}$, which translates to an
energy 21$\pm 6.5$ TeV, using the formula $\Delta=-2.1/E$[TeV] as before. The
uncertainty quoted is purely statistical and results in a relative error of
the calibration of about 30%. Clearly this result can be significantly
improved over more years at collecting statistics.

The density of Cherenkov photons produced during the shower development is
proportional to the primary particle energy. Therefore a good estimator of the
energy for WFCTA data is the sum of the number of photo electrons in the
camera image, $\sum$, once corrected for geometric effects. The quantity
$\sum$ is usually called the size of the image, in the traditional shower
energy reconstruction scheme ARGO-YBJ  for Cherenkov telescopes. The
correction on $\sum$ includes its dependence on the viewing angle $\beta$, the
space angle between the shower direction and the main axis of the telescope,
and the perpendicular distance, $R_{p}$, from the shower axis to the telescope
in shower-detector-plane (SDP). The $R_{P}$ correction takes into account the
fact that photon density decreases with increasing distance from the shower
core. The viewing angle $\beta$ correction is due to the weakening of the
shower image on the camera for off-axis showers. This is a combination of the
shadow of the container in which the telescope is installed and also the
reduced effective area of the mirror for off-axis showers. The correction has
been calculated with a detailed simulation of the WFCTA telescopes response
YLQ-composition-CPC  omitted here for the sake of brevity.

The shower energy can be estimated with corrected $\sum$ with a resolution of
20%, for energies below 100 TeV, and 15% above 100 TeV with a systematic shift
less than 5% YLQ-composition-CPC .

In order to correctly reconstruct the energy, a proper set of data has been
selected. For WCDA events, apart from the standard cut $20,000<N_{pe}<60,000$,
a requirement is added to have the shower core to be well inside the pond,
i.e., $|x_{core}|<$ 55 m and $|y_{core}|<$ 55 m . In addition, the WFCTA
events are required to be not too far, i.e., $R_{p}<$ 100 m, to have at least
2 triggered pixels, and to have the center of gravity of the images within
$5^{\circ}$ from the camera center. The distribution of energies, $E_{\tiny
WFCTA}$ reconstructed by $\sum$ for the sample selected with these cuts is
shown in Fig. 10. The mean value ( 21.9 TeV )of the distribution is shown by
the vertical red line. The WCDA-derived energies of commonly triggered events
with the same cuts can be reconstructed by the formula $E[GeV]$ =
$1.33~{}N_{pe}^{0.95}$, which is shown in Fig. 10 by the black line. The
median energy (23.6 TeV) of the distributions is shown by the black vertical
line.

For comparison, the energy scale with its statistical uncertainty ( 21$\pm$6.5
TeV) determined from the Moon shadow shift are also reported in Fig. 10. The
agreement between the mean of the energy reconstructed by the WFCTA telescopes
and the energy scale of WCDA-1 is evident and differences are of the order of
4%, well below the 30% uncertainty of the Moon shadow energy determined which
is dominated by the statistical error.

In order to evaluate the bias introduced by the requirement to have the shower
core well within WCDA, the energy scale has been estimated by the Moon shift
method adding the conditions $|x_{core}|<$ 55 m and $|y_{core}|<$ 55 m to the
previous ones,i.e., $N_{hit}>200$, $20,000<N_{pe}<60,000$ and
$22^{\circ}<\Theta_{z}<38^{\circ}$. The resulting deflection angle is
$\Delta=(0.13\pm 0.05)^{\circ}$, which translates into an energy of
16.2$\pm$6.2 TeV, according to the formula $\Delta=-2.1/E$(TeV), this is shown
as a blue square with its error, corresponding to a relative error of 38% in
Fig. 10. The change of the energy scale introduced by the further require of
having the shower core well within WCDA is about 30%.

Figure 9: The distributions of $N_{pe}$ (upper) and zenith angles (lower) for
the events commonly triggered by both WCDA-1 and WFCTA. The distribution of
$N_{pe}$ for WCDA-1 only (blue) is compared with that of common events used in
the Moon shadow measurements described above (red). In the selected range of
$N_{pe}$, indicated by the vertical lines, the selection efficiency is high
enough to have almost the same ratio, proton : Helium nuclei = 2 : 1. Figure
10: The distributions of reconstructed energies. The black line indicates the
distribution of energies reconstructed by WCDA-1 for the commonly triggered
events, using the estimator $E[GeV]=1.33\cdot N_{pe}^{0.95}$. The vertical
black line indicates the median energy (23.6 TeV) of the distribution. The red
line indicates energies reconstructed by WFCTA telescopes assuming a primary
composition ratio of protons : Helium nuclei = 1 : 1. The vertical red line
indicates the mean energy (21.9 TeV) of the distribution. The energy scale
with it uncertainty (21$\pm$6.5 TeV) measured using the Moon shadow shift with
available data (red square). The blue square indicates the energy scale
estimated only for events with a core located inside WCDA-1.

To extrapolate the WFCTA energy scale as determined in the overlapping region
with WCDA, the correspondence between the energy estimator of WFCTA and the
primary energy has been checked using simulations. Fig. 11 shows the
distribution of $E_{\tiny WFCTA}$ reconstructed by WFCTA as a function of
input primary energy $E_{0}$, spanning a range from 30 TeV to 10 PeV. In
addition, the energy scale determined from Moon shadow shifts in WCDA is
graphed as a black square together with the error bar at $21.0\pm 6.5$ TeV.
The uncertainty here is dominated by the statistical errors, as the systematic
error coming from composition uncertainty, about 3%, contribute marginally.

Figure 11: Distributions of $E$ reconstructed by WFCTA and primary energy
$E_{0}$ for primary composition ratio protons : Helium nuclei = 1 : 1 in the
energy range form 10 TeV to 10 PeV. The black square indicates the energy
scale (21 TeV) obtained from the Moon shadow shift and the mean energy as
reconstructed by WFCTA telescopes. The Color scale indicates the relative
number of events in each bin.

##  V Conclusion

In this work, we have shown how the shift of Moon shadow in cosmic rays, due
to the geomagnetic field, can be used to establish the energy scale of the
WCDA detector in the range from 6 TeV to 35 TeV, using the estimator
$E[GeV]=1.33\cdot N_{pe}^{0.95}$, based on the total number of photo-electrons
$N_{pe}$ measured in the pond. The uncertainty of the energy scale varies from
$\pm 12\%$ at 6.6 TeV to $\pm 50\%$ at 35 TeV and is dominated by statistical
errors. The systematic error coming from the assumption of the composition of
the primary ratio protons : Helium nuclei with respect with the direct
measurement by space born experiments CREAM-proton-helium  has been shown to
be about 3%.

Given the impossibility to measure the Moon shadow directly with WFCTA, a set
of events commonly triggered with WCDA-1 have been used to correlate the
energy scale of WCDA to the measurement using WFCTA telescopes. Unfortunately,
the two types of detector have different energy thresholds, and this overlap
occurs at upper limited energy range of WCDA, where the Moon shadow shift
method is less reliable. This results in a rather large uncertainty of about
$30\%$ on the ’common’ energy scale. As a matter of fact, this uncertainty is
largely dominated by the low statistics of event, given the small contribution
of the systematic error from the lack of information on the composition in
this energy. Other systematic uncertainties related to the hadronic
interaction model used in the simulation, at this energies, should be smaller
than the previous one.

In few years LHAASO will accumulate huge statistics, which will allow us to
reduce the uncertainty of the energy scale below 12%. For showers with energy
above 35 TeV, the energy reconstruction totally relies on the simulation of
the WCFTA, which indicates that the energy estimator $\sum$ will provide a
good energy resolution.

##  VI Acknowledgments

The authors would like to thank all staff members who work at the LHAASO site
all year around to keep the system running efficiently and smoothly, even in
the demanding conditions at a mean altitude of 4400 meters above sea level. We
are grateful to the Chengdu Management Committee of Tianfu New Area for the
constant financial support to the research with LHAASO data. This research
work is also supported by the National Key R&D program of China, with the
grant 2018YFA0404201, 2018YFA0404202 and 2018YFA0404203, the National Natural
Science Foundation of China, with NSFC grants 11635011, 11761141001, 11905240,
11503021, 11205126, 11947404, 11675187, U1831208, Schools of Science and
Technology Plan from SiChuan Province grant No.20SYSX0294, and Thailand
Science Research and Innnovation grant RTA6280002.

## References

  * (1) B. Bartoli et al. (ARGO-YBJ Collaboration), Phys. Rev. D 84 022003 (2011)
  * (2) F. Aharonianet al. arXiv:2101.03508v1 [astro-ph.IM]
  * (3) N.A. Tsyganenko,J. Geophys. Res.,100, 5599 (1995).
  * (4) Z. Cao (for LHAASO Coll.), Chin.Phys. C34 249-252 (2010);
  * (5) H. He (for LHAASO Coll.), Radiation Detection Technology and Methods 2-7, (2018)
  * (6) B. Bi et al.(for LHAASO Coll.), Nucl. Instrum. Methods Phys. Res. A899 (2018) 94-100.
  * (7) F. Aharonianet al. arXiv:2012.14622v1 [physics.ins-det]
  * (8) S. Zhang et al. (for LHAASO Coll.), Nucl.Instrum.Meth. A629 (2011) 57-65.
  * (9) M. Amenomori et al., Astrophys. J, 633, 1005-1012 (2005)
  * (10) Li, T. P., Ma, Y. Q. 1983, ApJ, 272, 317
  * (11) Q.An,et al.,Science Advances 5 (9) (2019) eaax3793
  * (12) Y.S.Yoon, et al., Astrophysical J., 839:5 (2017)
  * (13) A.M. Hillas, Annu.Rev.Astron.Astrophys. 22(1)425-444 (1984)
  * (14) Y.L.Qiao, et al.,Chinese physics c, 43(7) 075001 (2019)
  * (15) T.Pierog et al. arXiv:1306.0121[hep-ph] (2013)
  * (16) S.S. Ostapchenko, Phys.Rev. D83 014018(2011)