00footnotetext: Received 8 Jan. 2021

# Performance of LHAASO-WCDA and Observation of Crab Nebula as a Standard
Candle

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,⋆
caozh@ihep.ac.cn;linsj6@mail.sysu.edu.cn;gaobo@ihep.ac.cn;wuhr@ihep.ac.cn;zham@ihep.ac.cn;gmxiang<EMAIL_ADDRESS>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.
Wang1,3 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
26Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, 2 Dublin,
Ireland
27Max-Planck-Institut for Nuclear Physics, P.O. Box 103980, 69029 Heidelberg,
Germany
28 Dipartimento di Fisica dell’Università di Napoli “Federico II”, Complesso
Universitario di Monte Sant’Angelo, via Cinthia, 80126 Napoli, Italy.
29Institute for Nuclear Research of Russian Academy of Sciences, 117312
Moscow, Russia
30Moscow Institute of Physics and Technology, 141700 Moscow, Russia
31Département de Physique Nucléaire et Corpusculaire, Faculté de Sciences,
Université de Genève, 24 Quai Ernest Ansermet, 1211 Geneva, Switzerland
32Department of Physics, Faculty of Science, Mahidol University, 10400
Bangkok, Thailand

###### Abstract

The first Water Cherenkov detector of the LHAASO experiment (WCDA-1) has been
operating since April, 2019. The first 10 months of data have been analyzed to
test its performance by observing the Crab Nebula as a standard candle. The
WCDA-1 achieves the sensitivity of 65 mCU per year with a statistical
threshold of 5 $\sigma$. In order to do so, 97.7% cosmic ray background
rejection rate around 1 TeV and 99.8% around 6 TeV with an approximately
photon acceptance about 50% by using the $compactness$ of the shower
footprints to be greater than 10 as the discriminator between gamma induced
showers and the cosmic ray backgrounds. The angular resolution is measured
using the Crab Nebula as a point source about 0.45∘ at 1 TeV and better than
0.2∘ above 6 TeV with the pointing accuracy better than 0.05∘. They are all
matching the design specifications. The energy resolution is found 33% for
gamma rays around 6 TeV. The spectral energy distribution of the Crab Nebula
in the range from 500 GeV and 15.8 TeV is measured and in agreement with
results of other TeV gamma ray observatories.

###### keywords:

LHAASO-WCDA, Crab Nebula, angular resolution, spectral energy distribution

###### pacs:

P

ACS codes ( 95.85.Pw, 96.50.sd, 98.70.Sa)

00footnotetext: We are grateful to all source of financial support for the
corresponding research work. They are listed in the acknowledgement section.

##  1 Introduction

The Large High Altitude Air Shower Observatory (LHAASO) [2] is a complex of
extensive air shower (EAS) detector arrays, located at Mt. Haizi (29∘21’27.6”
N, 100∘08’19.6” E) and the altitude of 4410 m above sea level (a.s.l.), in an
area at the edge of the Qingzang plateau near Daocheng, Sichuan Province,
China. It consists of an extensive air shower detector array covering an area
of 1.3 km2 (KM2A) with 5195 scintillator counters (ED, 1 m2 active area) and
1188 muon detectors (MD, water Cherenkov detector with an area of 36 m2 buried
under 2.5 meters of dirt). In the center of the array are the Water Cherenkov
Detector Array (WCDA) covering 78,000 m2 and Wide Field-of-view air
Cherenkov/fluorescence Telescope Array (WFCTA) of 18 telescopes. LHAASO is
designed for detection of air showers induced by all kinds of cosmic ray
particles with energy ranging from a few tens of GeV to a few EeV, and for
very effectively identifying gamma rays out of the charged cosmic ray
background. LHAASO is, therefore, a multi-purpose facility for very high
energy gamma ray astronomical observation, cosmic ray energy spectrum
measurements for individual species and other wide-ranging topics associated
with the cosmic rays or air shower phenomena. For the gamma ray source survey,
the designed sensitivity is about 1% Crab Unit (CU). Its wide FoV allows
LHAASO to survey 1/7 of the northern sky at any moment for gamma ray sources.
With the operation in full duty cycle, LHAASO scans the entire northern
hemisphere every 24 hours. The spectra of all sources in its FoV will be
measured with high precision over a wide energy range from $10^{11}$ eV to
$10^{15}$ eV. These measurements, particularly above 100 TeV, will offer a
great opportunity for finding cosmic baryon PeVatrons and therefore
identifying cosmic ray origins among the gamma ray sources. The
Cherenkov/fluorescence telescopes of WFCTA are used to record the longitudinal
development of the air showers that trigger the arrays on the ground. The
combined detection of showers using all detectors in LHAASO enables the
identification of the cosmic ray species, at least for protons, H+He and iron
nuclei with an aperture greater than 4000 m${}^{2}\cdot$sr. The knees of the
spectra of cosmic rays are expected to be measured at energies above 100 TeV
for individual mass groups.

The Water Cherenkov Detector Array (WCDA), mainly used for surveying transient
phenomena and discovering new sources, has been built in 3 phases. The first
one, denoted as WCDA-1, is a square pond of 150 m$\times$150 m consisting of
900 detector units. It has been completed in April 2019 and in operation since
then. In the second phase, one more pond of the same size, referred to as
WCDA-2, has also been in operation since November 2019. The third pond,
WCDA-3, with a size of 300 m$\times$110 m, began operation at beginning of
2021. For gamma-ray induced showers, WCDA-2 and WCDA-3 have a threshold below
100 GeV because they use a 20” PMT in each unit, while 600 GeV is the
threshold for WCDA-1 because an 8” PMT is equipped in each unit. The threshold
could be treated as such an energy above which the detector has a sensitivity
$>1$ CU.

In this paper, the performance of WCDA-1 is thoroughly studied using the data
taken from April 2019 to March 2020. The experiment, including the detector
configuration and the calibration for both timing and charge measurements, is
described in the second section. The details about the observation of Crab
Nebula as a standard candle are described in the third section, including the
shower geometric reconstruction, the cosmic ray background suppression and
angular resolution. The details about air shower and detector response
simulation and the gamma-ray shower energy reconstruction are discussed in the
fourth section. The measurement of the spectral energy distribution of the
Crab Nebula is described along with the systematic uncertainties in the fifth
section. A short summary follows at the end of the paper.

##  2 The WCDA-1 Detector Array

WCDA-1 is composed of 900 detector units optically separated by non-reflecting
black plastic curtains. Each unit is a 5 m $\times$ 5 m cell equipped with two
upward-facing PMTs on the bottom at the center of the unit. The pond is filled
with purified water up to 4 m above the photo-cathodes of the PMTs, with a
total volume of about 100,000 cubic meters. This is the first completed
detector array out of three, as illustrated in the layout in Fig. 1. The inner
surface of the pond is covered with high density polyethylene (HDPE) film with
a thickness of 2 mm in order to guarantee a daily leakage rate less than 0.02%
of the total water capacity.

Figure 1: Schematic of the WCDA layout. The lower two ponds are WCDA-1 and
WCDA-2 from left to right, respectively, and the upper one is WCDA-3. The two
dots in each detector unit indicate the PMTs, while the dashed lines represent
the curtains between units. Two combinations of PMTs are used: 8” and 1.5”
PMTs in WCDA-1 and 20” and 3” PMTs in WCDA-2 and WCDA-3.

To maintain an attenuation length for near-ultra violet light longer than 15
meters, the water needs to be purified continuously. This is done with a
closed recycling loop. Water is sucked out from one side and filled back in on
the opposite side of the pond after purification with bacterium/organic-carbon
killing and filtering. The whole volume is recycled every 15 days.

In Fig. 1, two dots in each unit indicate two PMTs of 8” and 1.5”,
respectively, used in WCDA-1. The combination of the two covers a wide dynamic
range in number of photoelectrons (PE) from 1 to 200,000. This configuration
allows a better measurement of the particle density distribution in the shower
cores without saturation even for very energetic showers up to 10 PeV. It also
allows the shower core location to be reconstructed with a resolution better
than 3 m for 10 TeV and more energetic events. This greatly improves the
shower energy resolution for those events. It is useful for high precision
measurements of spectral energy distribution (SED) of gamma-ray sources in the
overlapping region with LHAASO-KM2A, and It is essential for measuring the
charged cosmic ray spectrum in the combined experiment with LHAASO-WFCTA. An
event detected by both large and small PMTs is shown in Fig. 2. It clearly
shows how the combination of the charge measurements (in units of PE) from the
two PMTs correctly measures the lateral distribution for high energy showers.

Figure 2: A typical high energy shower event detected by WCDA-1. The top left
panel shows the arrival time of the shower front at each cell. The units on
the colour scale are ns. The top right panel shows the number of
photoelectrons, NPE, recorded by the 8” PMT in each cell (bin) in which many
saturated PMTs can be seen. To improve the measurement of the core region, the
two maps are combined, as shown in lower left panel, where the lateral
distribution, measured in equivalent NPE , of the small PMTs (blue points
represent its original measurement and green points represent the charge
equivalent to the charge that would be measured by the 8” PMT) is overlaid to
that of the big PMTs (red points) to construct a distribution over the whole
dynamic range. The lower right panel shows the same event, where the
information of the 2 PMTs is combined.

The 8” PMTs are read out from both the anode and the last dynode as two
independent signals with a gain difference of a factor of about 48. Two 12-bit
ADC channels are used to digitize both signals and achieve a dynamic range
from 1 PE to 5000 PE with a non-linearity less than 1.5% and a charge
resolution of 35% for single PE and 2% for signals greater than 200 PE. The
anode signal is amplified and fed into a discriminator to measure the trigger
time. A timing resolution of 1 ns has been achieved, enough to reconstruct the
shower front conical structure. The 900 units in WCDA-1 are divided into 25
groups of 36 units (72 PMTs), each connected to a single electronic box, which
provides the power and hosts all electronics boards.

The White Rabbit Protocol (WR)[3] is used to distribute the clock and
guarantee the synchronization of each detector unit within 200 ps through a
fibre network over the entire LHAASO array. In the central timing station, a
rubidium oscillator synchronized with GPS/BeiDou Time (BDT) is the source of
the clock, which is distributed to the whole array through the WR network.
Data collected from each PMT are transmitted to the data station of WCDA-1,
through the same fiber network. The same is true for any other LHAASO
detector, where the same network is used for synchronization and data readout.
In WCDA-1, the correct synchronization of the 36 large PMTs, used for timing
measurements within each cluster, is achieved with a dedicated calibration
described below.

A hit is formed with the threshold of 1/3 PE for a 8” PMT and recorded with
the absolute time stamp from the WR system and charges simultaneously measured
by the 8” PMT and the 1.5” PMT. The stream of hits is transmitted to the data
acquisition (DAQ) computers and buffered for 30 seconds before being
discarded. The hits are sorted in time and fed into the triggering and event
building algorithm. This maximises the flexibility by allowing many trigger
algorithms to run in parallel. Built events are stored for further analysis in
specific data streams according to the trigger algorithms.

###  2.1 Timing and Charge Calibration

In the WCDA-1 architecture, 900 units are divided into 25 groups. The 36 8”
PMTs in each group are connected to an electronic box, which is synchronized
with the WR system. In order to calibrate the 36 PMTs in a group, an optical
fiber bundle is illuminated by a LED and used to distribute synchronized light
pulses to every PMT in the group. The LED is driven by a pulse generator
installed in the electronics box. Before routing the fibers to PMTs, the
arrival time of the light pulse generated by the LED through each fibre is
measured [4] in a lab with an accuracy of better than 100 ps. The time
differences between fibers are shown in Fig. 3 (b) as a distribution with a
Root Mean Square (RMS) of 163 ps among 1800 fibers. In order to check the
synchronization between the different groups, another bundle of 36 fibers is
also installed in each group. Twenty fibers of this bundle are used to
illuminate 20 PMTs belonging to the group, and the remaining 16 fibers
illuminate 4 PMTs in each of the 4 adjacent groups as shown in Fig. 3 (a). The
Root Mean Squared (RMS) time offsets measured between adjacent groups is 200
ps. This confirms the synchronization accuracy achieved by using the WR
system.

The calibration using the fiber network has been carried out on site several
times, including the first one during the ‘dry-run’ before water filling and
those in March, April, May, August and October 2019. The distribution of time
offsets measured in the last calibration run is shown in Fig. 3 (c). The RMS
of the offsets is 3.2 ns among the 900 PMTs. The calibration data and the
calibration system have proven to be very stable over time. The time offsets
are very stable with a maximal variation less than 0.2 ns among 900 PMTs over
all the calibrations over eight months.

Figure 3: (a) A sketch map of time cross calibration. (b) The distribution of
the time offset for 1800 fibers used in the calibration system of WCDA-1. (c)
The distribution of the time spread among the 900 units of WCDA-1.

As mentioned above, the 8” PMT is read out from both the anode and the last
dynode to cover a range from 1/3 PE to 5000 PE. The charge calibration of the
anode signal is straightforward by measuring the single PE peak position of
the ADC count distribution for each PMT. The resolution of the single PE
position, with 25 ADC counts as the typical value, is about 32%. The range
covered by the anode high gain channel is up to 130 PE with non-linearity less
than 5%. The details of the calibration are reported elsewhere [5] in which
the measurements of the linearity of every PMT have been summarized. The
calibration of the dynode signals is done by comparing them with the
corresponding anode signals in the overlapping region, namely above 25 PE. In
Fig.4 (a), the dynode signal is plotted versus the anode signal in terms of
ADC counts for a typical PMT. There exists saturation when the anode signal
exceeds 4096. The ratio between the two signals, denoted as $A/D$, could be
extracted by fitting in the region of the anode signals between 500 and 2300
counts. The ratio found, $A/D=53.9\pm 0.1$, is used as a calibration constant
for the PMT. A 2-D map of the 900 $A/D$ ratios of WCDA-1 is shown in the right
panel of Fig. 4. To investigate the calibration stability, the single PE peak
position and $A/D$ of all 900 PMTs have been monitored continuously. The
temporal variations were found to be within 0.1% and 0.05%, respectively,
during April and May 2019.

Figure 4: (a): An example of the amplitude, in ADC counts, of the last dynode
versus the anode signal for an 8” PMT. (b): Uniformity over the WCDA-1, as a
2D map, of the measured A/D ratios for all 900 PMTs.

The calibration of 1.5” PMTs has been done by matching the signals with the 8”
PMTs in the overlapping region, i.e., above 500 PE measured by the 8” PMTs.
Fig. 5(a) shows an illustrative example of how the ratio of the two signals is
extracted with a linear fit in the overlapping region, between 500 and 2300
ADC counts of the 8” PMT, where the non-linearity is below 3%. The ratio is
$8.18\pm 0.04$ for the pair of PMTs as a example. The distribution of the
ratios of all 900 pairs has a width of 17%. The mean value of 9.0 is related
the ratio of the geometrical areas of the photocathodes of the two PMTs and
their gains.

Figure 5: (a): Amplitude in ADC counts of dynode signal of 8” PMT versus anode
signal of 1.5” PMT in a typical PMT unit of the WCDA-1 detector. (b): This
measured A/D ratio for all the 900 PMTs.

The ‘spectral charge distribution’ of cosmic rays measured by each WCDA-1
unit, i.e., $Q{\rm d}N/{\rm d}Q$ where $Q$ is the charge measured by the PMT
in units of PE and $N$ is the trigger frequency of the unit by cosmic rays, is
found to have a clear peak as shown in Fig.6 for a typical unit. This is due
to the diffuse single cosmic ray muons which arrive at each unit with a
distribution that is uniform and nearly constant from a wide range of
directions. Studies show that the peak can be mainly attributed to muons
passing by the PMT with a distance of about 3 meters [6]. The measured charge
at the peak indicates a combination of the PMT sensitivity and the
transmission of water in the unit. This provides an ideal calibrating beam for
WCDA-1 units to measure the uniformity in response of the 900 units and
monitor the variation of the calibration with time, e.g., due to changes in
water transmission. In the left panel of Fig.7, the distribution of the charge
at peak of 900 units is shown in the 2-D map of WCDA-1. The average charge of
900 peaks was 11.6 PE on Nov. 22, 2019. A distribution with the RMS of 0.6 PE
measures the non-uniformity between the units in WCDA-1. The monitoring of the
average charge of the peaks shows that the water transmission was continuously
changing over 179 days since Sept. 5, 2019 when WCDA-1 began the stable
operation phase. In the right panel of Fig.7, the average of the charge at
peaks has been plotted as a function of time. All the differences between
units and time variations have been corrected with respect to a “standard
detector” at Nov. 22, 2019 to maintain the uniformity and stability of the
detector.

Figure 6: The “spectral charge distribution” of cosmic rays measured by a
typical WCDA-1 unit. The clear “muon peak” is due to the diffuse single muons
falling in the unit.

Figure 7: (a): Map of charges at the “muon peak” measured by all WCDA-1 units
on Nov. 22, 2019. The average charge is 11.6 PE and the spread of charges is
0.6 PE among the 900 units. (b): The variation of the average charge at the
muon peak over 179 days from Sept. 5, 2019 to Feb. 29, 2020. This variation
indicates the change of water transmission over this period of time. The
dashed vertical line marks the ‘standard’ day Nov. 22, 2019.

With such a well-calibrated uniform and stable detector, the trigger rates of
detector units due to the high energy cosmic rays, instead of the muons,
should be the same among the 900 units and also stable. This has been tested
at the threshold $Q_{th}=2$ PE. It is found that the spread over 900 units is
within 2.3%. The stability of the rate is 1.3% over 170 days since Sept. 1,
2019. This test also shows a spread of 5.9% and stability of 2.9% at
$Q_{th}=0.3$ PE. This will be discussed in the systematic analysis below.

##  3 Observation of Gamma Rays from the Crab Nebula for WCDA-1 Performance
Testing

WCDA-1 started operation on April 26, 2019. In a short test run at the
beginning, a large volume of monitoring data was collected to validate the
operation. A trigger algorithm was implemented to record cosmic ray air
showers by requiring at least 20 groups, a group being 3$\times$3 detector
units, with each group being hit simultaneously in a window of 300 ns. A hit
is formed if any one of the nine 8” PMTs in the cluster received a signal
greater than the threshold of 1/3 PE. The detector has been fully operational
for stable survey of the northern sky since then. The air shower event rate is
about 45 kHz, while the trigger rate of a single detector unit is 22.5 kHz on
average. Once a trigger is formed, charges recorded by both 8” and 1.5” PMTs
in 2 $\mu$s are read out together with the absolute trigger time to build a
complete shower event. Further reconstruction of the air shower for its
geometrical properties is carried out after the event building. In Fig.8, the
daily live time of WCDA-1 is plotted as a function of time since April, 26,
2019, the day WCDA-1 was turned on. $\thicksim 1.7\times 10^{9}$ shower events
are recorded and successfully reconstructed everyday during the stable
operation for 179 days from Sept. 5, 2019 to Feb. 29, 2020. The analysis
reported in this paper uses the data taken in this period of time. Afterwards,
the operation was switched to the combined mode with both WCDA-1 and WCDA-2
participated.

Figure 8: The daily live time of WCDA-1 as a function of days since April 26,
2019. Two red vertical lines indicate the stable operation period, i.e. from
Sept. 5, 2019 to Feb. 29, 2020, in which data are used for the analysis in
this paper. About 1.7$\times 10^{9}$ shower events are recorded and
successfully reconstructed everyday during this period of time.

###  3.1 Reconstruction of Shower Geometry

For an air shower event, the arrival direction and core location are
reconstructed in two steps.

The first step is to find the arrival direction by assuming a planar shower
front using $N_{hit}$ units that are registered in the event with their
locations $(X_{i}$,$Y_{i})$ and trigger time $t_{i}$. The charge $Q_{i}$ is
used to weight the $i$-th unit in a defined way. The arrival direction is the
normal vector of the plane. Units that have a trigger time too far ($>100$ ns)
from the plane are considered to be random triggers by noise and rejected in
further steps with the remaining $N_{hit}$ units. The fitting procedure is
iterated until the number $N_{hit}$ of firing PMTs becomes constant.

The second step is to find the shower core location $(X_{c},Y_{c})$ and refine
the arrival direction reconstruction by assuming a conical shower front
centered at $(X_{c},Y_{c})$. The core location is found by calculating the
charge-weighted averages
$X_{c}=\Sigma^{N_{10}}_{i=1}Q_{i}X_{i}/\Sigma^{N_{10}}_{i=1}Q_{i}$ and
$Y_{c}=\Sigma^{N_{10}}_{i=1}Q_{i}Y_{i}/\Sigma^{N_{10}}_{i=1}Q_{i}$, where
$N_{10}=0.1N_{hit}$ and $i$ runs over the units that record the highest
charges. The conical front is defined by the time offsets $t(R)=c\cdot R$ with
respect to the plane found in the previous step, where $R$ is the distance
from the core in the plane and $c$ is the conical coefficient, set to be 0.07
ns/m for gamma ray induced showers in the working energy range of WCDA-1. At
this point, units that have trigger time farther than 30 ns from the conical
front are rejected again as noise triggers. $N_{hit}$ is renewed by counting
the surviving hits.

$N_{hit}$ is selected as a shower energy estimator in the analysis. For
showers in the WCDA-1 working energy range, most of the hits have only one or
two PEs. This enables a simple shower energy reconstruction based on the
nearly linear relationship between $N_{hit}$ and the primary energy. For the
pure electromagnetic cascade induced by gamma rays, the Monte-Carlo simulation
of the shower development is reliable to establish the energy response
function. The distribution of $N_{hit}$ of the well reconstructed shower
events is shown in Fig.9 as the curve in red. There is a clear pile-up of
events with $N_{hit}>840$ that indicates the saturation of $N_{hit}$, which
has a maximum of 900. This will be discussed in detail in Sec. 4.3 .

Figure 9: Distribution of the number of triggered detector units $N_{hit}$
after filtering out random noise hits (in red) and after the $C>10$ gamma
selection cut (in black). This cut causes a reduction of 97.2% for events with
$N_{hit}<$ 550.

###  3.2 Cosmic Ray Background Suppression

To maximize the sensitivity to gamma ray showers, an effective cosmic ray
background suppression is essential.

For showers induced by cosmic rays, mainly protons, the secondary particles
are more spread out laterally than for gamma-ray-induced showers. In either
case, secondary particles, mainly gamma rays and electrons, rapidly induce
sub-showers in water and lose all their energy over a few tens of centimetres.
All relativistic charged particles in the sub-showers generate Cherenkov light
in water with the emitting angle of 41.3∘. The Cherenkov photons in each
detector unit are sampled by the PMTs, thus forming a footprint of the shower
in WCDA-1. A typical cosmic ray event with $N_{hit}=189$ is shown in the right
panel of Fig.10. Each dot indicates a registered unit with the size indicating
the corresponding measured charge. A clear dense area indicated by the dashed
circle marks the core location of the shower; however, many large hits are
also found throughout the footprint. Those ‘hot spots’ are mainly due to muons
produced in the air shower which can easily penetrate through the entire depth
of WCDA-1, producing Cherenkov light along the whole tracks, thus yielding
stronger signals. Since muons are produced by the decay of $\pi^{\pm}$
generated in the shower, they carry the transverse momentum of pions so that
they are distributed over a wider area in general, when they arrive to the
ground. Therefore, on top of the smooth lateral distribution yielded by the
electromagnetic component, there are isolated ‘hot spots’ outside the core
region.

In contrast, a typical gamma-ray-induced shower from the direction of the Crab
Nebula as shown in the left panel of Fig.10 has a more compact and rather
smooth distribution of charges in the footprint measured by WCDA-1. With the
similar $N_{hit}=165$, the event has a better defined single core region
because the shower is a pure electromagnetic cascade process. Such clear
difference between the gamma ray signal and the cosmic ray background allows a
way to suppress the flux of cosmic rays in a specific gamma ray source region,
in particular for a large detector such as WCDA-1 that contains almost all of
the secondary particles in showers.

To characterize the features of showers induced by different primary
particles, a parameter named compactness was invented by the HAWC
collaboration [14]. The compactness is defined as $C=N_{hit}/{\rm
Max}(Q_{i};r>R_{c})$, where the maximum function finds the brightest unit
outside an area of radius $R_{c}$ centred at the shower core. The lack of hot
spots in the shower, in case of gamma rays, makes the compactness $C$
obviously larger than that of most background cosmic rays with the same
$N_{hit}$. For WCDA-1 at the altitude of 4410 m, $R_{c}$ of 45 meters is found
to be optimal, because the electromagnetic component of showers becomes dim
enough beyond $R_{c}$ to make any muon signal outstanding. By cutting on
$C>10$, the cosmic ray background can be significantly reduced. This is shown
in Fig.9, where the red data points show the $N_{hit}$ distribution for all
triggers (after filtering out random noise hits as described above) while the
black data points are for triggers with $C>10$, rejecting much of the cosmic-
ray background. In the range of $N_{hit}<$ 550, about 97.2% of cosmic ray
events are rejected. The rejection power gradually improves with the event
size and reaches 99.8% for large showers having 700 hits. To find a suitable
threshold for the parameter $C$, the Crab Nebula can be used as a standard
candle by maximizing the significance of its detection.

Figure 10: A high-confidence gamma-ray candidate event from the direction of
the Crab Nebula with $N_{hit}$=165 and a compactness value of $C=39$ (left
panel) and a similar background cosmic-ray event with $N_{hit}$=189 and
$C=1.0$ (right panel). Dashed circles indicate the reconstruction
uncertainties of the shower cores.

###  3.3 The Crab Nebula Detection and WCDA-1 Resolution and Sensitivity

The Crab Nebula is one of the strongest and most widely studied gamma-ray
sources. Its position, angular extension, flux and Spectral Energy
Distribution (SED) have been studied in detail. The coordinates of the Crab
Nebula in the equatorial system (RA, dec) are $(83.63^{\circ},22.02^{\circ})$.
The angular extension is much less than 0.07∘ in the TeV regime. Events
measured by WCDA-1 within an area of 6∘ $\times$ 6∘, centered at the Crab
nominal position, have been collected and filled in a grid of
$0.1^{\circ}\times 0.1^{\circ}$, after converting the arrival direction of the
shower from local coordinates $(\theta,\phi)$, zenith and azimuth angles, and
trigger time to equatorial coordinates. The background has been suppressed as
described above. Here, we tuned the threshold to be $C>10$.

In order to calculate the excess of signals from the Crab Nebula direction,
the remaining background of the cosmic rays is estimated using two different
methods for cross-checking. The first one is the so-called ‘direct
integration’ method, developed by the MILAGRO collaboration [15]. This method
is based on the assumption that acceptance of the detector, in local
coordinates, is independent of the trigger rate over a period of several hours
and that the cosmic ray background is isotropic. The last assumption is well
justified by the low-energy range considered in this analysis. The acceptance
map is properly normalized and transformed into an efficiency map, which is
applied to the events that passed the cuts, as measured at time $t$.

An alternative method is the so-called ‘equal zenith angle method’, developed
by the ASγ collaboration [16]. This method also relies on the assumption of
isotropy of the cosmic ray background. For each bin $(t,\theta,\phi)$ where
the source is supposed to be, the background is integrated over all other bins
of the same size at nearby azimuthal angles and the same zenith angle as the
source bin. In this analysis, the background is integrated over 6 azimuthal
bins around the one with the source. Therefore, the observation time ratio
between the source and background is 1/6. The underlying assumption for this
method is that the detector has a constant detection efficiency for showers
arriving at the same zenith angle. The bin size is defined by the Point Spread
Function (PSF) of the detector, which is well described by a 2-dimensional
Gaussian functional form with the angular resolution defined as $\Psi_{68}\sim
1.51\sigma=0.84^{\circ}$, where $\sigma$ is the standard deviation of the
Gaussian function and $\Psi_{68}$ is the angle of a cone with 68% of the
probability of the distribution contained, for the lowest energies in this
analysis. The cross check of background estimation between the two methods
showed no noticeable difference. Both the event map and background map in the
grid have to be smoothed using the PSF. The excess map is then obtained by
subtracting the background map from the event data map.

Figure 11: The 2-dimensional significance map around the position of The Crab
Nebula (left) and the distribution of the significance among $3600\times 900$
bins (right). Measured background cosmic ray events well distributed along a
standard normal distribution and the gamma ray signals are clearly distributed
over a wide region from $\sim 4\sigma$ to 77 $\sigma$.

The significance map around the Crab Nebula and the distribution of the events
as a function of the significance are shown in Fig. 11. The significance is
estimated following the method used in[7]. A clear Crab image with a
significance of 77.4$\sigma$ has been achieved. In order to calculate the
integrated sensitivity of WCDA-1, all events that have $N_{hit}>100$ collected
up to Feb. 29, 2020 are used including those before Sept. 5, 2019. An
integrated sensitivity of 65 milli-Crab-Unit (mCU) per year above 1 TeV (see
below) is confirmed by the observation. The position of the image center at
R.A.=$83.65^{\circ}\pm 0.04^{\circ}$ and Dec=$22.05^{\circ}\pm 0.04^{\circ}$
indicates that the pointing accuracy of WCDA-1, as a gamma ray telescope, is
better than $0.05^{\circ}$.

The PSF width of WCDA-1 depends on the shower size, i.e. the shower arrival
direction is better measured for larger showers since the shower front is
smoother with larger charge being detected at more sampling points. Therefore,
all events collected in the period from Sept. 5, 2019 to Feb. 29, 2020 around
the Crab Nebula are grouped into six bins of $N_{hit}$, namely
[60,100],[100,200],[200,300],[300,400],[400,500] and [500,800]. Fig.12 shows
the excesses in the six bins as a function of the space angle between the
events and the location of the Crab Nebula. This gives a clear measurement of
the PSF in terms of the space angle of the cone in which 68% of events are
contained as a function of shower size, namely 0.84∘,0.45∘,0.39∘,0.29∘,0.21∘
and 0.20∘ in the six bins, respectively. Events with $N_{hit}>800$ give a
measure of PSF of 0.20∘. Events with $N_{hit}<60$ are too small to be
reconstructed reliably and are not taken into account in the analysis reported
in this paper.

Figure 12: Distribution of excess events as a function of the space angle with
respect to the Crab Nebula direction in the six bins of $N_{hit}$. The solid
curves indicate 2D-Gaussian fits to the excesses, which result in the angular
resolution of 0.84∘,0.45∘,0.39∘,0.29∘,0.21∘ and 0.20∘ in the six bins,
respectively.

##  4 Monte Carlo Simulation for Gamma Shower Energy Reconstruction and WCDA-1
Acceptance Estimation

###  4.1 Air Shower Simulation

The air shower events are generated using CORSIKA v75000 [8]. Hadronic models
EPOS-LHC [10] and FLUKA [9] are selected for the interactions above and below
100 GeV, respectively. All shower particles have been tracked down to the
threshold energy of 50 MeV for hadrons and muons, and of 0.3 MeV for pions,
photons and electrons [11]. Five primary cosmic ray mass groups, i.e., H, He,
C-N-O, Mg-Al-Si and Fe, are simulated in the background estimation. The Poly-
gonato model [12] is used for the fluxes and spectral indices of the species.
The energy range of the incident cosmic rays is assumed to be from 10 GeV to 1
PeV for protons and helium and from 100 GeV to 1 PeV for the remaining three
groups. Cosmic rays are assumed to isotropically arrive in the zenith angle
range from 0∘ to 60∘. The energy spectrum of gamma rays is sampled from a
power law SED with the index of 2.71 for the signals from the Crab Nebula.
Every air shower event generated using CORSIKA is initiated at 100 random
locations in an area of $2000\times 2000\rm\;m^{2}$ centered at the center of
WCDA-1 and followed by the detailed simulation of the response of WCDA-1. More
than 1.5$\times 10^{8}$ comic ray events and 2.7$\times 10^{8}$ gamma ray
events have been generated in total.

###  4.2 Detector Response Simulation for WCDA-1

The code G4WCDA based on GEANT4 [13] has been developed to simulate the water
Cherenkov detector response in the LHAASO experiment. In order to investigate
the detector performance and correctly estimate the acceptance of WCDA-1, the
detailed geometry and structure of the pond, including its roof, steel
structure, beams, columns, concrete walls, detector suspension structure and
water transparency are all taken into account in G4WCDA. Since about 300
Cherenkov photons are produced by an electron or a muon per centimetre of
their tracks in water, most of the CPU time is consumed to simulate the huge
number of Cherenkov photons for their propagation, including scattering and
absorption in water. The propagation process is completed in two sequential
steps as sketched in Fig. 13 to improve the simulation efficiency. In the
first step, photons and secondary particles in sub-showers are traced without
taking into account the attenuation due to scattering or absorption. Only for
those that entered into a small volume of $80\times 80\times 80~{}{cm^{3}}$
around the PMTs in each detector unit, all parameters of the photons and
secondary particles are collected and stored in a ROOT file. In the second
step, the attenuation is taken into account for the stored photons. The
surviving photons and the secondary particles and newly generated photons are
further traced until they reach the photocathodes of PMTs. This approach has
proven very effective in coping with various detector operation conditions,
especially those associated with the water transmission which was changing in
the beginning phase of operation of WCDA-1.

Figure 13: Sketch of the two step simulation procedure in G4WCDA. Left panel:
Cherenkov photons/secondary particles were traced in water until entering into
the box in the vicinity of PMTs. Right panel: Only photons or secondary
particles that survived the attenuation cut were further traced until reaching
the photocathodes of PMTs.

The simulation is tested for its validity by comparing the distributions of
several shower parameters with the measured in the observation. In Fig. 14(a)
and (b), comparisons of the zenith ($\theta$) and azimuth ($\varphi$) angular
distributions between the simulated and measured showers with
$\theta<50^{\circ}$ are shown. Good agreement is evident, indicating that the
attenuation of showers in the atmosphere as well as the uniform response of
the detector array are well described. Besides the variables related to the
arrival direction, two shower energy-related variables are also compared as
shown in panels (c) and (d) of Fig.14, i.e., the total charge $Q$ and
$N_{hit}$ recorded by PMTs, respectively.

The shower core reconstruction described in Sec.3.1 has been verified using
the simulated events. The resolution is found to be 15 m, i.e. a circle with
radius of 15 m contains 68% of the content of the distribution of distances
between the thrown and reconstructed core locations, for showers having more
than 60 hits in the shower front.

In order to estimate the acceptance of WCDA-1 for gamma ray showers, the ratio
of gamma ray showers surviving the $C$-cut for cosmic ray background
suppression has to be estimated correctly. The measured and simulated $C$
distributions are compared in the left panel of Fig.15. The agreement between
data and the simulation indicates that the algorithm is valid to be used for
the gamma-ray-induced showers from the direction of the Crab Nebula. In the
right panel, the $C$ distribution of selected events is compared with that for
simulated showers with an assumption of a simple power law gamma ray spectrum
of $E^{-2.71}$. With a cut, e.g., $C>10$, a certain number of gamma ray
showers are cut while some cosmic ray background events still survive, thus a
signal-to-noise ratio yields a significance of the detection for the Crab
Nebula. In this way, the threshold of $C$ is tuned to maximize the
significance. Simultaneously, the efficiency of gamma ray detection or in
other words the effective area is also estimated in each bin of $N_{hit}$, as
listed in Table 1.

Figure 14: Comparison between simulated (MC) and measured distributions (a)
the distribution of zenith angle, (b) distribution of azimuth angle, (c)
distribution of total charge in PE and (d) the distribution of number of hits
$N_{hit}$ in the shower front. The vertical axes are event rates in arbitrary
units. In d), simulation results using two composition models due to
Horandel[12] in green and Gaisser et al.[22] in blue are compared with the
data.

Figure 15: Left panel: Distributions of compactness for data (in blue) and
simulated cosmic ray background events (in red). Right panel: The
distributions of compactness for selected gamma-like events and simulated
events assuming a simple power-law gamma-ray spectrum $\propto E^{-2.71}$.

###  4.3 Energy Reconstruction for Gamma Ray Induced Showers

As mentioned above, the number of hits found within 30 ns of the reconstructed
shower front, $N_{hit}$, is selected as the shower energy estimator. For
gamma-induced pure electro-magnetic showers, the simulated showers are used to
establish the correlation between the primary energy and $N_{hit}$. In Fig.16,
primary energy distributions for gamma ray showers in the 6 bins of $N_{hit}$
are plotted. Broad and quite overlapping distributions are observed. This
indicates that the energy resolution of the detection with WCDA-1 is not
perfect, in particular for showers whose major part may not be contained in
WCDA-1. With WCDA-2 and WCDA-3 merged in, the detector sensitive area will be
large enough to allow more strict event filtering, and the resolution could be
improved for the low energy showers in particular. Here, the medians of the
distributions are used as the measure of the gamma ray energy for showers in
the corresponding bins of $N_{hit}$ as listed in Table 1. The energy
resolution could be defined well as a symmetric Gaussian function above 6 TeV
where the resolution is about 33%. For lower energies, the resolution of
$log_{10}E$ is better defined by a Gaussian distribution, e.g. with a width of
0.5 decades around 1 TeV . For showers having more than 800 detector units
registered, the pile-up effect as shown in Fig. 9 makes the differential flux
measurement impossible due to the natural saturation of $N_{hit}$ at 900. An
integrated flux above the energy corresponding to $N_{hit}=800$ could be in
principle reported, however, the statistical significance is too weak,
$<2\sigma$, to be reliable for the flux measurement here. Moreover, including
a large number of units with too much charge recorded would require more
careful tuning in the simulation to better simulate the response of the
detector, as indicated in Fig. 14 (d). The high energy events include both
gamma-ray and cosmic ray induced showers, thus requiring a better energy
estimator than $N_{hit}$.

Figure 16: Normalized distributions of primary energy of gamma-ray showers from a source with spectral index of -2.71. The color scheme for curves is illustrated in the legend in the figure corresponding to 6 bins of $N_{hit}$. Table 1: Summary of data used in the measurement of SED of the Crab Nebula over $3.57\times 10^{6}$ seconds.  |  Nhit | $E_{med}$ | Excess | Background | Significance | Differential Flux
---|---|---|---|---|---|---
|  | $({\rm TeV})$ |  |  |  ($\sigma$) | $(cm^{-2}s^{-1}TeV^{-1})$
(a) | 60 - 100 | 0.58 | 1438.2 | 24885.8 | 9.1 | $(1.66\pm 0.20)\times 10^{-11}$
(b) | 100 - 200 | 1.1 | 1082.7 | 5202.3 | 15.0 | $(2.89\pm 0.23)\times 10^{-11}$
(c) | 200 - 300 | 2.4 | 456.2 | 1376.8 | 12.3 | $(4.74\pm 0.48)\times 10^{-12}$
(d) | 300 - 400 | 3.9 | 161.2 | 335.8 | 8.8 | $(1.12\pm 0.17)\times 10^{-13}$
(e) | 400 - 500 | 5.9 | 60.3 | 77.7 | 6.8 | $(3.54\pm 0.74)\times 10^{-13}$
(f) | 500 - 800 | 12.1 | 82.7 | 45.3 | 12.3 | $(6.91\pm 1.0)\times 10^{-14}$

##  5 Spectral Energy Distribution of the Crab Nebula

The SED of the Crab Nebula is measured using data collected from Sep. 5, 2019
to Feb. 29, 2020, for a live time of $3.57\times 10^{6}$ seconds. The criteria
for events to be used in the analysis are as follows.

  * •

the core is within a square of 80m$\times$80m centered at the WCDA-1 center.

  * •

$60\leq N_{hit}<800$.

  * •

The Crab Nebula was at least 45∘ above the horizon.

  * •

The compactness is larger than 10.

All the criteria helped to maximize the detection significance. The
corresponding significance, number of events exceeding the background and
remaining cosmic rays in each bin are listed in Table 1. Under the cuts,
WCDA-1 has an effective area that varies from about 400 $m^{2}$ at 600 GeV to
a constant of 3200 $m^{2}$ above 4.5 TeV as plotted in Fig. 17. Importantly,
this indicates that above 4.5 TeV the SED measurement has a constant
efficiency of about 50% mainly due to the gamma selection requiring $C>10$.

Figure 17: The effective area of WCDA-1 for the Crab detection with the cuts
applied in the analysis according to the simulation for gamma rays. The red
line indicates the physical area of the WCDA-1 used in the analysis.

The SED is determined by minimizing the $\chi^{2}$ function

$\displaystyle\chi^{2}=\sum_{i=1}^{6}\frac{\left[N_{i}^{obs}-N_{i}^{exp}(\phi_{0},\alpha,\beta)\right]^{2}}{(\sigma_{i}^{obs})^{2}}$
(1)

where $N_{i}^{obs}$ is the number of events exceeding the background in the
$i$-th bin of $N_{hit}$, $\sigma_{i}^{obs}$ is the statistical error of
$N_{i}^{obs}$ and $N_{i}^{exp}(\phi_{0},\alpha,\beta)$ is the expected number
of events according to the hypothesis of a log-parabolic model

$\displaystyle\phi(E)=\phi_{0}\left(\frac{E}{\rm 3\
TeV}\right)^{-\alpha-\beta~{}\rm ln\left(\frac{E}{\rm 3\ TeV}\right)}$ (2)

with three parameters ($\phi_{0}$, $\alpha$, $\beta$), where $E$ is shower
energy and $\phi$ is the photon flux from the Crab Nebula. This formula was
suggested as one of the best functional approximations for the Crab nebula
spectrum over a wide energy range[23]. The number of events for signals from
the Crab Nebula direction is shown in Fig.12 for each bin of $N_{hit}$ with
the background cosmic ray events subtracted. They can be fit quite well with a
2-dimensional Gaussian functional template.

Figure 18: The spectral energy distribution of the Crab Nebula in the energy
range from 800 GeV to 13 TeV measured by WCDA-1 of LHAASO is shown as the red
filled triangles. As a comparison, the SED measured by using KM2A of
LHAASO[24] above 10 TeV is also shown as the green filled circles. The error
bars are statistic errors and shade area represents the systematic error.
Measurements of other experiments[18][19][17][21][20] are also plotted for
comparison.

The SED measured using WCDA-1 is shown in Fig. 18 as the red points. The
corresponding differential flux in each bin is listed also in Table 1. he
spectral parameters obtained in the log-parabolic fitting are
$\phi_{0}=(2.32\pm 0.19)\times 10^{-12}cm^{-2}~{}TeV^{-1}~{}s^{-1}$,
$\alpha=2.57\pm 0.06$, and $\beta=0.02\pm 0.05$ with $\chi^{2}/ndf$=0.513,
respectively.

For comparison, the SED measured by other ground-based gamma ray experiments,
i.e., HAWC [17], Tibet ASγ [19] and ARGO-YBJ [18], are shown in the same
figure, together with the measurements by the IACT experiments, HESS [21] and
MAGIC[20].

The systematic errors are mainly caused by following issues. Firstly, WCDA-1
has been continuously improving in terms of the water transparency and
detector stability during the observation reported in this paper. This is
reflected by Fig.7 (b) that shows the measured average PMT charge generated by
single muons nearly monotonically increased with time except during major
maintenance in October 2019. This requires a correction of about $\pm$11% to
the charge measurements, thus introducing an uncertainty in the detection
efficiency for events with a small number of hits, namely about 8.8% at
$Q=0.3$ PE, including non-uniformity among 900 detector units and variation
over 170 operational days, tested by using cosmic rays which are uniformly and
hitting the detector from all directions.. Secondly, the shower and detector
simulation is an essential tool in estimation of the detector acceptance in
terms of effective area for gamma ray detection. There is no way to check its
validity directly using gamma rays since we do not have a known gamma ray beam
for the testing. However, it can be checked using cosmic rays. Assuming that
the central region of WCDA-1, e.g., a circle with radius of 30 m, is fully
efficient for large enough showers, e.g., showers with more than 500 hits, one
can calibrate the absolute flux of cosmic rays calculated using the simulation
codes in comparison with observations. Here, the details about the cosmic ray
composition and the models of interactions with the atmospheric nuclei in the
relevant energy range below 100 TeV have to be assumed as well. As shown in
Fig.14, there are discrepancies between simulated and measured distributions.
Here two composition assumptions are compared with data in the $N_{hit}$
distributions. This results in an uncertainty of the estimate of detector
efficiency associated with the assumptions, including the entangled issues of
interaction modeling and cosmic ray composition. The larger discrepancy of 16%
between the two assumptions is taken as a conservative estimate of the
corresponding uncertainty of the detector efficiency. The uncertainty of gamma
ray detection acceptance estimated by using the simulation tool is assumed to
be the same in this analysis. Thirdly, the gamma ray shower energy scale also
totally relies on the simulation. As shown in Fig.14 (c), the simulation of
the detector response near the threshold is slightly overestimated. In order
to estimate the associated uncertainty of the energy scale, energies of the
gamma events from the Crab Nebula direction have been reconstructed with
different thresholds of $Q_{th}$, 0.3 PE and 0.6 PE, respectively. A
difference of 3% is found between the average energies due to the two
reconstructions; thus a uncertainty of 8% is expected in the SED. The overall
systematic uncertainty could be as large as ${}^{+8}_{-24}$%. It is presented
in Fig.18 as the shaded area around the SED measured by WCDA-1.

##  6 Conclusions

In this paper, the performance of the first LHAASO water Cherenkov detector,
WCDA-1, is tested by observing the Crab Nebula as a standard candle of TeV
gamma ray astronomy. The Crab Nebula has been detected with a significance of
77.4 $\sigma$, corresponding to a sensitivity of 65 mCU, in agreement with the
design specification. The angular resolution for the gamma ray arrival
direction is found to be better than 0.4∘ above 3 TeV and the pointing
accuracy to be better than 0.05∘. This is a result of well-measured detector
orientation and timing calibration of every detector unit. Air shower and
detector simulations are well tuned to reproduce the measured cosmic ray data
and used for the gamma-ray shower energy reconstruction and detector
acceptance estimation. The spectral energy distribution of the Crab Nebula has
been measured in the TeV range from 0.5 TeV to 15.8 TeV showing a good
agreement with all other major gamma ray astronomy experiments.

###### Acknowledgements.

Acknowledgments The authors would like to thank all staff members who work at
the LHAASO site above 4400 meter above the sea level year round to maintain
the detector and keep the water recycling system, electricity power supply and
other components of the experiment operating smoothly. We are grateful to
Chengdu Management Committee of Tianfu New Area for the constant financial
support for research with LHAASO data. This research work is also supported by
the following grants: The National Key R&D program of China under grants
2018YFA0404201, 2018YFA0404202 and 2018YFA0404203, by the National Natural
Science Foundation of China (NSFC grants No.12022502, No.11905227,
No.U1931112, No.11635011, No.11761141001, No.Y811A35, No.11675187,
No.U1831208), and in Thailand by RTA6280002 from Thailand Science Research and
Innovation.

## References

  * [1]
  * [2] Z. Cao, Chin. Phys. C, 34 (2)(2010): 249–252; H. He (for LHAASO Coll.), Radiation Detection Technology and Methods 2-7, (2018)
  * [3] G. Gong, S. Chen, Q. Du, J. Li, Y. Liu, H. He, et al., ‘Sub-nanosecond timing system design and development for LHAASO Project’ in Proceedings of ICALEPCS2011, Grenoble, France, 2011.
  * [4] J.Y.Liu, et al, JINST, 12(2017):10021.
  * [5] X. Li, et al, Nucl. Inst. and Methods in Phys. Res. A, 912 (2018): 382-385
  * [6] Hui-Cai Li, et al, Chinese Physics C, Volume 41(2017) No.2 026002.
  * [7] B. Bartoli, et al, ARGO-YBJ Collaboration, ApJ, 779 (2013) 27.
  * [8] http://www-ik.fzk.de/corsika.
  * [9] G. Battistoni, et al., AIP Conference Proceedings 896, (2007) 31.
  * [10] T. Pierog, et al., Phys. Rev. C 92, (2015) 034906.
  * [11] H.C. Li, et al, Chin. Phys. C, 38 (1)(2014): 016002.
  * [12] J.R. Horandel, Astroparticle Physics, 19 (2003): 193–220.
  * [13] S. Agostinelli, et al, Geant4 Collaboration, Nucl. Inst. Methods Phys. Res. A, 506(2003): 250–303.
  * [14] A. U. Abeysekara, et al, HAWC collaboration, Astropart. Phys. 50-52 (2013), 26-32.
  * [15] MILAGRO Collaboration, APJ, 595 (2) (2003) 803-811.
  * [16] M. Amenomori, et al, AS-Gamma Collaboration, ApJ, 525 (1999), L93-L96.
  * [17] A. U. Abeysekara, et al, HAWC Collaboration, ApJ, 881 (2019) 134.
  * [18] B. Bartoli, et al, ARGO-YBJ Collaboration, ApJ, 798 (2015) 119.
  * [19] M. Amenomori, et al, AS-Gamma Collaboration, ApJ, 692 (2009) 61.
  * [20] J.Aleksić, et al, MAGIC Collaboration, Journal of High Energy astrophysics, V5-6 (2015) p30-38 ( arXiv:1406.6892v5 [astro-ph.HE]).
  * [21] F.Aharonian, et al, HESS Collaboration, A&A, 457 (2006) 899-915.
  * [22] T. K. Gaisser, T. Stanev and S. Tilav, arXiv:1303.3565v1 [astro-ph.HE].
  * [23] J.Aleksić, et al, MAGIC Collaboration, Astroparticle Physics,72 (2016) 76.
  * [24] F. Aharonian, et al (LHAASO Coll.), Chin. Phys. C45, in press (2020), arXiv: 2010.06205.


#  Wave function of ${}^{9}\mathrm{Be}$ in the three-body $(\alpha\alpha
n)$-model

S. A. Rakityansky${}^{\,a,b}$
aJoint Institute for Nuclear Research, Dubna, Russia
bDepartment of Physics, University of Pretoria, Pretoria, South Africa

###### Abstract

A simple analytic expression of the three-body wave function describing the
system $(\alpha\alpha n)$ in the ground state $\frac{3}{2}^{-}$ of
${}^{9}\mathrm{Be}$ is obtained. In doing this, it is assumed that the
$\alpha$ particles interact with each other via the $S$-wave Ali-Bodmer
potential including the Coulomb term, and the neutron-$\alpha$ forces act only
in the $P$-wave state. This wave function is constructed by trial and error
method via solving in this way a kind of inverse problem when the two-body
$\alpha\alpha$ potential is recovered from a postulated three-body wave
function. As a result, the wave function is an exact solution of the
corresponding three-body Schrödinger equation for experimentally known binding
energy and for the $\alpha\alpha$ potential whose difference from the Ali-
Bodmer one is minimized by varying the adjustable parameters which the
postulated wave function depends on.

##  1 Introduction

Beryllium isotope ${}^{9}\mathrm{Be}$ is of special interest from several
points of view. First of all, it is clusterized with very high probability in
two $\alpha$ particles and a neutron. As a result the low-lying spectrum of
this nucleus is well-reproduced within the three-body ($\alpha\alpha n$)-model
[1, 2, 3, 4, 5].

Secondly, neither $\alpha\alpha$ nor $\alpha n$ pairs can form bound states,
which means that this isotope is an example of a Borromean system with a
rather small binding energy of $1.5736\,\mathrm{MeV}$ [6].

Thirdly, due to such a weak three-body binding the valence neutron is loosely
bound as well. This means that the neutron should move rather far away from
the centre of mass, which makes the nucleus ${}^{9}\mathrm{Be}$ a transient
one between ordinary compact nuclei and the halo-nuclei.

And finally, the last but not least is the fact that ${}^{9}\mathrm{Be}$ plays
an important role in astrophysical processes, namely, in synthesis of heavy
elements in the universe. Since there are no stable nuclei with $A=5$ and
$A=8$, the two-body $pp$-chain reactions practically stop at the formation of
${}^{4}\mathrm{He}$ [7, 8]. The synthesis of more heavy elements requires a
bridge over these so-called $A=5,8$ mass gaps. These gaps can be crossed over
via the three-body fusion reactions. Among them the most well known is the
triple-alpha fusion, $\alpha\alpha\alpha\rightarrow{}^{12}\mathrm{C}+\gamma$,
which is usually associated with the Hoyle resonance (see, for example, Ref.
[9]). In a neutron-rich environment one of the alternative bridges is provided
by the radiative process $n(\alpha\alpha,\gamma){}^{9}\mathrm{Be}$ which
involves the Beryllium isotope that is considered in the present paper.

Of course, in addition to what was already said, the nucleus
${}^{9}\mathrm{Be}$ appears either in the initial or in the final states of
various nuclear reactions, such as
${}^{9}\mathrm{Be}(n,\gamma){}^{10}\mathrm{Be}$ [10],
${}^{9}\mathrm{Be}({}^{18}\mathrm{O},{}^{17}\mathrm{O}){}^{10}\mathrm{Be}$
[11], ${}^{8}\mathrm{Li}(p,\gamma){}^{9}\mathrm{Be}$ [12, 13], etc.

When theoretically describing the processes mentioned above, one usually needs
the wave functions of the quantum systems involved, and in particular the wave
function of ${}^{9}\mathrm{Be}$. There are many different approaches to
obtaining such a function.

Strictly speaking, this nucleus should be considered as a nine-body system.
Such an (ab initio) approach can be based either on the shell model or on the
expansion over the hyperspherical harmonics. However, as was mentioned before,
${}^{9}\mathrm{Be}$ is clusterized in two $\alpha$ particles and a neutron.
Thanks to the tight binding of ${}^{4}\mathrm{He}$ ($\sim 28\,\mathrm{MeV}$),
the $\alpha$ particles move inside ${}^{9}\mathrm{Be}$ like solid bodies
without internal excitations, and therefore the nine-body problem can be
reduced to an effective three-body one. To the best of the author’s knowledge
all the calculations that were done for ${}^{9}\mathrm{Be}$ during the last
few decades, exploited the $(\alpha\alpha n)$ cluster-representation. The
three-body problem can also be solved in different ways: either using exact
Faddeev equations (see, for example, Refs. [2, 3]), or with the help of
various approximate methods [4, 14, 15, 16].

The common feature of all the publications, where these three-body approaches
are realised, is that the wave function of ${}^{9}\mathrm{Be}$ is not given
(and in most cases cannot be given) in such a form that could be used by the
other people in their own calculations. In other words, if somebody wants to
use the same wave function, he or she has to repeat all the complicated
calculations described in these publications. However, for the majority of the
researchers this is a difficult obstacle. Indeed, the procedure of solving,
for example, the Faddeev equations involves many mathematical and numerical
tricks that are only known to those who specialize in the field of the few-
body problem. The shell model and the hyperspherical expansion are not much
easier. It is therefore desirable to have some parametrized analytic
expressions of the wave functions of various nuclei that are easy to use. This
is what the present paper is devoted to. Here such a parametrization is
obtained for the wave function of ${}^{9}\mathrm{Be}$ in the three-body
cluster model $(\alpha\alpha n)$.

The way of obtaining the parametrized three-body wave function is based on the
work by Belyaev et al. [17]. That parer deals with a three-body problem in
which two of the three pairwise potentials as well as the three-body wave
function are given while one of the two-body potentials is unknown. The
authors of Ref. [17] developed a method for recovering this unknown potential.

In the present work, this method is kind of “reversed”, i.e. the wave function
is postulated in a parametrized form and its parameters are optimized by
minimizing the difference between the recovered potential and the
corresponding exact one. For the system under consideration, $(\alpha\alpha
n)$, it is chosen the $\alpha\alpha$ potential for the role of recovered one.
The analytic form of the wave function and its parameters are chosen to fit
(as close as possible) the Ali-Bodmer potential [18].

##  2 Formalism

Consider a bound state of two $\alpha$ particles and a neutron with a negative
three-body energy $E_{b}$. It is known [6] that this system has only one bound
state, which is the nucleus ${}^{9}\mathrm{Be}$. Using the Jacobi coordinates
shown in Fig. 1, the three-body Schrödinger equation for the wave function
$\Psi(\vec{r}_{\alpha\alpha},\vec{r}_{n})$ describing ${}^{9}\mathrm{Be}$, can
be written as follows:

$\displaystyle
V_{\alpha\alpha}({r}_{\alpha\alpha})\Psi(\vec{r}_{\alpha\alpha},\vec{r}_{n})=$
$\displaystyle=$ $\displaystyle\left[E_{b}-V_{\alpha n}(\rho_{1})-V_{\alpha
n}(\rho_{2})+\frac{\hbar^{2}}{2\mu_{\alpha\alpha}}\Delta_{\vec{r}_{\alpha\alpha}}+\frac{\hbar^{2}}{2\mu}\Delta_{\vec{r}_{n}}\right]\Psi(\vec{r}_{\alpha\alpha},\vec{r}_{n})\
,$

where $V_{\alpha\alpha}$ and $V_{\alpha n}$ are the two-body potentials for
the $\alpha\alpha$ and $\alpha n$ pairs, $\mu_{\alpha\alpha}=m_{\alpha}/2$ and
$\mu=2m_{\alpha}m_{n}/(2m_{\alpha}+m_{n})$ are the reduced masses
corresponding to the respective Jacobi coordinates, the radial variables
$\rho_{1}$ and $\rho_{2}$ given by

$\rho_{1,2}=\sqrt{\displaystyle\frac{1}{4}r_{\alpha\alpha}^{2}+r_{n}^{2}\mp
r_{\alpha\alpha}r_{n}\cos\vartheta}\ ,$ (2)

are the distances between the neutron the the $\alpha$ particles, and the
Laplacians involve the derivatives with respect to the corresponding Jacobi
vectors.

Figure 1: Jacobi coordinates, $\vec{r}_{\alpha\alpha}$ and $\vec{r}_{n}$, that
specify a space configuration of the system $(\alpha\alpha n)$. The symbols
$\rho_{1}$ and $\rho_{2}$ denote distances between the neutron and the two
alpha particles.

Multiplying Eq. (2) by $\Psi^{\dagger}(\vec{r}_{\alpha\alpha},\vec{r}_{n})$
from the left and integrating over vector $\vec{r}_{n}$ as well as over the
spherical angles, $\Omega_{\alpha\alpha}$, of vector $\vec{r}_{\alpha\alpha}$,
one can find the potential $V_{\alpha\alpha}$, if the energy $E_{b}$, the
potential $V_{\alpha n}$, and the function
$\Psi(\vec{r}_{\alpha\alpha},\vec{r}_{n})$ are given,

$\displaystyle V_{\alpha\alpha}({r}_{\alpha\alpha})=E_{b}-$ $\displaystyle-$
$\displaystyle\frac{1}{D({r}_{\alpha\alpha})}\int
d\Omega_{\alpha\alpha}d^{3}r_{n}\Psi^{\dagger}(\vec{r}_{\alpha\alpha},\vec{r}_{n})\left[\\!\\!\phantom{\frac{|}{|}}V_{\alpha
n}(\rho_{1})+V_{\alpha n}(\rho_{2})-\right.$ $\displaystyle-$
$\displaystyle\left.\frac{\hbar^{2}}{2\mu_{\alpha\alpha}}\Delta_{\vec{r}_{\alpha\alpha}}-\frac{\hbar^{2}}{2\mu}\Delta_{\vec{r}_{n}}\right]\Psi(\vec{r}_{\alpha\alpha},\vec{r}_{n})\
,$

where

$D({r}_{\alpha\alpha})=\int
d\Omega_{\alpha\alpha}d^{3}r_{n}\Psi^{\dagger}(\vec{r}_{\alpha\alpha},\vec{r}_{n})\Psi(\vec{r}_{\alpha\alpha},\vec{r}_{n})$
(4)

is the relative $\alpha\alpha$ radial probability density.

It should be noted that the wave function
$\Psi(\vec{r}_{\alpha\alpha},\vec{r}_{n})$ is (by construction) the exact
solution of the Schrödinger equation with a given $E_{b}$ and with the
potential $V_{\alpha\alpha}$ which is obtained from Eq. (2). This means that,
if one chooses an approximate $\Psi$, it is still an exact solution of Eq. (2)
with the corresponding approximate $V_{\alpha\alpha}$. The difference of this
approximate $\alpha\alpha$ potential from the one, which is considered as the
exact potential, characterizes the quality of the chosen $\Psi$. If one
manages to find such a function $\Psi$ that the potential, which follows from
Eq. (2), is almost the same as the exact one, one actually finds a reliable
solution of Eq. (2) for a given energy $E_{b}$. This can be achieved by
varying the parameters of a postulated function $\Psi$ in order to minimize
the difference between the approximate and the exact $\alpha\alpha$
potentials.

##  3 Two-body potentials

The main input information that is needed in the present work is given by the
$\alpha\alpha$ and $\alpha n$ two-body potentials. As is explained in the next
section, these potentials are needed for the $S$-wave and $P$-wave states,
respectively. The first of them is taken from Ref. [18]. This is well known
Ali-Bodmer $\alpha\alpha$-potential:

$V_{\alpha\alpha}(r)=V_{R}e^{-(\beta_{R}r)^{2}}-V_{A}e^{-(\beta_{A}r)^{2}}+V_{c}(r)\
,$ (5)

which involves the repulsive and attractive Gaussian terms as well as the
repulsive electric (Coulomb-like) one,

$V_{c}(r)=\frac{4e^{2}}{r}\mathop{\mathrm{erf}}\left(\frac{\sqrt{3}}{2R_{\alpha}}r\right)\
.$ (6)

The parameters of this potential are: $V_{R}=1050$ MeV,
$\beta_{R}=0.8\,\mathrm{fm}^{-1}$, $V_{A}=150$ MeV,
$\beta_{A}=0.5\,\mathrm{fm}^{-1}$, $R_{\alpha}=1.44\,\mathrm{fm}$. The Coulomb
term (6) takes into account a non-zero size of the $\alpha$ particle. The
expression (6) is obtained under the assumption that its charge has a Gaussian
distribution in space with the rms-radius $R_{\alpha}$.

The $\alpha n$ potential in the present work is the same that was suggested by
Bang and Gignoux [19]. It is of the Woods-Saxon type with a spin-orbit term:

$V_{\alpha
n}(r)=\frac{W_{1}}{1+\exp\left[(r-R_{1})/d_{1}\right]}+\frac{\vec{\ell}\cdot\vec{s}}{r}\frac{d}{dr}\frac{W_{2}}{1+\exp\left[(r-R_{2})/d_{2}\right]}\
,$ (7)

where $W_{1}=-43.0\,\mathrm{MeV}$, $R_{1}=2.0\,\mathrm{fm}$,
$d_{1}=0.7\,\mathrm{fm}$, $W_{2}=40.0\,\mathrm{MeV\cdot\mathrm{fm}^{2}}$,
$R_{2}=1.5\,\mathrm{fm}$, $d_{2}=0.35\,\mathrm{fm}$. For the quantum state
with the total angular momentum $J=3/2$, orbital angular momentum $\ell=1$,
and the neutron spin $s=1/2$, the action of the operator
$\vec{\ell}\cdot\vec{s}$ is equivalent to multiplication by 1/2.

The potentials (5, 7) were used by many authors in a number of calculations
within the three-body models of ${}^{9}\mathrm{Be}$ and ${}^{6}\mathrm{Li}$
(see, for example, Refs. [2, 3, 19] where the corresponding three-body Faddeev
equations were solved). It turned out that these potentials reliably described
the $\alpha\alpha$ and $\alpha n$ interactions at low energies as well as
allowed one to well reproduce the binding energies and the other properties of
the nuclei. The potentials (5, 7) are therefore used in the present paper as
well.

##  4 Wave function of ${}^{9}\mathrm{Be}$

As was explained in the Introduction, the Beryllium isotope
${}^{9}\mathrm{Be}$ is clusterized in two $\alpha$ particles and a neutron.
This is a Borromean three-body system that has only one bound state with a
rather small binding energy $|E_{b}|=1.5736\,\mathrm{MeV}$ and with the spin-
parity $(3/2)^{-}$ [6]. Since the $\alpha$ particle spins are zero, the three-
body wave function should be symmetric with respect to the $\alpha$-$\alpha$
permutation. This means that their relative orbital angular momentum is even,
$\ell_{\alpha\alpha}=0,2,4,\dots$. The negative parity can therefore be
obtained if the valence neutron moves with an odd orbital angular momentum,
$\ell_{n}=1,3,5,\dots$, relative to the $\alpha\alpha$ pair.

As a reliable approximation, it is reasonable to assume that
$\ell_{\alpha\alpha}=0$ and $\ell_{n}=1$. All the higher partial waves are
ignored since the $(\alpha\alpha n)$ binding is weak and Borromean while high
values of $\ell_{\alpha\alpha}$ and $\ell_{n}$ introduce additional repulsive
centrifugal potentials into the three-body hamiltonianian. The minimal
repulsion is in the state $\ell_{\alpha\alpha}=0$, $\ell_{n}=1$.

One may wonder why $\ell_{n}=1$ and not zero, which would generate even less
repulsion. An intuitive explanation can be found by considering the naive
shell model as follows. There are five neutrons in this nucleus. Four of them
are sitting in the clusters. Each cluster has its own mean-field potential and
its own level structure. Two neutrons in the first $\alpha$ particle occupy
its lowest $S$-level and the other two neutrons occupy the lowest $S$-level in
the second $\alpha$ particle. What remains for the fifth neutron is the
$P$-level either in the first or in the second cluster. Therefore the valence
neutron moves in the $P$-wave state with respect to both $\alpha$ particles.
This is therefore the shell-model configuration with the minimal energy. And
since $\ell_{\alpha\alpha}=0$, the orbital angular momentum of the neutron
with respect of the centre of mass of the $\alpha\alpha$ subsystem is also
one, i.e. $\ell_{n}=1$.

Based on the above reasoning, one comes to the following structure of the wave
function of ${}^{9}\mathrm{Be}$:

$\Psi(\vec{r}_{\alpha\alpha},\vec{r}_{n})=R(r_{\alpha\alpha},r_{n})Y_{00}(\Omega_{\alpha\alpha})\mathcal{Y}_{1\frac{1}{2}}^{\frac{3}{2}J_{z}}(\Omega_{n})\
,$ (8)

where the spin-angular function of the neutron,
$\mathcal{Y}_{1\frac{1}{2}}^{\frac{3}{2}J_{z}}(\Omega_{n})$, defined as

$\mathcal{Y}_{\ell s}^{JJ_{z}}(\Omega)=\sum_{ms_{z}}C_{\ell
mss_{z}}^{JJ_{z}}Y_{\ell m}(\Omega)\chi_{s}(s_{z})\ ,$ (9)

couples its orbital angular momentum and the spin. Here $\chi_{s}$ is the spin
function of the neutron.

Substituting wave function (8) in Eqs. (2) and (4), the following
decomposition of the $\alpha\alpha$ potential is obtained:

$V_{\alpha\alpha}({r}_{\alpha\alpha})=E_{b}-\left\langle V_{\alpha
n}(\rho_{1})\phantom{\frac{1}{1}}\\!\\!\\!+V_{\alpha
n}(\rho_{2})\right\rangle-\left\langle\frac{\ell_{n}(\ell_{n}+1)}{r_{n}^{2}}\right\rangle+\left\langle
E_{\mathrm{kin}}\right\rangle\ ,$ (10)

where the averaging, denoted by the brackets $\langle\rangle$, is done over
all the configuration space variables except the radial variable
$r_{\alpha\alpha}$. These average values represent the contributions that are
originated from the two $\alpha n$ potentials,

$\left\langle V_{\alpha
n}(\rho_{1,2})\phantom{\frac{1}{1}}\\!\\!\\!\right\rangle=\frac{3}{4D(r_{\alpha\alpha})}\int_{0}^{\infty}\\!dr_{n}\int_{0}^{\pi}\\!d\vartheta\,r_{n}^{2}(\sin\vartheta)^{3}V_{\alpha
n}(\rho_{1,2})\left|R(r_{\alpha\alpha},r_{n})\right|^{2}\ ,$ (11)

from the centrifugal potential caused by $\ell_{n}=1$,

$\left\langle\frac{\ell_{n}(\ell_{n}+1)}{r_{n}^{2}}\right\rangle=\frac{\hbar^{2}}{\mu
D(r_{\alpha\alpha})}\int_{0}^{\infty}\\!dr_{n}\left|R(r_{\alpha\alpha},r_{n})\right|^{2}\
,$ (12)

and from the radial derivatives (“kinetic energy” terms) of the Laplacians
$\Delta_{\vec{r}_{\alpha\alpha}}$ and $\Delta_{\vec{r}_{n}}$,

$\displaystyle\left\langle E_{\mathrm{kin}}\right\rangle$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{2\mu_{\alpha\alpha}r_{\alpha\alpha}^{2}D(r_{\alpha\alpha})}\int_{0}^{\infty}\\!dr_{n}\,r_{n}^{2}R^{*}(r_{\alpha\alpha},r_{n})\frac{\partial}{\partial
r_{\alpha\alpha}}\left[r_{\alpha\alpha}^{2}\frac{\partial}{\partial
r_{\alpha\alpha}}R(r_{\alpha\alpha},r_{n})\right]$ (13) $\displaystyle+$
$\displaystyle\frac{\hbar^{2}}{2\mu
D(r_{\alpha\alpha})}\int_{0}^{\infty}\\!dr_{n}\,R^{*}(r_{\alpha\alpha},r_{n})\frac{\partial}{\partial
r_{n}}\left[r_{n}^{2}\frac{\partial}{\partial
r_{n}}R(r_{\alpha\alpha},r_{n})\right]\ ,$

where

$D(r_{\alpha\alpha})=\int_{0}^{\infty}\\!dr_{n}\,r_{n}^{2}\left|R(r_{\alpha\alpha},r_{n})\right|^{2}\
.$ (14)

When obtaining Eq. (11) it was used the fact that the left hand side of Eq.
(10) cannot depend on $J_{z}$. The simplest choice is $J_{z}=3/2$ which
implies that the sum in Eq. (9) only includes one term involving the spherical
harmonics
$Y_{11}(\vartheta,\varphi)=-(1/2)\sqrt{3/(2\pi)}\sin\vartheta\exp(i\varphi)$,
where the $z$-axis is chosen along $\vec{r}_{\alpha\alpha}$ and the polar
angle $\vartheta$ is between vectors $\vec{r}_{\alpha\alpha}$ and
$\vec{r}_{n}$ (see Fig. 1).

It should be noted that in Eqs. (11-13) the normalization of the radial wave
function $R(r_{\alpha\alpha},r_{n})$ is arbitrary. This is because in each of
these equations the square of $R$ is present in both the numerator and
denominator. This fact significantly simplifies the procedure for fitting the
potential $V_{\alpha\alpha}({r}_{\alpha\alpha})$, because one does not have to
care about proper nomalization of $R(r_{\alpha\alpha},r_{n})$ when varying its
parameters. The function $R(r_{\alpha\alpha},r_{n})$ only needs to be
normalized to unity at the final stage when all its parameters have been
optimized and established.

The choice of a functional form of $R(r_{\alpha\alpha},r_{n})$ is the most
important and difficult part of the described procedure. It was done by trial
and error approach. After many unsuccessful attempts, it was found such a form
of this radial wave function that Eq. (10) produced the potential curve whose
shape was more or less similar to the Ali-Bodmer potential. Finally, the
parameters of $R(r_{\alpha\alpha},r_{n})$ were adjusted by minimizing the
difference between the Ali-Bodmer potential and the one produced by Eq. (10).
This was done via the least-square method using the minimization code MINUIT
[21].

The resulting radial wave function of ${}^{9}\mathrm{Be}$ looks as follows:

$R(r_{\alpha\alpha},r_{n})=\frac{N}{r_{\alpha\alpha}^{2}r_{n}^{3}}\arctan\left[(a_{1}r_{\alpha\alpha})^{7/2}\right]\arctan\left[(a_{2}r_{n})^{7/2}\right]\left[1-e^{-(a_{3}r_{\alpha\alpha})^{2}}\right]^{2}\frac{e^{-\varkappa\zeta}}{\zeta^{5/2}}\
,$ (15)

where $\varkappa=\sqrt{2\mu|E_{b}|/\hbar^{2}}$ is the momentum corresponding
to the experimental binding energy,

$\zeta=\sqrt{\frac{\mu_{\alpha\alpha}}{\mu}r_{\alpha\alpha}^{2}+r_{n}^{2}}$
(16)

is the hyperradius, $a_{1}=0.37441\,\mathrm{fm}^{-1}$,
$a_{2}=0.070112\,\mathrm{fm}^{-1}$, $a_{3}=0.41953\,\mathrm{fm}^{-1}$, and the
normalization constant $N=927446.734\,\mathrm{fm}^{9/2}$.

A specific choice of the functional form of $R(r_{\alpha\alpha},r_{n})$ was
based on some general properties of such a function as well as on the author’s
intuition. Unfortunately, there is no universal recipe for constructing this
function.

The following general features of the wave function were taken into account.
First of all, the radial wave function must exponentially tend to zero at
large distances in all directions in the configuration space. It is known from
the hyperspherical expansion theory that for any three-body system such an
exponential diminishing is given by the factor
$e^{-\varkappa\zeta}/\zeta^{5/2}$. This is why this factor is present in the
expression (15). Due to the Pauli blocking, the $\alpha$ particles cannot
penetrate each other. For the same reason the neutron cannot move inside them.
Therefore the Pauli blocking requires vanishing of the wave function when
$r_{\alpha\alpha}$ and $r_{n}$ tend to zero. And finally, the wave function
should have maxima at some short distances along both radial variables,
$r_{\alpha\alpha}$ and $r_{n}$. These maxima correspond to the most probable
configuration of the system.

A maximum at some intermediate distance along the variable $x$ as well as
vanishing near $x=0$ can be implemented by using the function
$\arctan[(ax)^{n}]/x^{m}$ with $n>m$. There are two such factors in the
function (15), for $x=r_{\alpha\alpha}$ and for $x=r_{n}$. In principle, the
powers $n$ and $m$ could be made additional adjustable parameters. However, it
was decided not to increase the number of the free parameters, because a
multi-variable function usually has many local minima and threfore with too
many variable parameters the fitting of the potential could be more difficult.
The optimal values of $n$ and $m$ for the variables $r_{\alpha\alpha}$ and
$r_{n}$ were found manually by the trial and error method.

It turned out that the suppression of $R(r_{\alpha\alpha},r_{n})$ at
$r_{\alpha\alpha}\to 0$ that was provided by the corresponding factor
$\sim\arctan(.)$, was insufficient (the resulting $V_{\alpha\alpha}$ potential
was too soft at the origin). In order to increase the suppression, an
additional factor $[1-e^{-(a_{3}r_{\alpha\alpha})^{2}}]^{2}$ was introduced.

The wave function (8) with the radial part (15) is an exact solution of Eq.
(2) for the $\alpha\alpha$ potential shown in Fig. (2) by the solid curve. As
is seen, this reconstructed potential is very close to the Ali-Bodmer one. On
one hand it is a bit different, but on the other hand it corresponds to the
exact experimental binding energy of ${}^{9}\mathrm{Be}$.

Figure 2: The exact Ali-Bodmer potential (5) shown by the dashed curve, and
the $\alpha\alpha$ potential corresponding to the wave function (8) with the
radial part (15) (solid curve).

In order to check how good the reconstructed potential is in describing the
low-energy $\alpha\alpha$ scattering, the $S$-wave scattering phase shifts
were calculated and compared with the corresponding values for the Ali-Bodmer
potential. For obtaining the phase shifts the differential equations for the
Jost functions were numerically integrated (see Eqs. (8.57, 8.58) of the book
[22]). The results of these calculations are shown in Fig. 3, where the solid
curve corresponds to the reconstructed potential while the dashed curve
represents the Ali-Bodmer one. In Fig. 3 these phase shifts are compared not
only with each other but also with the corresponding experimental values. The
experimental data are available (see Refs. [23, 24]) as the so-called nuclear
parts of the phase shifts, i.e. the pure Coulomb phase shifts are subtracted.
In calculating the curves such a subtraction was done as well.

Figure 3: The $S$-wave phase shift for the Ali-Bodmer potential (5) shown by
the dashed curve, and the corresponding phase shift for the $\alpha\alpha$
potential reconstructed from the wave function (8) with the radial part (15)
(solid curve). The experimental data are taken from Refs. [23, 24].

The purpose of this comparison is to indirectly check the quality of the wave
function (8) with the radial part (15). Since the reconstructed potential is
not very much different from the Ali-Bodmer one and since it generates almost
the same phase shifts that are also close to experimental data, one can say
that in a possible application of this wave function it can be considered as a
reliable approximation of the exact solution of the three-body Schrödinger
equation (2) with the two-body potentials given in Sec. 3.

An additional test of the wave function that is obtained in the present work
can be done by calculating various space-distances in the bound $(\alpha\alpha
n)$ state which this wave function describes. It should be emphasized that
none of the root-mean-square (RMS) distances that are calculated and given in
Table 1, were fitted. The fit (by varying the parameters $a_{1}$, $a_{2}$, and
$a_{3}$) was only done for the $\alpha\alpha$ potential. The RMS distances are
just those that the resulting wave function gives.

The calculation of the RMS distances is reduced to one- or three-dimensional
integrals over $r_{\alpha\alpha}$, $r_{n}$, and $\vartheta$ (see Fig. 1):

$\langle
r_{\alpha\alpha}^{2}\rangle=\left\langle\Psi\left|r_{\alpha\alpha}^{2}\right|\Psi\right\rangle\
,$ (17) $\langle
r_{\mathrm{cm}-n}^{2}\rangle=\left\langle\Psi\left|\left(\frac{2m_{\alpha}}{2m_{\alpha}+m_{n}}r_{n}\right)^{2}\right|\Psi\right\rangle\
,$ (18) $\langle
r_{\mathrm{cm}-\alpha}^{2}\rangle=\left\langle\Psi\left|\left(\frac{1}{2}\vec{r}_{\alpha\alpha}-\frac{m_{n}}{2m_{\alpha}+m_{n}}\vec{r}_{n}\right)^{2}\right|\Psi\right\rangle\
.$ (19)

When finding the RMS charge and matter radii of ${}^{9}\mathrm{Be}$, the
corresponding non-zero radii of the $\alpha$ particle,
$R_{\mathrm{ch}(\alpha)}=1.67824\,\mathrm{fm}$ [25] and
$R_{\mathrm{mat}(\alpha)}=1.457\,\mathrm{fm}$ [26] should be taken into
account:

$\langle r_{\mathrm{ch}}^{2}\rangle=\langle
r_{\mathrm{cm}-\alpha}^{2}\rangle+R_{\mathrm{ch}(\alpha)}^{2}\ ,$ (20)
$\langle
r_{\mathrm{mat}}^{2}\rangle=\frac{2m_{\alpha}}{2m_{\alpha}+m_{n}}\left[\langle
r_{\mathrm{cm}-\alpha}^{2}\rangle+R_{\mathrm{mat}(\alpha)}^{2}\right]+\frac{m_{n}}{2m_{\alpha}+m_{n}}\langle
r_{\mathrm{cm}-n}^{2}\rangle\ .$ (21)

distance | calculated (fm) | measured (fm)
---|---|---
$\sqrt{\langle r_{\alpha\alpha}^{2}\rangle}$ | 3.46 |
$\sqrt{\langle r_{\mathrm{cm}-n}^{2}\rangle}$ | 4.65 |
$\sqrt{\langle r_{\mathrm{cm}-\alpha}^{2}\rangle}$ | 1.83 |
$\sqrt{\langle r_{\mathrm{ch}}^{2}\rangle}$ | 2.48 |  2.519 [27, 28]
$\sqrt{\langle r_{\mathrm{mat}}^{2}\rangle}$ | 2.70 |  $2.50\pm 0.01$ [29]
Table 1: Various RMS distances in ${}^{9}\mathrm{Be}$: between the $\alpha$
particles ($\sqrt{\langle r_{\alpha\alpha}^{2}\rangle}$), from c.m. to the
valence neutron ($\sqrt{\langle r_{\mathrm{cm}-n}^{2}\rangle}$), from c.m. to
the $\alpha$ particle ($\sqrt{\langle r_{\mathrm{cm}-\alpha}^{2}\rangle}$),
charge radius of the nucleus ($\sqrt{\langle r_{\mathrm{ch}}^{2}\rangle}$),
matter radius of the nucleus ($\sqrt{\langle r_{\mathrm{mat}}^{2}\rangle}$).
The experimental data are taken from Refs. [27, 28, 29].

Figure 4: Radial part (15) of the wave function of ${}^{9}\mathrm{Be}$.

As is seen, the charge radius, $\sqrt{\langle r_{\mathrm{ch}}^{2}\rangle}$,
turned out to be very close to the corresponding experimental value. The
matter radius, $\sqrt{\langle r_{\mathrm{mat}}^{2}\rangle}$, is a bit too
large as compared to the measurement. The difference, however, is not drastic.
A very close result, namely, the value of 2.68 fm for the matter radius of
${}^{9}\mathrm{Be}$ was obtained in Ref. [30] within the three-body model with
microscopic nonlocal interactions based on Volkov V2 potential. Also, it
should be noted that in contrast to a rather simple and accurate procedure for
experimental determination of the charge radius, the measurements of the
matter radii of the nuclei are always more difficult and are not without some
ambiguities [26].

A quasi-3D image of the radial part (15) of the wave function of
${}^{9}\mathrm{Be}$ is given in Fig. 4. This function has a maximum at
$r_{\alpha\alpha}=2.4051\,\mathrm{fm}$ and $r_{n}=1.4654\,\mathrm{fm}$. It is
also clearly seen a “ridge” along $r_{n}$ with $r_{\alpha\alpha}\sim
2\,\mathrm{fm}$. This extended “ridge” is the cause of rather large values of
$\sqrt{\langle r_{\mathrm{cm}-n}^{2}\rangle}$ and $\sqrt{\langle
r_{\mathrm{mat}}^{2}\rangle}$.

##  5 Conclusion

The wave function (8) with the radial part (15) is an exact bound-state
solution of the three-body integro-differential equation (2) which follows
from the equivalent Schrödinger equation (2). By construction, this solution
corresponds to the experimentally known binding energy, as well as to the
$\alpha n$ potential (7), and to the $\alpha\alpha$ potential that is almost
the same as the Ali-Bodmer one. It can therefore be considered as a reliable
approximation of the wave function of the nucleus ${}^{9}\mathrm{Be}$ in the
three-body model ($\alpha\alpha n$).

This wave function is obtained in the form of a compact analytical expression
that is easy to use. Among possible applications of this function could be,
for example, the constructing of various folding potentials as well as
estimating the cross sections of various nuclear reactions that involve
${}^{9}\mathrm{Be}$.

The method used in the present paper, is universal and can be applied for
constructing the ground-state wave functions of any other nuclear or some
atomic bound three-body systems. Of course it cannot compete with the
approaches based on Faddeev equations, hyperspherical expansions, etc., where
whole spectrum of the excited states can be described. However, the proposed
method has its own application niche where it may be preferable. In
particular, this method is suitable for the systems that have no excited
states, such as, for example, the nucleui ${}^{6}\mathrm{He}$ and
${}^{6}\mathrm{Li}$ if they are considered in the three-body model $(\alpha
NN)$. The proposed method can give an analytical expression for the wave
function which is easy to use when one needs a reliable estimate of something
and does not want to delve into complicated numerical calculations of the
rigorous three-body theory. In principle, the excited state wave functions can
be constructed within this method as well. This, however, will require
additional effort in making the excited states orthogonal to the ground state
and to each other.

The method can be generalized for constructing the wave functions that are
linear combinations of several mutually orthogonal components. For example,
the ground $0^{+}$-state of ${}^{6}\mathrm{He}$ in the $(\alpha nn)$-model is
mainly composed of the two components: the state $|\Psi_{0}\rangle$ with zero
$nn$-spin as well as zero orbital angular momenta associated with both $nn$
and $\alpha(nn)$ Jacobi coordinates (82.87%), plus the state
$|\Psi_{1}\rangle$ with the $nn$-spin one and the $P$-waves along both Jacobi
coordinates (13.96%) [31], i.e.
$|\Psi\rangle=C_{0}|\Psi_{0}\rangle+C_{1}|\Psi_{1}\rangle$. In such a case,
when going from Eq. (2) (which should be written for the $nn$-potential) to
Eq. (2), one should multiply Eq. (2) from the left either by
$\Psi_{0}^{\dagger}$ or by $\Psi_{1}^{\dagger}$ and integrate over all the
configuration-space variables except the $nn$-distance. As a result one
obtains a system of two coupled integro-differential equations where on the
left hand sides are the singlet (${}^{1}S_{0}$) and triplet (${}^{3}P_{1}$)
$nn$-potentials which can be fitted by adjusting the parameters of some
postulated $\Psi_{0}$ and $\Psi_{1}$.

Another possible generalization of the proposed method may consist in
considering four-, five-, and larger numbers of bodies in the bound system. In
doing this, one term of the Schrödinger equation that involves a chosen two-
body potential, can be moved to the left-hand side as is done in Eq. (2). Then
whole equation should be multiplied from the left by $\Psi^{\dagger}$ and
integrated over all the variables except the one which the chosen potential
depends on. After that the wave function should be postulated and its
parameters be adjusted in order to fit the potential. The main difficulty here
is that with increasing number of particles it is not a simple task to find
(to guess) an appropriate parametrization of the wave function because it
depends on many variables. Moreover, even if one manages to find such a
parametrization, the right-hand side of Eq. (2) will involve multi-dimensional
integrals which may cause some numerical difficulties when the system under
consideration consists of too many particles.

## References

  * [1] Y.C. Tang, F.C. Khanna, R.C. Herndon, K. Wildermuth, “$\mathrm{Be}^{9}$ in the cluster model”, Nucl. Phys., 35, 421-433 (1962)
  * [2] V.D. Efros, H. Oberhummer, A. Pushkin, I.J. Thompson, “Low-energy photodisintegration of ${}^{9}\mathrm{Be}$ and $\alpha+\alpha+n\leftrightarrow{}^{9}\mathrm{Be}+\gamma$ reactions at astrophysical conditions”, The European Physical Journal A - Hadrons and Nuclei, 1(4), 447-453 (1998)
doi:10.1007/s100500050079

  * [3] I. Filikhin, V.M. Suslov, B. Vlahovic, “$\mathrm{Be}^{9}$ Low-Lying Spectrum Within a Three-Cluster Model, Few-Body Syst., 50, 255-257 (2011)
https://doi.org/10.1007/s00601-010-0135-3

  * [4] J. Casal, M. Rodrıiguez-Gallardo, J.M. Arias, and I.J. Thompson, “Astrophysical reaction rate for ${}^{9}\mathrm{Be}$ formation within a three-body approach”, Phys.Rev. C, 90, 044304 (2014)
DOI: 10.1103/PhysRevC.90.044304

  * [5] Elena Filandri, Paolo Andreatta, Carlo A. Manzata, Chen Ji, W. Leidemann, and G. Orlandini, “Beryllium-9 in Cluster Effective Field Theory”, arXiv:2002.00780v1 [nucl-th] (2020)
  * [6] D.R. Tilley, J.H. Kelley, J.L. Godwin, D.J. Millener, J.E. Purcell, C.G. Sheu, and H.R. Weller, “Energy levels of light nuclei A = 8, 9, 10”, Nucl. Phys., A 745, 155-362 (2004)
  * [7] Claus E. Rolfs and William S. Rodney, “Cauldrons in the Cosmos”, The University of Chicago Press (1988)
  * [8] A. Aprahamian, K. Langanke, and M. Wiescher, “Nuclear structure aspects in nuclear astrophysics”, Prog. Part. Nucl. Phys., 54, 535-613 (2005)
  * [9] Dr. Md A Khan, “Nuclear Astrophysics - A Course of Lectures”, CRC Press, Taylor & Francis Group, International Standard Book Number-13: 978-1-138-58816-5, (2018)
  * [10] Peter Mohr, “Direct capture cross section of ${}^{9}\mathrm{Be}(n,\gamma){}^{10}\mathrm{Be}$”, Phys. Rev. C, 99, 055807 (2019)
  * [11] D. Carbone, M. Bondi, A. Bonaccorso, C. Agodi, F. Cappuzzello, M. Cavallaro, R. J. Charity, A. Cunsolo, M. De Napoli, and A. Foti, “First application of the $n-{}^{9}\mathrm{Be}$ optical potential to the study of the ${}^{10}\mathrm{Be}$ continuum via the $({}^{18}\mathrm{O},{}^{17}\mathrm{O})$ neutron-transfer reaction”, Phys. Rev. C, 90, 064621 (2014)
  * [12] Su Jun, Li Zhi-Hong, Guo Bing, Liu Wei-Ping, Bai Xi-Xiang, Zeng Sheng, Lian Gang, Yan Sheng-Quan, Wang Bao-Xiang, and Wang You-Bao, “Astrophysical Reaction Rates of the ${}^{8}\mathrm{Li}(p,\gamma){}^{9}\mathrm{Be}_{\mathrm{g}.s.}$ Direct Capture Reaction”, Chinese Physics Letters, 23 (1), 55-57 (2006)
  * [13] S.B. Dubovichenko, N.A. Burkova, A.V. Dzhazairov-Kakhramanov, “The role of resonances in the capture of ${}^{8}\mathrm{Li}(p,\gamma){}^{9}\mathrm{Be}$ on the reaction rate of the relevant astrophysical synthesis of ${}^{9}\mathrm{Be}$”, Nucl. Phys. A, 1000,121842 (2020) https://doi.org/10.1016/j.nuclphysa.2020.121842
  * [14] N. Itagaki, and K. Hagino, “Low-energy photodisintegration of ${}^{9}\mathrm{Be}$ with the molecular orbit model”, Phys. Rev. C, 66, 057301 (2002)
  * [15] P. Descouvemont, T. Druet, L.F. Canto, M.S. Hussein, “Low-energy ${}^{9}\mathrm{Be}+{}^{208}\mathrm{Pb}$ scattering, breakup, and fusion within a four-body model”, Phys. Rev. C, 91, 024606 (2015) DOI: 10.1103/PhysRevC.91.024606
  * [16] P. Descouvemont, and N. Itagaki, “${}^{9}\mathrm{Be}$ scattering with microscopic wave functions and the continuum-discretized coupled-channel method”, Phys. Rev. C, 97, 014612 (2018) DOI: 10.1103/PhysRevC.97.014612
  * [17] V.B. Belyaev, S.A. Rakityansky, I.M. Gopane, “Recovering the Two-Body Potential from a Given Three-Body Wave Function”, Few-Body Syst., 64, 4 (2023) https://doi.org/10.1007/s00601-022-01785-7
  * [18] S. Ali and A.R. Bodmer, “PHENOMENOLOGICAL alpha-alpha POTENTIALS”, Nucl.Phys.,80, 99-112 (1966)
  * [19] J. Bang, C. Gignoux, “A REALISTIC THREE-BODY MODEL OF ${}^{6}\mathrm{Li}$ WITH LOCAL INTERACTIONS”, Nucl. Phys. A, 313, 119-140 (1979)
  * [20] I.J. Thompson, B.V. Danilin, V.D. Efros, J.S. Vaagen, J.M. Bang, M.V. Zhukov, “Pauli blocking in three-body models of halo nuclei”, Phys.Rev. C, 61 024318 (2000)
  * [21] F. James and M. Roos, “Minuit - a system for function minimization and analysis of the parameter errors and correlations”, Comp. Phys. Comm. 10 (1975) 343-367; http://hep.fi.infn.it/minuit.pdf
  * [22] S. A. Rakitiansky, “Jost Functions in Quantum Mechanics: A Unified Approach to Scattering, Bound, and Resonant State Problems”, Springer (2022)
  * [23] S.A. Afzal, A.A.Z. Ahmad, S. Ali, “Systematic Survey of the $\alpha$-$\alpha$ Interaction”, Rev. Mod. Phys., 41(1), 247-273 (1969)
  * [24] W.S. Chien and Ronald E. Brown, “Study of the $\alpha$+$\alpha$ system below 15 MeV (c.m.)”, Phys. Rev. C10(5), 1767-1784 (1974)
  * [25] Wilfried Nortershauser, “Helium nucleus measured with record precision”, Nature, 589, 518-519 (2021)
  * [26] I. Tanihata, H. Savajols, R. Kanungo, “Recent experimental progress in nuclear halo structure studies”, Progress in Particle and Nuclear Physics, 68,215-313 (2013)
  * [27] I. Angeli, K.P. Marinova, “Table of experimental nuclear ground state charge radii: An update”, At. Data and Nucl. Data Tables, 99(1), 69-95 (2013)
https://doi.org/10.1016/j.adt.2011.12.006

  * [28] Tao Li, Yani Luo, Ning Wang, “Compilation of recent nuclear ground state charge radius measurements and tests for models”, At. Data and Nucl. Data Tables, 140, 101440 (2021)
https://doi.org/10.1016/j.adt.2021.101440

  * [29] I. Tanihata, H. Hamagaki, O. Hashimoto, Y. Shida, and N. Yoshikawa, “Measurements of Interaction Cross Sections and Nuclear Radii in the Light $p$-Shell Region”, Phys.Rev.Lett., 55(24), 2676-2679 (1985)
  * [30] M. Theeten, H. Matsumura, M. Orabi, D. Baye, P. Descouvemont, Y. Fujiwara, and Y. Suzuki, “Three-body model of light nuclei with microscopic nonlocal interactions”, Phys. Rev., C76, 054003 (2007)
  * [31] B.V. Danilin and M.V. Zhukov, S.N. Ershov, F.A. Gareev, and R.S. Kurmanov, J.S. Vaagen, J.M. Bang, “Dynamical multicluster model for electroweak and charge-exchange reactions”, Phys. Rev., C43(6), pp. 2835-2843 (1991)