Production of radioactive isotopes through cosmic muon spallation in KamLAND
S. Abe,1 S. Enomoto,1, 2 K. Furuno,1 Y. Gando,1 H. Ikeda,1 K. Inoue,1, 2 Y. Kibe,1 Y. Kishimoto,1 M. Koga,1, 2 Y. Minekawa,1
T. Mitsui,1 K. Nakajima,1, ∗ K. Nakajima,1 K. Nakamura,1, 2 M. Nakamura,1 I. Shimizu,1 Y. Shimizu,1 J. Shirai,1 F. Suekane,1A. Suzuki,1 Y. Takemoto,1 K. Tamae,1 A. Terashima,1 H. Watanabe,1 E. Yonezawa,1 S. Yoshida,1 A. Kozlov,2H. Murayama,2, 3 J. Busenitz,4 T. Classen,4, † C. Grant,4 G. Keefer,4 D.S. Leonard,4, ‡ D. McKee,4, § A. Piepke,2, 4
T.I. Banks,3 T. Bloxham,3 J.A. Detwiler,3 S.J. Freedman,2, 3 B.K. Fujikawa,2, 3 F. Gray,3, ¶ E. Guardincerri,3 L. Hsu,3, ∗∗K. Ichimura,3 R. Kadel,3 C. Lendvai,3 K.-B. Luk,3 T. O’Donnell,3 H.M. Steiner,3 L.A. Winslow,3, †† D.A. Dwyer,5C. Jillings,5, ‡‡ C. Mauger,5, §§ R.D. McKeown,5 P. Vogel,5 C. Zhang,5 B.E. Berger,6 C.E. Lane,7 J. Maricic,7 T. Miletic,7
M. Batygov,8 J.G. Learned,8 S. Matsuno,8 S. Pakvasa,8 J. Foster,9 G.A. Horton-Smith,2, 9 A. Tang,9 S. Dazeley,10, ¶¶K.E. Downum,11 G. Gratta,11 K. Tolich,11, ∗∗∗ W. Bugg,12 Y. Efremenko,2, 12 Y. Kamyshkov,12 O. Perevozchikov,12
H.J. Karwowski,13 D.M. Markoff,13 W. Tornow,13 K.M. Heeger,2, 14 F. Piquemal,15 J.-S. Ricol,15 and M.P. Decowski2, 16(The KamLAND Collaboration)
1

Research Center for Neutrino Science, Tohoku University, Sendai 980-8578, Japan
Institute for the Physics and Mathematics of the Universe, Tokyo University, Kashiwa 277-8568, Japan3
Physics Department, University of California, Berkeley and
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
4
Department of Physics and Astronomy, University of Alabama, Tuscaloosa, Alabama 35487, USA
5
W. K. Kellogg Radiation Laboratory, California Institute of Technology, Pasadena, California 91125, USA
6
Department of Physics, Colorado State University, Fort Collins, Colorado 80523, USA
7
Physics Department, Drexel University, Philadelphia, Pennsylvania 19104, USA
8
Department of Physics and Astronomy, University of Hawaii at Manoa, Honolulu, Hawaii 96822, USA
9
Department of Physics, Kansas State University, Manhattan, Kansas 66506, USA
10
Department of Physics and Astronomy, Louisiana State University, Baton Rouge, Louisiana 70803, USA
11
Physics Department, Stanford University, Stanford, California 94305, USA
12
Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA
13
Triangle Universities Nuclear Laboratory, Durham, North Carolina 27708, USA andPhysics Departments at Duke University, North Carolina Central University, and the University of North Carolina at Chapel Hill
14
Department of Physics, University of Wisconsin, Madison, Wisconsin 53706, USA
15
CEN Bordeaux-Gradignan, IN2P3-CNRS and University Bordeaux I, F-33175 Gradignan Cedex, France
16
Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands(Dated: February 25, 2010)
2

Radioactive isotopes produced through cosmic muon spallation are a background for rare event detectionin ν detectors, double-beta-decay experiments, and dark-matter searches. Understanding the nature of cosmogenic backgrounds is particularly important for future experiments aiming to determine the pep and CNO
solar neutrino fluxes, for which the background is dominated by the spallation production of 11 C. Data fromthe Kamioka Liquid scintillator Anti-Neutrino Detector (KamLAND) provides valuable information for better understanding these backgrounds, especially in liquid scintillator, and for checking estimates from current
simulations based upon MUSIC, FLUKA, and Geant4. Using the time correlation between detected muonsand neutron captures, the neutron production yield in the KamLAND liquid scintillator is measured to be
Yn = (2.8 ± 0.3) × 10−4 µ−1 g−1 cm2 . For other isotopes, the production yield is determined from the observed time correlation related to known isotope lifetimes. We find some yields are inconsistent with extrapolations based on an accelerator muon beam experiment.
PACS numbers: 96.50.S-,96.60.-j,25.20.-x,25.30.Mr∗ Present

address: Center of Quantum Universe, Okayama University,
Okayama 700-8530, Japan† Present address: Department of Physics, University of California, Davis,
California 95616, USA‡ Present address: Department of Physics, University of Maryland, College
Park, Maryland 20742, USA§ Present address: Department of Physics, Kansas State University, Manhattan, Kansas 66506, USA¶ Present address: Department of Physics and Computational Science, Regis
University, Denver, Colorado 80221, USA∗∗ Present address: Fermi National Accelerator Laboratory, Batavia, Illinois

60510, USA†† Present address: Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA‡‡ Present address: SNOLAB, Lively, ON P3Y 1M3, Canada§§ Present address: Subatomic Physics Group, Los Alamos National Laboratory, Los Alamos, NM 87545, USA¶¶ Present address: Lawrence Livermore National Laboratory, Livermore,
California 94550, USA∗∗∗ Present address: CENPA, University of Washington, Seattle, Washington
98195, USA2
I.

INTRODUCTION

Cosmic ray muons and their spallation products are potential sources of background for neutrino detectors, double-betadecay experiments, and dark-matter searches, even when thedetectors are deployed underground. Characterizing cosmicray-muon-induced backgrounds, particularly the secondary
neutrons and radioactive isotopes produced by muon-initiatedspallation processes, is essential for interpreting these experiments.
Liquid-scintillator detectors such as KamLAND, Borexino [1–3], CANDLES IV [4, 5], SNO+ [6, 7], LENS [8],and LENA [9] are designed to detect low-energy phenomena. In organic liquid scintillator (LS), energetic muons andsubsequent showers interact mostly with 12 C, the most abundant nucleus heavier than 1 H in the LS, generating neutrons
and isotopes by electromagnetic or hadronic processes. The
muon-initiated spallation of carbon targets is a matter of primary interest.
Isotope production by muon-initiated spallation has beenstudied by an earlier experiment [10] using the CERN Super Proton Synchrotron (SPS) muon beam. The energy dependence was studied with 100 and 190 GeV incident muons.
The production yield at other energies is estimated from thisdata by extrapolation, assuming a power-law dependence on
the muon energy. Direct measurements of the production
yield by underground detectors such as LSD [11], LVD [12],and Borexino [13] were compared to calculations exploiting
simulations based on MUSIC [14], FLUKA [15, 16], andGeant4 [17, 18]. Particular attention was paid to neutron
production since isotope production measurements are difficult with the small scintillator masses used in these detectors.
KamLAND, owing to its larger mass – ∼1 kton of LS – does
not suffer from this difficulty and is well placed to study avariety of isotopes of interest.
This paper presents the neutron and isotope production
rates in KamLAND from muon-initiated spallation based
upon data collected from 5 March 2002 to 12 May 2007. The
results are compared to simulations and other experiments.
These comparisons provide important information for validating Monte Carlo simulations.

II.

DETECTOR DESCRIPTION AND PERFORMANCE

KamLAND is located under the peak of Ikenoyama (Ike
Mountain, 36.42◦ N, 137.31◦ E), and the vertical rock overburden is approximately 2700 meters water equivalent (m.w.e.).
A schematic diagram of KamLAND is shown in Fig. 1.
KamLAND consists of an active detector region of approximately 1 kton of ultra-pure LS contained in a 13-mdiameter spherical balloon made of 135 µm thick transparent
nylon/EVOH (ethylene vinyl alcohol copolymer) composite
film and supported by a network of Kevlar ropes. In addition to providing containment for the LS, the balloon protects
the LS against the diffusion of ambient radon from the surrounding components. The total volume of LS in the balloon is 1171 ± 25 m3 as determined with flow meters during

FIG. 1: Schematic diagram of the KamLAND detector at the
Kamioka underground laboratory.

TABLE I: Elemental composition of the KamLAND liquid scintillator. The hydrogen-to-carbon ratio was verified by elemental analysis (±2%). The liquid scintillator contains traces of nitrogen andoxygen from its fluor, PPO (C15 H11 NO), and dissolved gases. Measurements in similar liquid hydrocarbons [19, 20] indicate that theamount of dissolved nitrogen saturates at about 200 to 500 parts per
million (ppm). The range of nitrogen composition values given inthis table reflects the extreme cases of zero dissolved nitrogen gas
to full saturation. The dissolved oxygen content of liquid scintillator
taken from the center of KamLAND was measured to be less than
3 ppm, which is insignificant compared to the oxygen contribution
from PPO.
Number of Targets
Element

Stoichiometry

HydrogenCarbonNitrogen
Oxygen1.97≡1
1 × 10−4 to 6 × 10−4
1 × 10−4(per kiloton)
8.47 × 1031
4.30 × 1031
5 × 1027 to 3 × 1028
5 × 1027

the initial filling. The LS comprises of 80% dodecane, 20%
pseudocumene (1,2,4-Trimethylbenzene) by volume, and1.36 ± 0.03 g/liter of the fluor PPO (2,5-Diphenyloxazole).
The density of the LS is 0.780 g/cm3 at 11.5◦ C. The calculated elemental composition of the LS is given in Table I.
A buffer comprising of 57% isoparaffin and 43% dodecaneoils by volume fills the region between the balloon and the surrounding 18-m-diameter spherical stainless-steel outer vesselto shield the LS from external radiation. The specific gravity of the buffer oil (BO) is adjusted to be 0.04% lower thanthat of the LS. An array of photomultiplier tubes (PMTs),
1325 specially developed fast PMTs masked to 17-inch diameter and 554 older 20-inch diameter PMTs reused from the
Kamiokande experiment [21], are mounted on the inner surface of the outer containment vessel, providing 34% photocathode coverage. During the period from 5 March 2002 to
27 February 2003 the photo-cathode coverage was only 22%,
since the 20-inch PMTs were not operated. A 3 mm thickacrylic barrier at 16.6-m-diameter helps prevent radon emanating from the PMT glass from entering the BO. The innerdetector (ID), consisting of the LS and BO regions, is surrounded by a 3.2 kton water-Čerenkov detector instrumented
with 225 20-inch PMTs. This outer detector (OD) absorbs
γ-rays and neutrons from the surrounding rock and enables
tagging of cosmic-ray muons.
The KamLAND front-end electronics (FEE) system
is based on the Analog Transient Waveform Digitizer(ATWD) [22] which captures PMT signals in 128 10-bit digital samples at intervals of 1.5 ns. Each ATWD captures threegain levels of a PMT signal to obtain a dynamic range fromone photoelectron (p.e.) to 1000 p.e. Each ATWD takes 27 µs
to read out, so two are attached to each PMT channel to reduce dead time. The FEE system contains discriminators set
at 0.15 p.e. (∼0.3 mV) threshold which send a 125 ns long
logic signal to the trigger electronics. The trigger electronics
counts the number of ID and OD PMTs above the discriminator threshold with a sampling rate of 40 MHz and initiates
readout when the number of 17-inch ID PMTs above the discriminator threshold (N17 ) exceeds the number corresponding
to ∼0.8 MeV deposited energy. The trigger system also issues independent readout commands when the number of OD
PMTs above threshold exceeds a preset number.
The energy can be estimated from the Nmax parameter, defined as the maximum value of N17 in a 200 ns period following the trigger command. However, the offline analysis takes
full advantage of the information stored in the digitized PMT
signals by identifying individual PMT pulses in the waveform
information that is read out. The time and integrated area
(called charge) are computed from the individual pulses. For
each PMT, the average charge corresponding to a single p.e.
is determined from single-pulse waveforms observed in low
occupancy events. The ID PMT timing is calibrated with lightpulses from a dye laser (∼1.2 ns pulse width), injected at thecenter of the detector through an optical fiber. The verticesof spatially localized low-energy (<30 MeV) events are estimated by comparing calculated time-of-flights of optical photons from the hypothetical vertex to the measured arrival times
at the PMTs in KamLAND.
The reconstructed energies of events are calibrated with γ
sources: 203 Hg, 68 Ge, 65 Zn, and 60 Co; and with n+γ sources:241
Am+9 Be and 210 Po+13 C [23]. These are deployed at various positions along the vertical axis of the detector and occasionally off the vertical axis within 5.5 m from the detector
center [24]. Such calibrations cover energies between 0.28
and 6.1 MeV. The energy calibration is aided with studies
of background contaminants 40 K and 208 Tl, 212 Bi-212 Po and214
Bi-214 Po sequential decays, 12 B and 12 N spallation products, and γ’s from thermal neutron captures on 1 H and 12 C.
The visible energy (Evis ) of an event is computed from themeasured light yield. Specifically, Evis is the number of detected p.e. after corrections for PMT variation, dark noise,solid angle, shadowing by suspension ropes, optical transparencies, and scattering properties in the LS. The relationship between Evis and the deposited energy (Edep ) of γ’s,
e± ’s, protons, and α’s is non-linear and modeled as a combination of Birks-quenched scintillation [25, 26] and Čerenkov

Latitude (deg)3

Altitude (m)
1600

36.5

1400
1200
36.45
1000
800
36.4

600
400
200

36.35
0
137.25137.3

137.35Longitude (deg)FIG. 2: Ikenoyama topological profile [27]. The black point near thecenter is the location of KamLAND.

radiation. The scale is adjusted so that Evis is equal to Edep
1for the 2.225 MeV γ-ray from neutron
p capture on H. The observed energy resolution is ∼7.4%/ Evis (MeV)
p for the period without the 20-inch PMTs, and ∼6.5%/ Evis (MeV)
for the rest of the data.
The calibration sources are also used to determine systematic deviations in position reconstruction by comparison withthe source’s known position. This comparison gives an average position reconstruction uncertainty of less than 3 cm forevents with energies in the range 0.28 to 6.1 MeV.

III.

COSMIC RAY MUONS

A digital map [27] of the topological profile of Ikenoyama
is shown in Fig. 2. The vertical overburden at KamLAND isapproximately 1000 meters of rock and the minimum overburden corresponding to a nearby valley is approximately 900
meters.
Cosmic ray muons are identified either by the large amountof scintillation and Čerenkov light detected by the ID PMTs,
or by the Čerenkov light detected by the OD PMTs. Muons

4
× 106
Light Yield (p.e.)Light Yield (p.e.)

107
106
105

(a)

1

0.5

0
0

0

200
400
600
Impact Parameter (cm)800

FIG. 3: Correlation between the light yield in the inner detector andthe shortest distance from the muon track to the KamLAND detector
center (impact parameter). The vertical dashed red line represents
the boundary between the LS and BO at 650 cm. Events located to
the left of this boundary and above the horizontal blue dashed line
(LID = 4 × 104 p.e.) are designated as LS muons.

crossing the ID (ID muons) are selected by requiring that oneof the following conditions was satisfied:
1. LID ≥ 10000 p.e. (∼30 MeV)
2. LID ≥ 500 p.e. and NOD ≥ 5
where LID is the total light yield measured by the ID 17-inch
PMTs, and NOD is the number of OD PMTs with signals
above threshold. Approximately 93% of the ID muons satisfy the first selection criterion. The ID muon track is reconstructed from arrival times of the first-arriving Čerenkov
or scintillation photons at the PMTs. Since for relativisticmuons the wavefront of the scintillation light proceeds at theČerenkov angle, and since muons generate enough light togenerate photoelectrons in every PMT, by restricting the fit to
the first-arriving photons both Cerenkov and scintillation photons can be treated identically. The observed muon track isthen established by minimizing time-of-flight deviations fromhypothetical muon tracks. The fit converges for 97% of all
ID muon events. The majority of the events that are not reconstructed are believed to be multiple muons, or muons accompanied by large electromagnetic or hadronic showers forwhich the tracking model is not valid.
Muons passing only through the BO produce mostlyČerenkov light, whereas muons passing through the LS generate both Čerenkov and scintillation light. Figure 3 shows
the correlation between the light yield (LID ) and the shortest distance between the reconstructed muon track and thecenter of KamLAND (impact parameter). The boundary at
650 cm between the BO and LS regions is evident. The correlations between LID and the reconstructed muon track length
in the BO and LS regions (LBO and LLS , respectively) areplotted in Figs. 4a and 4b. A linear trend, corresponding tominimum ionizing muons, is apparent in both distributions.
The slope of each line is the light yield per unit length in
the respective material. In BO, where the light is predom-

Light Yield (p.e.)

104

6

200

× 104

400 600 800 1000 1200 1400
Track Length in LS (cm)
(b)

4
2
0
0

200

400 600 800 1000 1200 1400
Track Length in BO (cm)FIG. 4: Correlation between the total light yield measured by the
ID 17-inch PMTs (LID ) and the muon track length for (a) muonsthat pass through both the LS and BO regions and (b) muons thatpass through only the BO region. The red lines shows the fitted
light yield per unit length of hdLS /dXi = 629 ± 47 p.e./cm andhdLČ /dXi = 31 ± 2 p.e./cm.

inantly Cerenkov, the light yield per unit length is found tobe hdLČ /dXi = 31 ± 2 p.e./cm; the fit was restricted to path
lengths above 700 cm since fits at shorter path lengths arecomplicated by the presence of PMTs which may obstructsome of the emitted light. In the LS we obtain
LID − LBO hdLČ /dXi
dLS=
= 629 ± 47 p.e./cm, (1)dX
LLS
where dLS /dX includes the Čerenkov light created in the LS.
The muons in Fig. 4 generating light yields above the baseline
linear trend are likely to involve secondary particles. We define an excess light yield parameter ∆L,




dLČ
dLS
∆L = LID − LBO
− LLS
(2)
dX
dX
for the purpose of describing showering muons associatedwith secondary particles.
The LS muon rate is estimated by selecting muons with
LID > 4 × 104 p.e. and impact parameter < 650 cm. The
light yield cut has a negligible inefficiency for LS muons,
while the impact parameter cut eliminates LS muons thatreconstruct outside the balloon because of the resolution of
the fitter algorithm. This provides the lower limit on the
LS muon rate. An upper limit is established by removing the impact parameter cut and increasing the LID cut to> 105 p.e. to eliminate muons that pass through the BO without transversing the LS. This cut again has a small inefficiency for LS muons but does not eliminate muons which

5

350

TABLE II: Calculated Jµ and E µ values for various rock types using the simulation method from Ref. [33]. The range in values
for Jµ and E µ corresponds to varying the specific gravity of therock from 2.65 g/cm3 to 2.75 g/cm3 . The Generic Skarn is defined
to be 70% Granite and 30% Calcite (by weight) with the following properties: <Z>= 10.22, <A>= 20.55 u, and radiation length
X0 = 25.411 g/cm2 .

Events/10 cm

300
250
200
150
100
50
0
0

500
1000
Track Length in LS (cm)FIG. 5: Distribution of muon track lengths in the KamLAND LS. The
dashed blue line shows the expected linear distribution for a sphere of
nominal volume, and the average track length in this case is 867 cm.
The average muon track length measured in the detector is 878 cm,
in agreement with the value 874 ± 13 cm calculated using the proper
LS shape and muon angular distributions.

shower in the BO. Averaging the two measurements yields
Rµ = 0.198 ± 0.014 Hz, where the error estimate is dominated by systematic uncertainties due to temporal variations,the difference between the two measurements, and the presence of muon bundles, in which multiple simultaneous muons
may be tagged as single muon events. Using the parameterization of Ref. [28], the effect of muon bundles was estimatedto be <5%. The muon rate corresponds to an integrated muon
intensity of Jµ = 5.37 ± 0.41 m−2 h−1 , where the error includes uncertainties in the shape of the LS volume but is dominated by the uncertainty in the muon rate. Using a cut of∆L > 106 p.e., which is equivalent to a ∼3 GeV threshold,the rate of showering muons in the LS is ∼0.03 Hz. It is possible that some atmospheric neutrino interactions leak into theLS muon sample. However, an estimate with the neutrino
flux from Ref. [29] and cross sections from Ref. [30] gives
less than 4 × 10−5 Hz from atmospheric neutrinos. The muontrack length distribution is shown in Fig. 5. The measured average track length is Lµ = 878 cm, in agreement with the calculated value of Lµ = 874 ± 13 cm where the non-sphericalcorrections to the balloon shape and the muon angular distributions were taken into account.
The production yield of radioactive isotopes from muoninitiated spallation is energy dependent. KamLAND does notmeasure the muon energy, so it is estimated from simulation.
There are previous estimates of the mean energy (E µ ) ranging
from 198 GeV [31] to 285 GeV [10, 32] at the Kamiokande,
KamLAND, and SuperKamiokande detectors. The authors ofRef. [33] made a detailed determination of Jµ and E µ using
the MUon SImulation Code (MUSIC) [14] to transport muons,generated according to the Modified Gaisser Parameterization [34] sea-level muon flux distribution, through a digitalprofile of Ikenoyama. Although this calculation reproduces
the zenith and azimuthal angular distributions observed by
KamLAND, as shown in Fig. 10 of Ref. [33], it overestimatesJµ by 14% relative to this work. The calculation assumed a

Ikenoyama Rock Model
Inishi RockStandard RockGeneric Skarn
This MeasurementJµ (m−2 h−1 )5.66 to 6.71
4.95 to 5.83
4.90 to 5.82
5.37 ± 0.41

E µ (GeV)
262 to 268
256 to 262
254 to 260—

homogeneous rock entirely of the Inishi type (see Table II of
Ref. [33] for the chemical composition), but other rock types,
such as granite, limestone, and several types of metamorphic
rock, are common in Ikenoyama in unknown quantities [35].
Table II shows the result of calculations of Jµ and E µ using the same simulation method of Ref. [33] but for Standard
Rock [36, 37] and Generic Skarn, a generic mixture of rocktypes found in a skarn-type mine like that of Ikenoyama. Here,
Generic Skarn is defined to be 70% Granite and 30% Calcite by weight. The Jµ values range from 4.90 m−2 h−1 , for2.75 g/cm3 specific gravity Generic Skarn, to 6.71 m−2 h−1 ,for 2.65 g/cm3 specific gravity Inishi Rock, while E µ varies
from 254 GeV to 268 GeV. The value of Jµ for 2.70 g/cm3
specific gravity Standard Rock is 5.38 m−2 h−1 , in excellent
agreement with our measured value. The value of E µ for thisrock is 259 GeV. We take E µ = 260 ± 8 GeV, where the uncertainty is chosen to cover the full range for the various rocktypes.

IV.

SPALLATION NEUTRON YIELD

Most of the neutrons produced in the KamLAND LS capture on hydrogen or carbon atoms. The capture cross section
varies inversely with respect to velocity, and the mean neutron
capture time (τn ) is constant with respect to energy. The capture time (t) distribution is exponential, P (t) ∝ e−t/τn . A calculation using the elemental composition of the KamLAND
LS given in Table I and the thermal neutron capture cross sections from Ref. [38] gives τn = 206 µs. This calculation indicates that 99.5% of the neutrons capture on 1 H, while theremainder capture mostly on 12 C. The probability for captureon the other isotopes in the LS, such as 13 C, is 2 × 10−4 or
less.
Neutrons produced by muon-initiated spallation in the LS
can be identified by the characteristic capture γ-rays. Figure 6
shows the Evis distributions in signal (150 ≤ ∆T < 1000 µs)
and background (4150 ≤ ∆T < 5000 µs) coincidence windows following some muons, where ∆T ≡ (t − tµ ) is thetime elapsed since the muon’s passage. The Evis distribution clearly shows peaks from neutron captures on 1 H(2.225 MeV) and 12 C (4.9 MeV) which are not evident in the

6

n1H

105

104
103

Events/10 µs

Events/50 keV

105

n12C

2

10

104

10
103

1
1

2

3
4
5
6
Visible Energy (MeV)

7

8

FIG. 6: Visible energy spectra of the events following muons inthe signal time window 150 ≤ ∆T < 1000 µs (cross) comparedwith background in the time window 4150 ≤ ∆T < 5000 µs (dash).
The low-energy tails in the n1 H and n12 C peaks are caused by
the electronics effect discussed in the text. In the time window
1300 ≤ ∆T < 2150 µs (solid) there is a small n1 H peak without
the low-energy tail.

background. Figure 7 shows the ∆T distribution for events
within the LS volume and with 1.8 < Evis < 2.6 MeV, which
includes the single 2.225 MeV γ-ray emitted by neutron capture on 1 H. For ∆T < 1000 µs, there is a clear deviationfrom the exponential distribution due to the overload that large
muon signals produce on individual electronics channels andto the dead time in the system arising from the very high event
multiplicity following the muon. Both effects intervene in
events that are quite different to those which KamLAND wasdesigned to record.
The number of neutrons produced by muon-initiated spallation in the LS is established by a binned maximum likelihood
fit [39] to the data in Fig. 7 using the functionr(t) =Nn −(t−tµ )/τne
+ rB ,
τn(3)

where Nn is the total number of neutron captures associated with the selected events, τn is the mean neutron capture time, and rB is the background rate, which is assumed to be approximately constant in the region of interest (∆T < 2500 µs) due to the low muon rate (∼ 0.2 Hz).
To avoid the electronics-induced distortions, the fit is restricted to the region ∆T ≥ 1300 µs. The parameters Nn
and rB are free, but the mean capture time is constrained toτn = 207.5 ± 2.8 µs with a gaussian penalty function. This
mean capture time is determined from two independent measurements of τn : a 241 Am+9 Be calibration source, and ananalysis of a sample of neutrons generated by clean muons.
These clean muons are identified by a multiplicity parameterηT defined to be the number of trigger commands that follow the muon within a 10 ms period. Fits to the ∆T distribution of the subset of neutrons selected with various limits

500

1000

1500
∆T (µs)

2000

2500

FIG. 7: Time difference between a muon and the neutron capture
event, with the cuts 1.8 < Evis < 2.6 MeV. The red line shows
the fit restricted to the region ∆T ≥ 1300 µs. The fit results in
Nn = (4.2 ± 0.3) × 106 with χ2 /d.o.f.= 98/118.

TABLE III: Summary of the dominant contributions to the neutron
detection efficiency.
EffectNeutron Eliminating Reactions (e.g. (n, p))
Neutron Captures on 1 H
LS-BO Boundary
Electronics Dead Time Effects
Combined Efficiency

Value
96.3 ± 3.7%
99.5 ± 0.1%
93.3 ± 2.0%
> 98%
89.4 ± 3.8%

on ηT demonstrate that muons with ηT < 10 give unbiased
fit residuals. The value of τn is 207.5 ± 0.3 µs from the cleanmuon sample and 205.2 ± 0.5 µs from the 241 Am+9 Be source
data. The observed 2.3 µs discrepancy between these values
is not completely understood, but is suspected to be causedby neutrons from the 241 Am+9 Be source that capture on thestainless-steel source capsule. In this analysis, we use thevalue τn = 207.5 µs from the non-showering muon sample
with an uncertainty of ±2.8 µs covering both measurements.
The fit shown in Fig. 7 in the region ∆T ≥ 1300 µs results in
Nn = (4.2 ± 0.3) × 106 and χ2 /d.o.f.= 98/118.
The actual number of neutrons Nn produced by muoninitiated spallation is related to the fit result Nn by an efficiency n ,
Nn =

Nn
,
n

(4)

that accounts for other neutron-eliminating nuclear reactions
like 12 C(n, p) and a net neutron loss at the LS-BO boundary.
This efficiency is calculated using the MUSIC-based muon
simulation described in Sec. III and the Geant4-based Monte
Carlo of KamLAND described in Sec. VI A. The muon simulation generates muons with a 3-momentum distribution appropriate to the KamLAND site inside Ikenoyama. Muon

7

Yn =

Nn
,
Rµ TL ρLµ(5)

where Rµ = 0.198 ± 0.014 Hz is the measured rate of
LS muons, TL = 1.24 × 108 s is the detector live time,ρ = 0.780 ± 0.001 g/cm3 is the density of the LS, andLµ = 874 ± 13 cm is the calculated mean muon track
length.
The resulting neutron production yield is
Yn = (2.8 ± 0.3) × 10−4 µ−1 g−1 cm2 , with 64 ± 5% of theneutrons produced by events classified as showering muons.

Events/ms

5000
4000
Events/0.5 MeV

transport is then modeled through KamLAND. In this simulation neutrons are created and destroyed; neutrons that survive to thermalization are tracked until they are captured. Energy scale nonlinearities and the finite γ-ray resolution are included. Two competing effects of the LS-BO boundary aretaken into account. The first is for neutrons produced bymuons in the LS that leak out, leading to an undercountingof neutrons. The second is for neutrons produced outside of
the LS that leak in, primarily from the BO, leading to an overcount. The leak-out fraction (11.7 ± 1.9%) is larger than theleak-in fraction (5.4 ± 0.4%) resulting in a net 6.7 ± 2.0%
correction. The component of the efficiency related to theelectronics dead time for high-multiplicity events is measured
by comparing the number of recorded waveforms with thenumber of trigger commands. The correction is measured
to be less than 2% for ∆T ≥ 1300 µs. In counting neutrons,
we use the definition by which the (n, 2n) reaction generates
one new neutron and the (n, n0 ) reaction generates no newneutrons. This method gives an efficiency n = 89.4 ± 3.8%
(broken down in Table III).
Using Eqn. 4 we find Nn = (4.7 ± 0.4) × 106 , which isthen used to extract the neutron production yield:

3000

103
102
200
400
∆T (ms)

2000

Total
12
B

1000

12
9

0
4

6

8
10 12 14 16
Visible Energy (MeV)N
Li, 8He, 9C

18

20

FIG. 8: Background-subtracted Evis spectrum above 4 MeV.
Signal and background events are taken from 2 ≤ ∆T < 60 ms
and 502 ≤ ∆T < 560 ms, respectively.
The production rateof 12 B (τ = 29.1 ms, Q = 13.4 MeV) is estimated to be
54.8 ± 1.5 kton−1 day−1 , from fitting Eqn. 6 to the ∆T distribution shown in the inset. This fit has a χ2 /d.o.f.= 509/495. The
production rate 12 N (τ = 15.9 ms, Q = 17.3 MeV) is estimated tobe 2.2 ± 0.5 kton−1 day−1 , from a similar fit to the ∆T distribution
for events with 14 ≤ Evis < 20 MeV, where the higher Evisthreshold is imposed to exclude 12 B.

where Rµ , TL , ρ, and Lµ are defined as in Eqn. 5. The spallation production rate for isotope i is equal to
Ri =

Ni
,
ρVT TL(8)

where VT = 1171 ± 25 m3 is the target volume.
V.

SPALLATION ISOTOPE YIELDThe method for determining the yields of spallationgenerated isotopes is similar to the neutron analysis describedin Sec. IV. Spallation-generated isotopes are identified bytheir decay time relative to their creation, and by their decayenergy. The decay time, ∆T ≡ t − tµ , is calculated for each
event relative to all previous muons. Usually, several different
isotopes decay in a given time window. The number of eachisotope produced, Ni , is obtained from a binned maximum
likelihood fit [39] to the ∆T distribution, using the function:r(t) =

X Nii

τie−(t−tµ )/τi + rB ,(6)

where Ni is related to Ni by an event selection efficiency(Ni = Ni /i ), τi is the mean lifetime, and rB is a constantrate of uncorrelated background events that are in accidental
coincidence with muons.
In terms of Ni , the spallation production yield for isotope i
is equal to
Yi =

Ni
,
Rµ TL ρLµ

(7)

A.

12

B and 12 N

12B (τ = 29.1 ms, Q = 13.4 MeV) [40] β − decay andN (τ = 15.9 ms, Q = 17.3 MeV) [40] β + decay candidate events are selected via cuts on Evis and ∆T . The inset in Fig. 8 shows the distribution of ∆T for events with
4 ≤ Evis < 20 MeV and 2 ≤ ∆T < 500 ms. A binned maximum likelihood fit to the ∆T distribution using Eqn. 6, and including the long-lived isotopes 8 He, 9 Li, 9 C, and 12 N as contaminants, yields N (12 B) = (5.94 ± 0.10) × 104 events with
χ2 /d.o.f.= 509/495. Longer lived isotopes give roughly constant decay rates on this time scale and fit out as a componentof rB . A similar fit to the ∆T distribution for events with
14 ≤ Evis < 20 MeV, where the higher Evis threshold is
imposed to exclude 12 B, gives N (12 N) = (2.8 ± 0.3) × 102
events. A comparison with the predicted Evis spectrum isshown in Fig. 8. The Evis spectra are predicted from theallowed 12 B and 12 N β ± decay spectra, taking into accountthe KamLAND detector response, normalized to the observed
N (12 B) and N (12 N). The detector response model includes
the energy nonlinearities and boundary effects described inSec. II.
12

8

3000
Events/1.0 MeV

2500

10

× 103

103

9

Events/10 cm

Events/0.1 s

3500

8

2000

2

1500

4 6 8 10∆T (s)

Total
8

1000

Li
B

500

10-10
6

8
10 12 14 16
Visible Energy (MeV)18

20

FIG. 9: Background-subtracted Evis spectrum above 4 MeV, where
signal and background events are taken from 0.6 ≤ ∆T < 4.0 s
and 10.6 ≤ ∆T < 14.0 s, respectively.
The production rates
of 8 Li (τ = 1.21 s, Q = 16.0 MeV) and 8 B (τ = 1.11,s,Q = 18.0 MeV) are estimated to be 15.6 ± 3.2 kton−1 day−1
and 10.7 ± 2.9 kton−1 day−1 , respectively, from simultaneously
fitting the Evis spectrum and the ∆T distribution shown in the inset.
The fit to the ∆T distribution has a χ2 /d.o.f.= 95/91.

The inefficiency in identifying 12 B and 12 N candidates isdominated by the Evis cut. The efficiencies calculated by integrating the predicted Evis spectra over the selection window
give (12 B) = 82.9 ± 0.7% and (12 N) = 9.3 ± 1.6%, where
the errors come from the uncertainty in the detector response.
Using Eqn. 7, the resultant isotope production yields arecalculated to be Y (12 B) = (42.9 ± 3.3) × 10−7 µ−1 g−1 cm2and Y (12 N) = (1.8 ± 0.4) × 10−7 µ−1 g−1 cm2 . The production rates are R(12 B) = 54.8 ± 1.5 kton−1 day−1 andR(12 N) = 2.2 ± 0.5 kton−1 day−1 .

B.

8Li and 8 B

Li (τ = 1.21 s, Q = 16.0 MeV) [41] β − decay and 8 B(τ = 1.11 s, Q = 18.0 MeV) [41] β + decay candidate events
are selected according to Evis and ∆T . The inset in
Fig. 9 shows the distribution of ∆T for all preceding muonsfor events with 4 ≤ Evis < 20 MeV and 0.6 ≤ ∆T < 10 s,
where the lower limit on ∆T is chosen to exclude 9 C,
8
He, 9 Li, and other isotopes with shorter lifetimes. Isotopes with longer lifetimes are roughly constant over theselected range of ∆T . To avoid a large accidental coincidence background from uncorrelated muons, all showering muons and a portion of non-showering muons whose
track is within 3 m of a 8 Li/8 B candidate are selected. Figure 9 shows the Evis distribution of 8 Li/8 B candidate eventswith 0.6 ≤ ∆T < 4.0 s after subtracting background estimated from the range 10.6 ≤ ∆T < 14.0 s. N (8 Li) andN (8 B) are determined from a simultaneous binned maximum
likelihood fit to the ∆T distribution and a chi-square fit to the
8

10
1

8

0
4

102

200 400 600 800 1000 1200 1400
Distance from Muon Track (cm)FIG. 10: Impact parameter distribution for events identified as 12 B
(blue points) and neutrons (black histogram) produced by muoninduced spallation. Only non-showering muons (∆L < 106 p.e.)were used to make this plot. The efficiencies for the ∆L < 3 m
track cut (represented by the vertical dashed line) are evaluated to be
91.6 ± 4.3% from 12 B events, and 95.9 ± 0.7% from neutron events.

Evis distribution. The expected Evis spectra for KamLAND
are calculated by convolving the β ± spectra from Ref. [42, 43]with KamLAND’s detector response. For the fit to the Evis
distribution, the energy scale parameters are constrained to anallowed region determined by a prior fit to γ-ray calibration
data (described in Sec. II) and the 12 B Evis distribution in
Fig. 8.
The efficiency is given by the product of E for the Evis selection and S for the muon-8 Li/8 B spatial correlation. By
integrating the expected spectra over 4 ≤ Evis < 20 MeV,E (8 Li) and E (8 B) are estimated to be 77.6 ± 0.9% and88.4 ± 0.7%, respectively. For showering muons, S is 100%
because no correlation requirement is imposed. For nonshowering muons, S is estimated from the 12 B/12 N analysis (Sec. V A), and the systematic error arising from to variations between isotopes is estimated with a FLUKA simulation. Figure 10 shows the impact parameter (∆L) distribution
for the 12 B/12 N candidate events for non-showering muons(∆L < 106 p.e.). We find that 91.6 ± 4.3% of the candidates
are within 3 m of the muon track. This fraction is the value ofS (12 B) for ∆L < 3 m.
In order to obtain S (8 Li) and S (8 B), an additional correction for the difference between the muon-8 Li/8 B and themuon-12 B spatial correlations is applied. This correction isderived from the FLUKA simulation described in Sec. VI B.
The simulation does not include the uncertainties in the muontrack and the isotope decay vertex reconstruction; it is onlyused to study the isotope dependance of S . The range ofvalues of S from FLUKA for different spallation isotopes isused to estimate the systematic error that should be added toS (12 B) in order to obtain a common S for all spallation isotopes. The resulting value, S = 91.6 ± 8.4%, is used for estimating the 8 Li and 8 B (and, later the 9 C, 8 He, and 9 Li) yields.
Combiningtheabove
analyses,
weobtain

5
00

250

10
5
00

5
10 15
8
Li kton-1day-1

Events/0.1 s

10

300

(b)

15

5
10 15
Li kton-1day-1

8

Events/1.0 MeV(a)

B kton-1day-1

15

8

8

B kton-1day-1

9

200

103
102
0

150

2

4 6 8 10∆T (s)

Total
9
Li

100

8

He

50
8

FIG. 11: Allowed regions for the production rates of Li and B fromthe combined fits of the energy spectra and the ∆T distributions at
1σ, 2σ, 3σ C.L. for (a) showering muons (∆L > 106 p.e.) and (b)non-showering muons (∆L < 106 p.e.) with ∆L < 3 m track cut.
The black points indicate the best-fit parameters.

C.

8

He and 9 Li

4

6
8
10
12
Visible Energy (MeV)

14

16

FIG. 12: Background-subtracted Evis spectrum above 1 MeV, where
signal and background events are taken from 0.002 ≤ ∆T < 1 s
and 5.002 ≤ ∆T < 6 s, respectively. Detection of a neutron capture following a beta decay is required to select events from theβ − + n decay mode. The branching ratios for 9 Li and 8 He are
50.8 ± 0.9% and 16 ± 1%, respectively. The production rates
of 8 He (τ = 171.7 ms, Q = 10.7 MeV) and 9 Li (τ = 257.2 ms,
Q = 13.6 MeV) are estimated to be 1.0 ± 0.5 kton−1 day−1 and
2.8 ± 0.2 kton−1 day−1 , respectively, from simultaneously fitting theEvis spectrum and the ∆T distribution shown in the inset. The fit tothe ∆T distribution has a χ2 /d.o.f.= 95/97.

8

He kton-1day-1

2.5

8

He (τ = 171.7 ms, Q = 10.7 MeV) [41] and 9 Li
(τ = 257.2 ms, Q = 13.6 MeV) [41] β − decay candidate
events are selected according to the cuts 1 ≤ Evis < 13 MeV
and ∆T < 10 s, and by the detection of a neutron following the β − decay event. 8 He decays to neutron-unstable
excited states of 8 Li with a 16 ± 1% branching ratio [41],and 9 Li decays to neutron-unstable excited states of 9 Be
with a 50.8 ± 0.9% branching ratio [41]. The neutron is
identified by the 2.225 MeV γ-ray from radiative captureon 1 H (1.8 ≤ Evis < 2.6 MeV). The γ-ray is required
to be within 200 cm and 1.0 ms of the 8 He/9 Li β − decay
candidate. Finally, the 8 He/9 Li analysis is performed using a
5.5-m-radius spherical fiducial volume to reduce the numberof accidental coincidences between the β − decay candidate
and external γ-ray backgrounds near the balloon.
The inset in Fig. 12 shows the ∆T distribution for the
events that satisfy the criteria outlined above. Figure 12 also
shows the residual Evis distribution corresponding to the subtraction of a background spectrum in the 5.002 ≤ ∆T < 6 swindow from a signal in the 0.002 ≤ ∆T < 1 s window. The
expected Evis distributions for 8 He and 9 Li are calculated byincorporating the KamLAND response and adjusting for theenergy deposited by the thermalizing neutron from 8 He or 9 Li
decay. N (8 He) and N (9 Li) are determined from a simultaneous binned maximum likelihood fit to the ∆T distribution anda chi-square fit to the Evis distribution. For the fit to the Evis
distribution, the uncertainty in the energy scale parameters aretreated in the same manner as the 8 Li/8 B analysis described inSec. V B.

2

2

2.5
(a)

1.5
1
0.5
00

He kton-1day-1

Y (8 Li) = (12.2 ± 2.6) × 10−7 µ−1 g−1 cm2
and
Y (8 B) = (8.4 ± 2.4) × 10−7 µ−1 g−1 cm2 .
The isotope
production rates are R(8 Li) = 15.6 ± 3.2 kton−1 day−1 andR(8 B) = 10.7 ± 2.9 kton−1 day−1 . The contour plots in
Fig. 11 show the correlation between 8 Li and 8 B. Due to theirsimilar lifetimes, 8 Li and 8 B are identified primarily by theirenergy spectra.

0

8

8

1
2
3
9
-1
-1
Li kton day

(b)

2
1.5
1
0.5
00

1
2
3
-1
-1
Li kton day

9FIG. 13: Allowed regions for the production rates of 9 Li and 8 Hefrom the combined fits of the energy spectra and the ∆T distributionsat 1σ, 2σ, 3σ C.L. for (a) showering muons (∆L > 106 p.e.) and (b)non-showering muons (∆L < 106 p.e.) with ∆L < 3 m track cut.
The black points indicate the best-fit parameters.

The procedure for calculating the 8 He and 9 Li selection efficiency is the same as for the 8 Li/8 B efficiency
analysis (Sec. V B), except for the correction for theneutron detection requirement, which is calculated with
the Geant4-based Monte Carlo simulation described inSec. VI A. The resultant efficiencies (8 He) = 14.9 ± 1.0%
and (9 Li) = 46.1 ± 1.1% include the appropriate β-nbranching fractions. Since a reduced volume is used inthis analysis, the 1.6% fiducial volume uncertainty from
Ref. [44] is included in the above efficiencies. The re-

10

100
50

6000

× 102

5000

6
5

Events/5 min

150

8
7

2

4 6 8 10∆T (s)
Total
9
C
8
Li, 9Li
8
B

Events/0.05 MeV

Events/0.1 s

Events/0.5 MeV

200

4000

× 103
5
4
3
100 200∆T (min)

3000

300
11

C

2000
1000

0
12

13

14 15 16 17 18
Visible Energy (MeV)19

20

FIG. 14:Background-subtracted Evis spectrum above
12 MeV, where signal and background events are taken from
0.2 ≤ ∆T < 0.6 s and 10.2 ≤ ∆T < 10.6 s, respectively. The
production rate of 9 C (τ = 182.5 ms, Q = 16.5 MeV) is estimatedto be 3.8 ± 1.5 kton−1 day−1 from a simultaneous fit to the Evis
spectrum and the ∆T distribution shown in the inset. The fit to the∆T distribution has a χ2 /d.o.f.= 97/96.

sulting yields are Y (8 He) = (0.7 ± 0.4) × 10−7 µ−1 g−1 cm2and Y (9 Li) = (2.2 ± 0.2) × 10−7 µ−1 g−1 cm2 . The production rates are R(8 He) = 1.0 ± 0.5 kton−1 day−1 andR(9 Li) = 2.8 ± 0.2 kton−1 day−1 .
The contour plots in
Fig. 13 show the correlation between 9 Li and 8 He.

D.

9

C

The inset in Fig. 14 shows the ∆T distribution for all
events with visible energy 12 ≤ Evis < 20 MeV. The analysis region (0.2 ≤ ∆T < 0.6 s) contains events from 9 C(τ = 182.5 ms, Q = 16.5 MeV) [41] β + decay. N (9 C) is determined from a simultaneous binned maximum likelihood
fit to the ∆T distribution and a chi-square fit to the Evis
distribution. The uncertainty in the energy scale parameters are treated in the same manner as the 8 Li/8 B analysis described in Sec. V B. In this fit, 8 Li, 8 B, and 9 Li aretreated as possible contaminants, the amounts are constrained
to the values obtained in the previously described analyses.
This constraint includes the correlation between 8 Li and 8 B
shown Fig. 11. By integrating the theoretical 9 C Evis spectrum, we obtain the efficiency for the 12 ≤ Evis < 20 MeV
cut of (9 C) = 7.2 ± 1.0%. Combining this with the above
results gives Y (9 C) = (3.0 ± 1.2) × 10−7 µ−1 g−1 cm2 andR(9 C) = 3.8 ± 1.5 kton−1 day−1 .

E.

11

C

The production of 11 C (β + -decay, τ = 29.4 min,
Q = 1.98 MeV) [40] through muon-initiated spallation

0
1.4

1.5

1.6
1.7
1.8
Visible Energy (MeV)

1.9

2

FIG. 15: Background-subtracted Evis spectrum above 1.4 MeV,
where signal and background are selected from 5 ≤ ∆T < 90 min
and 185 ≤ ∆T < 270 min, respectively. The production rateof 11 C (τ = 29.4 min, Q = 1.98 MeV) is estimated to be1106 ± 178 kton−1 day−1 by fitting Eqn. 6 to the ∆T distribution shown in the inset. The fit to the ∆T distribution has aχ2 /d.o.f.= 76/57.

is usually accompanied by a neutron, allowing identificationby the triple coincidence of the primary muon, the spallation
neutron, and the subsequent β + [13, 45]. The 11 C β + decays
are selected in the range 1.4 ≤ Evis < 2.0 MeV and arepreceded by a detected muon that is accompanied by at least
one neutron capture, identified by the 2.225 MeV γ-ray from
capture on 1 H. The γ-ray is required to be in the time window
10 ≤ ∆T < 2500 µs relative to the muon. In order to reducethe background, a 7-m-diameter fiducial volume is used.
To avoid inefficiencies from run boundaries and the long
lifetime of 11 C, the first 5 hours of the typically 24-hour longrun are not used in the selection of the 11 C candidates. The
number of muon-11 C coincidences is extracted from the ∆T
distribution for all events that meet the criteria, shown in theinset of Fig. 15.
The efficiency determination takes into account the visibleenergy range for 11 C β + -decay, 22.7 ± 3.6%, and the previously discussed neutron detection efficiency (Sec. IV). The
efficiency also takes into account a correction for 11 C production modes, designated invisible modes, which do not produce neutrons [13, 45]. To measure this correction, muon11
C event pairs were selected with and without the neutron
requirement for a subset of the data where the 11 C candidate is required to be within 50 cm of the muon track; restricting the study to a subset of the data mitigated the reduced signal-to-background ratio associated with relaxing
the neutron requirement. The number of muon-11 C coincidences in each case was extracted from a fit of Eqn. 6 tothe corresponding ∆T distribution. The visible mode efficiency, vis , taken as the ratio of the number of muon-11 Cpairs with one or more neutrons to the number of muon11
C pairs without the neutron requirement, is 88.4 ± 2.4%.

11

50
102
10
100 200∆T (s)

100

Events/4 s

150

40
Events/0.5 MeV

Events/5 s

Events/0.1 MeV

200

300
10

C

30

70
60
50
40
30
50

20
10

100 150∆T (s)
11
Be

0

50

-10

0
2

6

2.5
3
3.5
Visible Energy (MeV)

4

FIG. 16: Background-subtracted Evis spectrum above 2 MeV,
where signal and background are selected using 10 ≤ ∆T < 90 s
and 190 ≤ ∆T < 270 s, respectively. The production rate of 10 C
(τ = 27.8 s, Q = 3.65 MeV) is 21.1 ± 1.8 kton−1 day−1 , as determined by fitting Eqn. 6 to the ∆T distribution shown in the inset.
The fit to the ∆T distribution has a (χ2 /d.o.f.= 80/56).

Applying the correction for post muon electronics effects,
the visible mode fraction is 96.3 ± 2.0%, consistent with
Ref. [45], which obtains vis = 95.6% for 285 GeV muons.
Due to the relatively long 11 C lifetime, we have also considered the effect of diffusion. An analysis of 222 Rn that was
accidentally introduced into the center of KamLAND during the deployment of a calibration device shows that the
diffusion speed is approximately 1 mm/h. From this study
the effect of 11 C diffusion on the efficiency is estimatedto be less than 0.5%. Combining this with the above results gives Y (11 C) = (866 ± 153) × 10−7 µ−1 g−1 cm2 andR(11 C) = 1106 ± 178 kton−1 day−1 .

F.

10

C

As with 11 C, the production of 10 C(β + -decay, τ = 27.8 s,
Q = 3.65 MeV) [41] through muon-initiated spallation is usually accompanied by a neutron, so the selection criterion requiring a triple coincidence of the primary muon, the neutron,and 10 C candidate is used. The neutron is identified by the2.225 MeV n+1 H capture γ-ray. The number of 10 C candidates is determined from fitting Eqn. 6 to the ∆T distribution for all events identified as 10 C, shown in the inset in
Fig. 16. 11 Be is a potential background for this 10 C analysis, but the correction to the 10 C yield is estimated to beless than 1% to due to the low 11 Be production rate and theneutron coincidence requirement. The efficiency for the visible energy cut 2.0 ≤ Evis < 4.0 MeV is 73.5 ± 3.2%. The
visible mode efficiency is vis = 90.7 ± 5.5% after correcting for the electronics effects following muons. The final
efficiency is 89.6 ± 5.5%. The resulting isotope yield is81012Visible Energy (MeV)

14

FIG. 17: Background-subtracted Evis spectrum above 5.5 MeV for
showering muons (∆L > 106 p.e.), where signal and background are
selected by 8 ≤ ∆T < 60 s and 408 ≤ ∆T < 460 s, respectively.
The production rate of 11 Be (τ = 19.9 s, Q = 11.5 MeV) for showering muons is estimated to be 1.0 ± 0.2 kton−1 day−1 by fitting
Eqn. 6 to the ∆T distribution shown in the inset (χ2 /d.o.f.= 37/41).

Y (10 C) = (16.5 ± 1.9) × 10−7 µ−1 g−1 cm2 and the production rate is R(10 C) = 21.1 ± 1.8 kton−1 day−1 .

G.

11

Be

The 11 Be β − decay (τ = 19.9 s, Q = 11.5 MeV) [40]events are selected according to 5.5 ≤ Evis < 16.0 MeV. The
Evis cut efficiency is 63.4%. The inset in Fig. 17 showsthe ∆T distributions for the events only after showeringmuons. For non-showering muons, a tighter muon track cut∆L < 1 m is applied in order to reduce the background rate.
The track cut efficiency is estimated from the 12 B candidates
using an analysis similar to Sec. V B. The resulting isotope
yield is Y (11 Be) = (1.1 ± 0.2) × 10−7 µ−1 g−1 cm2 and production rate is R(11 Be) = 1.4 ± 0.3 kton−1 day−1 .

VI.

MONTE CARLO SIMULATIONSThe Geant4 and FLUKA simulations are used to reproduce
the measurements from KamLAND. While Geant4 is used
only to simulate neutron production, both neutron and lightisotope production are tested with FLUKA.

A. Geant4

Geant4 is a widely used toolkit for performing particle
tracking simulations on an event-by-event basis. A description of the available physics processes included is given in
Refs. [17] and [18]. Here we compare the Geant4 (version
9.1) prediction for neutron yield by spallation with the results

90
80
70
60
50
40
30
20
10
0

Measurement
MC (tracked muon)
MC (untracked muon)

2.5
2Events/m 3

Events/7 MeV

12Measurement
MC (tracked muon)
MC (untracked muon)

1.5
1
0.5

0

100

200 300 400 500
Visible Energy (MeV)

600

700

FIG. 18: Evis distribution of the prompt events from candidates
identified as neutrons produced by muon-induced spallation in thematerial outside of the KamLAND ID.

obtained in Sec. IV. We use the physics list QGS BIC, developed by the Geant4 group to support the binary cascade(BIC) model at lower energies (below 10 GeV for p and n, andbelow 1.2 GeV for π). This treatment is also appropriate forthe simulation of interactions of nucleons and ions. At higher
energies, a quark-gluon string (QGS) model is applied for the
hadronic interactions. Neutron elastic and inelastic interactions below 20 MeV are described by a high-precision datadriven model (NeutronHP). The G4EmExtraPhysics
physics list is also used to model the photo-nuclear and muonnuclear interaction processes, which dominate the neutron
production by muons in the simulation.
To estimate the neutron production yield as a function ofmuon energy in Geant4, mono-energetic muons of several
energies are injected at the center of a generic hydrocarbon
block of thickness 40 m. The region more than 10 m awayfrom the edges of the block is analyzed to avoid boundary effects. As shown in Fig. 20, the neutron production yields predicted by Geant4 are systematically lower than experiment,with the exception of point (F) from LVD [12]. These results
are consistent with previous work [46].
A Monte Carlo simulation based upon Geant4 andMUSIC (described in Sec. III) is used to study neutrons produced by muon-induced spallation in the material outside of
the KamLAND ID. Some of these neutrons have sufficient energy to enter the ID where they thermalize and capture. They
can be identified by the coincidence of a prompt signal (forexample from n + p elastic scattering) and a delayed signal
from the capture γ-ray. This is the same inverse beta decay
reaction signature used for ν e detection, ν e + p → e+ + n,
where the e+ is the prompt signal and the γ-ray from neutron capture is the delayed signal and therefore a potential
background. These neutrons are also a background for dark
matter experiments that employ nuclear recoils as a detection
method.
The primary purpose of this Monte Carlo simulation is0
300 350 400 450 500 550 600 650 700
Prompt Event Radius (cm)FIG. 19: Radial distribution (from the center of KamLAND)
of the prompt event from candidates identified as neutrons produced by muon-induced spallation in the material outside of the
KamLAND ID. Assuming an exponential distribution, a fit to themeasured data (red dashed lines) yields an attenuation length of
70 ± 2 g/cm2 . A similar fit to the Monte Carlo (black dotted lines)
yields 69 ± 2 g/cm2 .

to estimate the rate of untagged fast neutrons, i.e. neutrons
produced by muons where the muon is undetected by the
KamLAND ID. A few of these muons are detected by the
OD, either from the Čerenkov radiation produced by the muon
itself (tracked muons), or by accompanying electromagnetic
and hadronic showers that enter the OD (untracked muons).
A prompt signal in coincidence with these tracked and untracked muons followed by a delayed capture gamma signal
identifies candidates for untagged fast neutrons.
Fast neutrons generated by untracked muons are produced
primarily in the rock surrounding the OD, whereas fast neutrons from tracked muons are primarily produced in the waterof the OD. It is shown in the Monte Carlo that the tracked anduntracked muons give distinguishable OD visible energy distributions. A comparison between Monte Carlo and measurement of the distribution of the number of OD PMTs with signals above threshold for tracked and untracked muons reveals
a deficiency of fast neutrons from untracked muons, consistent
with the underproduction of neutrons by Geant4 in concrete
reported in Ref. [47].
Figure 18 shows the Evis distribution of the prompt eventsfrom the Monte Carlo simulation of tracked and untracked
muons compared to the measured data. The measured data
and the Monte Carlo simulation corresponds to an equal live
time exposure of 1368 days. Figure 19 shows the radial distribution of the prompt events relative to the KamLAND center.
Both the measured data and the Geant4 simulation exhibit an
exponential attenuation of the neutrons as they penetrate further into the detector, with the simulation yielding an attenuation length of 69 ± 2 g/cm2 , consistent with the measurement
70 ± 2 g/cm2 .

13

Neutron Production Yield (µ-1g-1cm2)

10-3

simulation. The results, including the primary production process and power law exponent, are summarized in Table IV.
For some isotopes the primary production process is muchlarger than any secondary processes, as in the case of 12 B. For
other isotopes the primary production process is only slightly
larger than the secondary processes, as is the case for 9 Li. The
isotopes produced primarily by γ interactions, 11 C and 10 C,
show the weakest dependence on muon energy. In comparison, 12 N and 13 N, where the primary production mechanism
is by p interactions, show the strongest dependence on muonenergy.

(G)this meas.
(E)
(D)

10-4

(F)
(B)
(C)
(A)

10-50

100

200
300
400
Muon Energy (GeV)

500

FIG. 20: Neutron production yield in liquid scintillator as a functionof muon energy. The red star point shows the KamLAND result,
Yn = (2.8 ± 0.3) × 10−4 µ−1 g−1 cm2 , for 260 ± 8 GeV. Other
points show the results from experiments at (A) 20 m.w.e. [48],(B) 25 m.w.e. [49], (C) 32 m.w.e. [50], (D) 316 m.w.e. [49], (E)570 m.w.e. [51], (F) 3000 m.w.e. [12], and (G) 5200 m.w.e. [11].
The blue square and green triangle show the Geant4 and FLUKA
Monte Carlo predictions, respectively, from monochromatic muon
beams.

The use of a mono-energetic µ− beam overestimates the
production of neutrons and light isotopes. Simulations using a mono-energetic µ+ beam and a beam with the energy
spectrum from Ref. [33] were also run. The simulations showthat the production yield for µ+ relative µ− is on average
0.96±0.01 for the light isotopes. This reduction is expected
since µ− may capture creating spallation products while µ+
may not. This ratio, combined with the µ+ to µ− ratio at
KamLAND, leads to a correction to the flux of 0.98±0.06
for light isotopes and 0.981±0.005 for neutrons. The reduced
production yield due to averaging over the muon spectrum is
on average 0.92±0.02 for the light isotopes, which is slightly
higher than the correction factor suggested by Ref. [10]. The
results of the FLUKA simulations for KamLAND presented inTable V include these two corrections.

VII.
B.

DISCUSSION

FLUKA

FLUKA is a mature code that models nuclear and particle physics processes from thermal neutrons to heavy ion
collisions [15, 16]. It has been used previously to modelmuon-initiated spallation in liquid scintillator [31, 46, 52,53]. We use FLUKA Version 2006.3b to model neutron
and light isotope production from muon-initiated spallationin KamLAND. A 40 m radius by 40 m high cylinder of
KamLAND liquid scintillator is used in the simulation; theconcentric inner cylinder of 20 m radius and 20 m length isused for analysis.
To estimate neutron production yield as a function of muonenergy in FLUKA, mono-energetic beams of µ− ranging from10 GeV to 350 GeV were simulated as in Ref. [31, 46, 52, 53].
Care is taken not to double count neutrons involved in reactions like (n, 2n). The results of this simulation are includedin Fig. 20. The neutron production yield of this FLUKA simulation is 10% lower than previous work by [31, 46, 52, 53],
but the power law dependance on muon energy (Eµα , where
α = 0.77) is consistent. Different scintillator compositions
were studied, but could not explain the deficit. This deficit isinsignificant compared to the discrepancy between these simulations and the data.
The production of light isotopes was studied in the same

The isotope production yields from muon-initiated spallation in a liquid scintillator target were investigated at CERN
by earlier experiment [10] using the SPS muon beam withmuon energies of 100 GeV and 190 GeV. Based on those crosssection measurements and the predicted muon energy spectrum at the KamLAND site, we calculated the isotope production rates by extrapolation, assuming a power-law of themuon energy, Eµα . The mean muon energy at KamLAND is
260 ± 8 GeV, so KamLAND provides data to test the extrapolation method at this energy.
The production yields for the isotopes from muon-initiatedspallation in KamLAND are provided in Table V. On average,
the yields from the showering muons (∼15% of all muons),
whose excess light yield parameter (∆L, Eqn. 2) is greaterthan 106 p.e. (∼3 GeV), constitute 70 ± 2% of the yield fromall muons. The production yield for 11 C is the largest, and its
measured yield is larger than the FLUKA calculation by a factor of ∼2. The Borexino collaboration also reported a similar
discrepancy [2, 3], that is consistent with what is observed in
KamLAND. Some measured production yields, such as 8 Li
and 10 C, deviate significantly from estimates based on themuon beam experiment, indicating that perhaps estimation by
extrapolation is not sufficient. All isotope yields are consistent
within an order of magnitude.

14
TABLE IV: Simulation of neutron and light isotope production in the KamLAND LS by muon-initiated spallation with FLUKA: the isotope
production yield by mono-energetic (260 GeV) µ− ; the ratio of the production yields by mono-energetic µ+ compared to µ− ; the ratio of
the production yields by a µ− spectrum that matches Fig. 6 (KamLAND curve) in Ref. [33] compared to mono-energetic µ− ; the power law
exponent for the production yield Y (Eµ ) ∝ Eµα from a fit to the yields from mono-energetic µ− with 10 ≤ Eµ ≤ 350 GeV; and the primaryprocess for producing the isotope. The uncertainties are statistical.
Simulated Production Yield
−7

(×10
n
C
7
Be
10
Be
12
B
8
Li
10
C
6
He
8
B
9
Li
9
C
12
N
11
Be
8
He
13
B
15
O
13
N
11

−1 −1

2

µ g cm )2344 ± 4460.8 ± 1.7
116.8 ± 0.9
44.63 ± 0.53
30.85 ± 0.44
23.42 ± 0.39
21.13 ± 0.37
13.40 ± 0.29
6.40 ± 0.20
3.51 ± 0.15
1.49 ± 0.10
0.86 ± 0.07
0.94 ± 0.08
0.35 ± 0.05
0.31 ± 0.04
0.05 ± 0.02
0.06 ± 0.02

Ratio of Simulated Production Yields forµ+ /µ−
0.969 ± 0.002
0.971 ± 0.005
0.986 ± 0.011
0.960 ± 0.018
0.970 ± 0.021
0.927 ± 0.026
0.982 ± 0.025
0.916 ± 0.035
0.996 ± 0.045
0.856 ± 0.074
0.850 ± 0.114
0.963 ± 0.128
0.842 ± 0.145
0.964 ± 0.200
1.020 ± 0.197
1.250 ± 0.379
1.500 ± 0.272

Spectrum / Mono-Energetic
0.912 ± 0.003
0.913 ± 0.006
0.945 ± 0.011
0.891 ± 0.019
0.936 ± 0.022
0.936 ± 0.025
0.915 ± 0.027
0.918 ± 0.035
0.915 ± 0.050
0.842 ± 0.078
0.949 ± 0.102
1.006 ± 0.120
0.804 ± 0.161
0.576 ± 0.372
1.062 ± 0.176
1.635 ± 0.234
1.190 ± 0.401

Power Law Exponent

Primary Process

0.779 ± 0.001
0.703 ± 0.002
0.684 ± 0.004
0.825 ± 0.007
0.828 ± 0.009
0.821 ± 0.010
0.810 ± 0.010
0.818 ± 0.013
0.804 ± 0.019
0.801 ± 0.026
0.772 ± 0.039
0.921 ± 0.045
0.753 ± 0.051
0.926 ± 0.078
0.742 ± 0.075
0.793 ± 0.244
1.120 ± 0.220

π − +1 H,12 C
12C(γ, n)12C(γ, nα)12C(n,3 He)
12
C(n, p)12
C(n, pα)12C(π + , np)
12C(n, 2p3 He)
12C(π + ,2 H2 H)
12C(π − ,3 He)
12C(π + ,3 H)
12C(p, n)12
C(n, 2p)12C(π − ,3 H)
13C(n, p)16
O(γ, n)13
C(p, n)TABLE V: Summary of the neutron and isotope production yields from muon-initiated spallation in KamLAND. The results of the FLUKA
calculation shown in this table include corrections for the muon spectrum and the µ+ /µ− composition of the cosmic-ray muon flux.

n
B
12
N
8
Li
8
B
9
C
8
He
9
Li
11
C
10
C
11
Be
6
He
7
Be
12

Lifetime in
KamLAND LS
207.5 µs
29.1 ms
15.9 ms
1.21 s1.11 s182.5 ms
171.7 ms
257.2 ms
29.4 min27.8 s19.9 s1.16 s76.9 day

Radiation Energy
2.225 MeV (capt. γ)
13.4 MeV (β − )17.3 MeV (β + )16.0 MeV (β − α)18.0 MeV (β + α)16.5 MeV (β + )10.7 MeV (β − γn)13.6 MeV (β − γn)1.98 MeV (β + )3.65 MeV (β + γ)
11.5 MeV (β − )3.51 MeV (β − )0.478 MeV (EC γ)

VIII.

Spallation Production Yield (×10−7 µ−1 g−1 cm2 )Hagner, et al. [10]FLUKA calc.
this measurement—2097 ± 132787 ± 311—
27.8 ± 1.9
42.9 ± 3.3—
0.77 ± 0.08
1.8 ± 0.4
1.9 ± 0.8
21.1 ± 1.4
12.2 ± 2.6
3.3 ± 1.0
5.77 ± 0.42
8.4 ± 2.4
2.3 ± 0.9
1.35 ± 0.12
3.0 ± 1.2o0.32 ± 0.05
0.7 ± 0.4
1.0 ± 0.3
3.16 ± 0.25
2.2 ± 0.2421 ± 68416 ± 27866 ± 153
54 ± 1219.1 ± 1.3
16.5 ± 1.9< 1.1
0.84 ± 0.09
1.1 ± 0.2
7.5 ± 1.5
12.08 ± 0.83—
107 ± 21105.3 ± 6.9—

SUMMARY

We have analyzed KamLAND data to measure production yields of radioactive isotopes and neutrons
through muon-initiated spallation in liquid scintillator.
The neutron production yield is evaluated to be
Yn = (2.8 ± 0.3) × 10−4 µ−1 g−1 cm2 , which is higher thanFraction from showering µ
this measurement
64 ± 5%68 ± 2%77 ± 14%65 ± 17%78 ± 23%91 ± 32%
76 ± 45%77 ± 6%62 ± 10%76 ± 6%74 ± 12%——the expectation from Monte Carlo simulations based on
Geant4 and FLUKA. Some isotope production yields arefound to be inconsistent with extrapolations – based on a
power law dependance with respect to muon energy – ofresults from muon beam experiments.

15
IX.

ACKNOWLEDGMENTS

The KamLAND experiment is supported by the COE program under grant 09CE2003 of the Japanese Ministry of Education, Culture, Sports, Science and Technology; World Premier International Research Center Initiative (WPI Initiative),
MEXT, Japan; and under the United States Department of[1] G. Bellini, J. Benziger, S. Bonetti, M. B. Avanzini, B. Caccianiga, L. Cadonati, F. Calaprice, C. Carraro, A. Chavarria, F. Dalnoki-Veress, et al. (Borexino Collaboration), astroph/0808.2868v1.
[2] C. Arpesella, G. Bellini, J. Benziger, S. Bonetti, B. Caccianiga, F. Calaprice, F. Dalnoki-Veress, D. D’Angelo, H. de Kerret,
A. Derbin, et al. (Borexino Collaboration), Phys. Lett. B 658,
101 (2008).
[3] C. Arpesella, H. O. Back, M. Balata, G. Bellini, J. Benziger,
S. Bonetti, A. Brigatti, B. Caccianiga, L. Cadonati, F. Calaprice,
et al. (Borexino Collaboration), Phys. Rev. Lett. 101, 091302(2008).
[4] Y. Hirano, T. Kishimoto, I. Ogawa, R. Hazama, S. Umehara,
K. Matsuoka, G. Ito, and Y. Tsubota, J. Phys. Conf. Ser. 120,052053 (2008).
[5] T. Kishimoto, talk at the International Workshop on ”Double
Beta Decay and Neutrinos”, June 11-13, 2007, Osaka, available
online at http://dbd07.phys.sci.osaka-u.ac.jp/ (2007).
[6] K. Zuber, Am. Inst. Phys. Conf. Proc. 942, 101 (2007).
[7] M. C. Chen, Earth, Moon, and Planets 99, 221 (2006).
[8] C. Grieb, J. M. Link, M. L. Pitt, R. S. Raghavan, D. Rountree,
and R. B. Vogelaar, hep-ex/0705.2769v1.
[9] T. M. Undagoitia, F. von Feilitzsch, M. Göger-Neff, L. Oberauer, W. Potzel, A. Ulrich, J. Winter, and M. Wurm, J. Phys.
Conf. Ser. 120, 052018 (2008).
[10] T. Hagner, R. von Hentig, B. Heisinger, L. Oberauer,
S. Schönert, F. von Feilitzsch, and E. Nolte, Astropart. Phys.
14, 33 (2000).
[11] M. Aglietta et al., Nuovo Cimento Soc. Ital. Fis. C 12, 467(1989).
[12] M. Aglietta, E. Alyea, P. Antonioli, G. Badino, G. Bari,
M. Basile, V. Berezinsky, F. Bersani, M. Bertaina, R. Bertoni,
et al., Phys. Atomic Nuclei 66, 123 (2003).
[13] H. Back, M. Balata, G. Bellini, J. Benziger, S. Bonetti, B. Caccianiga, F. Calaprice, D. D’Angelo, A. de Bellefon, H. de Kerret, et al. (Borexino Collaboration), Phys. Rev. C 74, 045805(2006).
[14] P. Antonioli, C. Ghetti, E. V. Korolkova, V. A. Kudryavtsev, and
G. Sartorelli, Astropart. Phys. 7, 357 (1997).
[15] A. Fassò, A. Ferrari, S. Roesler, P. Sala, G. Battistoni, F. Cerutti,
E. Gadioli, M. Garzelli, F. Ballarini, A. Ottolenghi, et al., hepph/0306267.
[16] A. Ferrari, P. R. Sala, A. Fassó, and J. Ranft, FLUKA: A
multi-particle transport code (program version 2005) (CERN,
Geneva, 2005).
[17] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo,
P. Arce, M. Asai, D. Axen, S. Banerjee, G. Barrand, et al., Nucl.
Instr. Meth. A 506, 250 (2003).
[18] J. Allison, K. Amako, J. Apostolakis, H. Araujo,
P. Arce Dubois, M. Asai, G. Barrand, R. Capra, S. Chauvie, R. Chytracek, et al., IEEE Trans. Nucl. Sci. 53, 270

Energy (DOE) grants DEFG03-00ER41138 and DE-AC0205CH11231; as well as other DOE grants to individual institutions. We are grateful to the Kamioka Mining and Smelting
Company for supporting the activities in the mine. We would
also like to thank V. A. Kudryavtsev and M. Spurio for helpfulcomments.

(2006).
[19] R. Battino, T. R. Rettich, and T. Tominaga, J. Phys. Chem. Ref.
Data 13, 563 (1984).
[20] P. J. Hesse, R. Battino, P. Scharlin, and E. Wilhelm, J. Chem.
Eng. Data 41, 195 (1996).
[21] H. Kume, S. Sawaki, M. Ito, K. Arisaka, T. Kajita,
A. Nishimura, and A. Suzuki, Nucl. Instr. Meth. Phys. Res. 205,
443 (1983).
[22] S. Kleinfelder, IEEE Trans. Nucl. Sci. 50, 955 (Aug. 2003).
[23] D. W. McKee, J. K. Busenitz, and I. Ostrovskiy, Nucl. Instr.
Meth. A 587, 272 (2008).
[24] B. E. Berger, J. Busenitz, T. Classen, M. P. Decowski, D. A.
Dwyer, G. Elor, A. Frank, S. J. Freedman, B. K. Fujikawa,
M. Galloway, et al. (KamLAND Collaboration), JINST 4,
P04017 (2009).
[25] J. B. Birks, Proc. Phys. Soc. A64, 874 (1951).
[26] J. B. Birks, The Theory and Practice of Scintillation Counting(Pergamon, London, 1964).
[27] Digital Map 50 m Grid (Elevation), Geographical Survey Institute of Japan (1997), unpublished.
[28] Y. Becherini, A. Margiotta, M. Sioli, and M. Spurio, Astropart.
Phys. 25, 1 (2006).
[29] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, and
T. Sanuki, Phys. Rev. D 75, 043006 (2007).
[30] D. Casper, Nucl. Phys. B Proc. Suppl. 112, 161 (2002).
[31] D.-M. Mei and A. Hime, Phys. Rev. D 73, 053004 (2006).
[32] C. Galbiati and J. F. Beacom, Phys. Rev. C 72, 025807 (2005).
[33] A. Tang, G. Horton-Smith, V. A. Kudryavtsev, and A. Tonazzo,
Phys. Rev. D 74, 053007 (2006).
[34] S. Eidelman, K. Hayes, K. Olive, M. Aguilar-Benitez, C. Amsler, D. Asner, K. Babu, R. Barnett, J. Beringer, P. Burchat,et al., Phys. Lett. B 592, 1 (2004).
[35] Tech. Rep., Kamioka Mining & Smelting Company (1977), internal report.
[36] D. E. Groom, N. V. Mokhov, and S. I. Striganov, Atomic Data
and Nuclear Data Tables 78, 183 (2001).
[37] P. H. Barrett, L. M. Bollinger, G. Cocconi, Y. Eisenberg, andK. Greisen, Rev. Mod. Phys. 24, 133 (1952).
[38] S. F. Mughabghab, M. Divadeenam, and N. E. Holden, Neutron
Cross Sections, Volume 1, Neutron Resonance Parameter and
Thermal Cross Sections, Part A Z = 1 − 60 (Academic Press,
New York, 1981).
[39] S. Baker and R. D. Cousins, Nucl. Instr. Meth. Phys. Res. 221,437 (1984).
[40] F. Ajzenberg-Selove, Nucl. Phys. A 506, 1 (1990).
[41] D. R. Tilley, J. H. Kelley, J. L. Godwin, D. J. Millener, J. E.
Purcell, C. G. Sheu, and H. R. Weller, Nucl. Phys. A 745, 155(2004).
[42] W. T. Winter, S. J. Freedman, K. E. Rehm, and J. P. Schiffer,Phys. Rev. C 73, 025503 (2006).
[43] M. Bhattacharya, E. G. Adelberger, and H. E. Swanson, Phys.

16
Rev. C 73, 055802 (2006).
[44] S. Abe, T. Ebihara, S. Enomoto, K. Furuno, Y. Gando,
K. Ichimura, H. Ikeda, K. Inoue, Y. Kibe, Y. Kishimoto,et al. (KamLAND Collaboration), Phys. Rev. Lett. 100, 221803(2008).
[45] C. Galbiati, A. Pocar, D. Franco, A. Ianni, L. Cadonati, and
S. Schönert, Phys. Rev. C 71, 055805 (2005).
[46] H. Araújo, V. Kudryavtsev, N. Spooner, and T. Sumner, Nucl.
Instr. Meth. A 545, 398 (2005).
[47] M. G. Marino, J. A. Detwiler, R. Henning, R. A. Johnson, A. G.
Schubert, and J. F. Wilkerson, Nucl. Instr. Meth. A 582, 611(2007).
[48] R. Hertenberger, M. Chen, and B. L. Dougherty, Phys. Rev. C
52, 3449 (1995).

[49] L. B. Bezrukov, V. I. Beresnev, G. T. Zatsepin, O. G. Ryazhskaya, and L. N. Stepanets, Sov. J. Nucl. Phys. 17, 51 (1973).
[50] F. Boehm, J. Busenitz, B. Cook, G. Gratta, H. Henrikson, J. Kornis, D. Lawrence, K. B. Lee, K. McKinny, L. Miller, et al.,Phys. Rev. D 62, 092005 (2000).
[51] R. I. Enikeev, G. T. Zatsepin, E. V. Koroĺkova, V. A.
Kudryavstev, A. S. Maĺgin, O. G. Ryazhskaya, and F. F.
Khaĺchukov, Sov. J. Nucl. Phys. 46, 883 (1987).
[52] Y.-F. Wang, V. Balic, G. Gratta, A. Fassò, S. Roesler, and
A. Ferrari, Phys. Rev. D 64, 013012 (2001).
[53] V. A. Kudryavtsev, N. J. C. Spooner, and J. E. McMillan, Nucl.
Instr. Meth. A 505, 688 (2003).