{'fraction_non_alphanumeric': 0.027625117542576533, 'fraction_numerical': 0.025911607982446974, 'mean_word_length': 4.74088291746641, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmogenic radionuclide production in NaI(Tl) crystals

The production of long-lived radioactive isotopes in materials due to the exposure to cosmic rays on Earth surface can be an hazard for experiments demanding ultra-low background conditions, typically performed deep underground. Production rates of cosmogenic isotopes in all the materials present in the experimental set-up, as well as the corresponding cosmic rays exposure history, must be both well known in order to assess the relevance of this effect in the achievable sensitivity of a given experiment. Although NaI(Tl) scintillators are being used in experiments aiming at the direct detection of dark matter since the first nineties of the last century, very few data about cosmogenic isotopes production rates have been published up to date. In this work we present data from two 12.5 kg NaI(Tl) detectors, developed in the frame of the ANAIS project, which were installed inside a convenient shielding at the Canfranc Underground Laboratory just after finishing surface exposure to cosmic rays. The very fast start of data taking allowed to identify and quantify isotopes with half-lives of the order of tens of days. Initial activities underground have been measured and then production rates at sea level have been estimated following the history of detectors; values of about a few tens of nuclei per kg and day for Te isotopes and 22Na and of a few hundreds for I isotopes have been found. These are the first direct estimates of production rates of cosmogenic nuclides in NaI crystals. A comparison of the so deduced rates with calculations using typical cosmic neutron flux at sea level and a carefully selected description of excitation functions will be also presented together with an estimate of the corresponding contribution to the background at low and high energies, which can be relevant for experiments aiming at rare events searches.


Introduction
Production of radioactive isotopes by exposure to cosmic rays is a well known issue [1]. Several cosmogenic isotopes are important because of their broad applications in extraterrestrial and terrestrial studies dealing with astrophysics, geophysics, paleontology and archaeology. For experiments searching for rare phenomena like the nuclear Double Beta Decay (DBD) or the interaction of Weakly Interacting Massive Particles (the so-called WIMPs, which could be filling the galactic dark matter halo) requiring detectors working in ultra-low background conditions, long-lived radioactive nuclei induced in materials of the set-up can be a sensitivity limiting hazard. Operating in deep underground locations, using active and passive shields and selecting carefully radiopure materials for the building of the detectors and shields are common practices in this kind of experiments to reduce very efficiently the background entangling the searched signal [2,3]. Cosmogenic radionuclides produced during fabrication, transport and storage of components may be even more important than residual contamination from primordial nuclides and become very problematic, depending on the target. For instance, the poor knowledge of cosmic ray activation in detector materials is highlighted in [4] as one of the three main uncertain nuclear physics aspects of relevance in the direct detection approach pursued to solve the dark matter problem.
One of the most relevant processes in cosmogenic activation is the spallation of nuclei by high energy nucleons, but other reactions like fragmentation, induced fission or capture can be very important for some nuclei. Cosmogenic activation is strongly dependent on the nucleon flux, neutron to proton ratio, and energies available. At sea level, for instance, the flux of neutrons and protons is virtually the same at a few GeV, but the proton to neutron ratio decreases significantly at lower energies because of the charged particles absorption in the atmosphere. For example, at 100 MeV this ratio is about 3% [1]. Summarizing, at sea level, nuclide production is mainly dominated by neutrons at low energies, whereas if materials are flown at high altitude cosmic flux is much greater, energies at play are larger and activation by protons can not be neglected.
In principle, cosmogenic activation can be kept under control by minimizing exposure at surface and storing materials underground, avoiding flights and even using shields against the hadronic component of cosmic rays during surface operation. But since these requirements usually complicate the preparation of experiments (for example, while crystal growth and detectors mounting steps) it would be desirable to have reliable tools to quantify the real danger of exposing the different materials to cosmic rays. Direct measurements, by screening of exposed materials in very low background conditions as those achieved in underground laboratories, and calculations of production rates and yields, following different approaches, have been made for several materials in the context of dark matter, double beta decay and neutrino experiments (see for instance [5] and references therein). Many different studies are available for germanium [6]- [12] and interesting results have been derived also for copper [12,13], stainless steel [14], tellurium [15,16] or neodymium [17,18]. However, systematic studies for other targets are missed.
For NaI, production of 129 I, 125 I and 22 Na was analyzed in [19] in the context of the DAMA/LIBRA experiment and the presence of cosmogenic isotopes has been reported by ANAIS [20,21], DM-Ice [22] and KIMS [23] collaborations, working also with this material. In this work, cosmogenic yields induced on the Earth surface in two large NaI(Tl) crystals produced by Alpha Spectra 1 for the ANAIS (Annual modulation with NaI Scintillators) experiment have been quantified. The ANAIS project [24] is intended to search for dark matter annual modulation with 250 kg of ultrapure NaI(Tl) scintillators at the Canfranc Underground Laboratory (LSC) in Spain; the aim is to provide a model-independent confirmation of the annual modulation positive signal reported by the DAMA collaboration [25] using the same target and technique. Several prototypes have been developed and operated in Canfranc [26][27][28]; results presented here come from the so-called ANAIS-25 set-up [20]. The ultra-low background conditions of ANAIS and a prompt start of data taking after moving the detectors underground were essential to carry out the analysis of the cosmogenic activation. The structure of the paper is the following. Section 2 describes the experimental set-up and the measurements taken, while the data analysis performed and the production rates at sea level derived for the cosmogenic isotopes identified in the measurements are presented in section 3. A calculation of the corresponding production rates has been attempted using carefully selected excitation functions and the history of the raw material the detectors consist of; such calculations are summarized in section 4, together with the comparison with the experimental rates and corresponding discussion. Finally, conclusions are presented in section 5.

Experimental set-up and measurements
Two NaI(Tl) detectors built in collaboration with Alpha Spectra form the ANAIS-25 set-up [20]. Each detector consists of a cylindrical 12.5 kg NaI(Tl) crystal, with 4.75" diameter and 11.75" length. They were grown with selected ultrapure NaI powder and encapsulated in OFHC copper with two synthetic quartz windows allowing the coupling to photomultiplier tubes (PMTs) in a second step. Two Hamamatsu R12669SEL2 PMTs were used for one crystal and Hamamatsu R11065SEL PMTs for the other; we will refer to both detectors  as D0 and D1, respectively, in the following. Only white Teflon was used as light diffuser, wrapping the crystal, inside the copper encapsulation. A Mylar window allows to calibrate at low energy both detectors. ANAIS-25 detectors are being operated at hall B of LSC under 2450 m.w.e., inside a shielding consisting of 10 cm archaeological lead plus 20 cm low activity lead, all enclosed in a PVC box tightly closed and continuously flushed with boil-off nitrogen. Figure 1 shows a picture of detector D0 and a drawing of the ANAIS-25 set-up is depicted in figure 2.
Each PMT charge output signal is separately processed [20]; each one is divided into a trigger signal, a signal going to the digitizer, and several ones differently amplified /attenuated, and fed into QDC (charge-to-digital converter) module channels to be integrated in a 1 µs window. Low energy (LE) and high energy (HE) ranges, below and above ∼200 keV respectively, have been considered in the data analysis presented here. The building of the spectra was done by software (off-line) by adding the signals from both PMTs. Data triggering is done in logical OR mode between the two detectors and in logical AND between the two PMT signals from each detector.
The main goal of ANAIS-25 set-up was to assess the total background of the detectors and, in particular, to determine precisely the 40 K content of the crystals, since it was found to be the dominant background in previous ANAIS prototypes [26]. The technique used to estimate the bulk 40 K activity in the crystals is the measurement in coincidence: one detector measures the total energy (3.2 keV) released by the X-ray/Auger electrons emissions of argon following the K-shell electron capture (EC) of 40 K, while the other one measures the fully absorbed energy (1460.8 keV) of the high energy gamma escaping from the former. Averaging the results for the two crystals, an activity of 1.25±0.11 mBq/kg (41.7±3.7 ppb of natural potassium) was derived [27]; this value is one order of magnitude better than that found in other crystals, but still a factor of 2 higher than the required purity. Background assessment is in progress, as far as cosmogenic contribution to the background is still decaying. On the other hand, an outstanding light collection efficiency has been obtained in ANAIS-25 modules, with measured values of 16.13±0.66 and 12.58±0.13 phe/keV for D0 and D1 detectors respectively [20]; this results in a very good energy resolution which has helped in the identification of cosmogenically produced isotopes.
ANAIS-25 crystals were built in Alpha Spectra facilities in Colorado and then shipped by boat from the US to Spain, avoiding air travel. They arrived at LSC in the evening of 27 th November 2012 and were immediately stored underground; data taking started only three days later, after coupling PMTs in the LSC clean room, placing the detectors inside the shielding and preliminary testing of electrical connections and general performance. This prompt commissioning of ANAIS-25 allowed to observe short-lived isotopes and to properly quantify induced activation in those long-lived. Data analyzed for this study correspond to three data sets. Data sets I and II span altogether for 210 days, from the beginning of the data taking to the end of June 2013; we will refer to set I for data obtained until February 2013 and to set II for those taken from March to June 2013. They were acquired in slightly different gain conditions, which implied that some signals present in the LE or HE spectra in one set were not available in the other. Set III includes 88 days of data from March to June 2014, once the contribution of most of the cosmogenic nuclei is significantly reduced, allowing to identify longer-living products. Figure 3 compares the spectra, both in the low and high energy regions, evaluated in the beginning of data set I and about fifteen months afterwards, during data set III, over 13.0 and 18.5 days, respectively. Several lines can be clearly attributed to cosmogenic activation, as figure 3 shows. Iodine and tellurium isotopes identified are presented in table 1; recommended values for half-lives by Decay Data Evaluation Projet (DDEP) [29] 2 have been considered for all the analyzed nuclides. Although not directly seen in figure 3, cosmogenically produced 22 Na, having a longer mean life than I and Te products, has been identified and quantified using coincidence spectra along data set III.
There are other isotopes, also cosmogenically produced, that have not been identified in this analysis, although they are probably present in ANAIS-25 crystals. It is worth to mention the cases of 129 I and 123 Te.
129 I can be produced by uranium spontaneous fission and by cosmic rays. Its concentration is strongly affected by the ore material exposure either to cosmic rays or to high uranium content environment. It presents 100% β − decay to the excited level of 39.6 keV of the daughter nucleus with a half-life T 1/2 = (16.1 ± 0.7)10 6 y [29], being hence the expected signal in large NaI(Tl) crystals a continuous beta spectrum starting in 39.6 keV. This signal is above the region of interest for dark matter searches, but it is important for background understanding [26]. The long lifetime and the difficulty to disentangle the signal from other  emissions are the reasons why we cannot quantify the amount of 129 I in ANAIS-25 crystals. In [19] the estimated fraction of this isotope was determined to be 129 I/ nat I = (1.7 ± 0.1)10 −13 ; however, strong variability of this concentration is expected in different origin ores.
In the case of 123 Te, we expect the presence of this isotope as daughter of the metastable state 123m Te, effectively identified, and probably also directly produced by cosmic rays by similar reaction mechanisms. This isotope would decay 100% following EC to the ground state of the daughter nucleus, having a half-life larger than 10 13 y [30]. The signature in ANAIS-25 detectors would be a peak corresponding to the full absorption of the binding energy of L-shells of Sb. For the moment, this signal cannot be resolved from other contributions from Te/I isotopes still present, having L-shell EC in higher or lower extent, and then quantification is not possible in ANAIS-25 data.  Table 1. Main cosmogenically induced isotopes identified in the data of ANAIS-25 modules. Halflives and decay mechanisms are indicated. Fourth column shows the energy of main gamma emissions or, in the case of the metastable states, the energy of the corresponding excited level. Last column gives the considered product of intensity I and detection efficiency ǫ for the analyzed distinctive signal of each isotope (see text). Decay information is taken from [29] for 125 I, 127m Te, 123m Te and 22 Na and from [30] for the other isotopes.

Data analysis and estimate of production rates at sea level
We analyze in the following the distinctive signatures considered for the identification of each cosmogenically induced isotope listed in table 1, as well as the estimate of the corresponding detection efficiency ǫ in ANAIS-25 set-up. When required, the efficiencies have been obtained by GEANT4 simulation (version 9.4.p01) considering the complete geometry of the ANAIS-25 set-up and isotopes homogeneously distributed in the crystals; this kind of Monte Carlo simulations was validated with a similar set-up, but slightly different detector geometry, in [26] by reproducing measurements with calibration sources, finding an accuracy at the level of 10%, although a Compton/peak suppression is systematically observed in the simulations. For the metastable isotopes of Te, full energy depositions corresponding to the excited level energy have been considered. Activation has been assessed independently for each crystal, and whenever possible, results have been properly combined.
• For 125 I, the peak at 67.3 keV has been considered. Since this isotope decays with 100% probability by electron capture to an excited state of the Te daughter nucleus, this peak is produced by the sum of the energy of the excited level (35.5 keV) and the binding energy of the K shell of Te (31.8 keV). Efficiency for the detection of the energy deposition is expected to be ∼100%, hence the product Iǫ is taken to be equal to the K-shell EC probability, Iǫ =0.8011 [29].
• For 127m Te, the peak at 88.3 keV, corresponding to the energy of the metastable state, has been analyzed. The branching ratio of the internal transition is I =97.27% [29] and 100% detection efficiency has been assumed. Signal from D1 detector could not be evaluated along data set II since a worse resolution hindered good fitting and therefore in this case final results come only from D0 detector.
• For 125m Te, the peak at 144.8 keV has been studied. Internal transition occurs with 100% probability and detection efficiency has been evaluated by simulation giving ǫ =0.968. In this case, it has been possible to quantify the peak from both the low and high energy spectra, although not for all data sets and detectors; therefore, final results have been obtained only from D0 in the low energy spectrum.
• For 123m Te, the peak at 247.6 keV in high energy spectra has been analyzed. The branching ratio of internal transition is 100%. The loss of efficiency due to the emission with 84% probability of a 159 keV photon [29] which can escape detection has been evaluated by GEANT4 simulation giving an overall detection efficiency Iǫ =0.972.
• Similarly, for 121m Te, the peak at 294.0 keV in high energy spectra has been analyzed and a correction factor to the detection efficiency has been estimated by means of GEANT4 simulation to account for the probability of losing the 212 keV gamma ray which is emitted in 81% of the decays [30]. The branching ratio of internal transition is here 88.6% [30] and then the final value considered for Iǫ is 0.842.
• A coincidence technique analysis has been applied to quantify 126 I and 121 Te, which decay by electron capture to excited states of the daughter nuclei. The identifying signature is the number of coincidences produced by the full absorption of high energy gamma emissions (see table 1) in one detector and that corresponding to the binding energy for K shell electrons in Te or Sb (31.8 or 30.5 keV, respectively) in the other detector. Just as an example, figure 4 shows the coincidence plot between detectors D0 and D1, considering their high and low energy spectra respectively, obtained in the beginning of data set I. By selecting a window from 15 to 45 keV in D1 data, we obtain the D0 coincident spectrum shown in figure 4: peaks from 121 Te and 126 I are clearly observed and are used to determine the number of events attributable to each isotope in a given time period. The product Iǫ for these coincidence signals has been evaluated by simulating the decays of each isotope in the crystals with GEANT4, obtaining Iǫ =0.0035 for 126 I and Iǫ =0.011 for 121 Te. These two isotopes have the shortest half-lives of all the analyzed products. The coincidence signal has been properly identified only within the first 42 days of available data for 126 I. However, for 121 Te the signal has not faded out since this isotope is being produced as result of the decay of the longer lived 121m Te, also cosmogenically induced; therefore, this point has been considered in the analysis of the time evolution of the coincidence signature of 121 Te (see section 3.1.2).
• Presence of 22 Na, having a longer lifetime, is more subtle. The analyzed signal has been the integral number of events from 2300 to 2900 keV (corresponding to full absorption of its positron and gamma emissions) in the spectrum obtained summing D0 and D1 energies, for coincidence events leaving 511 (or 1275) keV at any detector. This is the best signature found for this isotope to avoid the interference of backgrounds. The product Iǫ for these signals has been evaluated by GEANT4 simulation of 22 Na decays in the crystals, obtaining Iǫ =0.0050 for the 511 keV peak and Iǫ =0.0038 for the 1275 keV signal.
In order to quantify the cosmogenic production of every isotope identified in ANAIS-25 measurements, the cosmogenically induced activity, which is the initial activity A 0 corresponding to the moment of storing crystals deep underground at LSC, has been first deduced (section 3.1). For that, the time evolution of the identifying signature produced by each I and Te isotope has been studied considering data sets I and II; as already commented, special analysis is required for 121 Te. For 22 Na, the integrated signal along data set III has been  . Top: Coincidence plot between D0 and D1 in the beginning of data set I. Coincidences due to both 126 I and 121 Te between high energy photons absorbed in D0 and the binding energies of Te or Sb registered in D1 are clearly seen, for K shell (around 31 keV) and even for L shell (around 5 keV). D0 energies below 80 keV have been cut because the high rate of coincidences there prevents from seeing 126 I and 121 Te coincidences. Bottom: Spectrum of the high energy D0 data in coincidence with a window in D1 data around K shell binding energies of Te and Sb; peaks at 507.6-573.1 and 666.3 keV are singled out and can be quantified.
used. Finally, an attempt has been made to estimate production rates R p at sea level from those A 0 activities (section 3.2), taking into account the history of crystals.

I and metastable Te isotopes
The evolution in time of the counting rates R for the different identifying signatures described above has been analyzed and is shown in the plots of figure 5 for I and metastable Te isotopes 3 . Single exponential decays following the radioactive decay law are clearly observed in all of them. Fits of R vs time have been carried out considering not only the sum of the signals from each detector (shown for most of the products) but also independent signals from D0 and D1, and for 125m Te, even signals from the low and high energy spectra have been independently analyzed; in all cases compatible results have been obtained. From the fit parameters the half-lives T 1/2 and initial activities underground A 0 have been derived for each isotope, taking into account the described intensity I and efficiency ǫ for the considered signals, and are summarized in table 2. Error bars plotted in figure 5 represent for Y-axis the uncertainty of peak areas, which has been taken into account in the fits; for X-axis, they only indicate the time range signals have been evaluated over. As shown in table 2, the half-lives deduced from the fits agree with the known values within uncertainties, confirming the choice of signatures for each isotope.
The initial activities A 0 have been also deduced following a different approach, taking into consideration the decay of isotopes during the time intervals over which the identifying signals are integrated. For an isotope with decay time constant λ, the expected number of events C from time t i to t f is: An average of the A 0 values deduced from eq. (3.1) for each period of time has been obtained for each isotope, weighted with uncertainties coming from peak area estimates; results are shown in table 2 too. The two methods give compatible results, although slightly lower values are systematically obtained from the second one; only a larger difference is found for 126 I due to its shorter half-life. The advantage of deriving A 0 by fitting ln R vs time is that no assumption on the nature of the isotope is made, while the value from the average of the initial activities obtained from time intervals takes into account the decay of isotopes during those intervals, which is not completely negligible. Figure 6 shows the evolution in time of the activity of 121 Te deduced from its identifying signal adding D0 and D1 data. A single exponential decay with the half-life of 121 Te is not observed in this case because the isotope is being produced by the decay of 121m Te, also cosmogenically induced and having a longer half-life. A different fitting function is needed to derive the initial activity underground and to estimate the production rate at sea level for 121 Te. The time dependence of the activity of 121 Te during the data taking in ANAIS-25 can be written as the addition of two terms, corresponding to the decay of the cosmogenically produced 121 Te and to the production, and subsequent decay, of 121 Te from 121m Te:

121 Te
with subindex g(m) referring to the ground (metastable) states of the isotope and being N 0 the initial number of nuclei (therefore those cosmogenically induced) and P IT =0.886 the probability of internal transition of the metastable state. To obtain the initial number of nuclei, and then the corresponding activities, for both the ground and metastable states, several fits to the dependence given by eq. (3.2) for the signal from 121 Te have been tried and are also plotted in figure 6. One fit has been performed letting four parameters free: the activities A 0,g = N 0,g λ g and A 0,m = N 0,m λ m and the half-lives T 1/2,g and T 1/2,m . Other fit 122±27 10.5±5.4 2 parameters 59.3±3.9 134±13 --1 parameter -110±12 -- Table 3. Initial activities underground and half-lives derived for the metastable and ground states of 121 Te by fitting the signal plotted in figure 6 to the dependence given by eq. (3.2) with different sets of free parameters (see text).
has been made with two free parameters, fixing the two half-lives according to the known values reported in table 1. The last fit has been carried out with only one free parameter, fixing in addition the activity of the metastable state to the average of the two estimates given in Te activity obtained from the 1-parameter fit (which takes into account all known information) will be used to estimate the production rate.

22 Na
For 22 Na, having a much longer mean life than I and Te isotopes, the evolution in time of the identifying signal has not been studied; signal integrated along the available data in set III (88.12 days) has been considered instead. Figure 7 (top) shows the energy spectra of D0 and D1 registered for coincidences between them for data set III. Two complementary signatures have been considered: events selected by a 511 keV (1275 keV) energy deposition in one detector (fixing a window at ±1.4σ) have been used to build total deposited energy spectra and then, only the region between 2300 and  Figure 7. Top: energy spectra of D0 (black squares) and D1 (red circles) registered for coincidences between them for data set III. Bottom: spectrum of the summed energy of D0 and D1 for events in a ±1.4σ window around the 511 (left) and 1275 keV (right) peaks in any of the coincidence spectra. 22 Na signature has been evaluated from the region between 2300 and 2900 keV (corresponding to full absorption of its positron and gamma emissions). Measured spectra are compared with simulations normalized with the deduced 22 Na activity. 2900 keV (corresponding to full absorption of its positron and gamma emissions) has been used for the rate estimate. Figure 7 (bottom) presents these total energy spectra (summing D0 and D1 energy) obtained for 511 or 1275 keV peak. The effect of coincidences not due to 22 Na has been taken into account subtracting their contribution in the 2300-2900 keV region, evaluated from events in a window on the left of the 511 keV peak and on the right of the 1274 keV peak. From the net number of events registered in the 2300-2900 keV region region and the estimated Iǫ value reported at table 1 the activity of 22 Na at measuring time has been evaluated and then the initial activity underground A 0 deduced, using its half-life from table 1. Table 4 presents the results, for both 511 and 1275 keV windows; since they are compatible, they have been properly averaged. In figure 7 (bottom) the measured spectra of the summed energy of D0 and D1 are compared with the corresponding simulations normalized with the deduced 22 Na activity; the spectra shapes follow reasonably the expectation, somehow validating the method applied to single out 22 Table 4. Measured activity of 22 Na from data set III and the deduced initial activity underground A 0 , considering the identifying signal for coincident events leaving 511 or 1275 keV at one detector (see text). These activities independently estimated have been properly averaged to derive the final result.

Estimate of production rates at sea level
ANAIS-25 crystals were built in Alpha Spectra facilities in Grand Junction, Colorado. The history of the starting material used to produce the NaI powder from which crystals are grown is unknown; this exposure along geological time scales must be relevant only for very longlived isotopes like 129 I or 123 Te, for those having shorter lifetimes saturation would have been probably reached independently from previous exposure history while purifying, growing the crystals and building the detectors (except maybe for 22 Na). Colorado is placed at a quite high altitude; consequently, a very important correction factor f to the cosmic neutron flux at New York City coordinates is reported at [31] for Colorado locations due to altitude and geomagnetic rigidity (4.11 and even 12.86 for Denver and Leadville, respectively). Afterwards, detectors were transported to LSC by boat and by road, being the travel duration about one month, most of the time at sea level. Therefore, the initial activities underground derived in section 3.1 cannot be directly considered the saturation activities (i.e. production rates, R p ) at sea level. However, in the following a method to estimate the production rates in our crystals, based on a few reasonable assumptions, is proposed. Considering that saturation activity is reached at a given place characterized by a correction factor f to the cosmic neutron flux at sea level, the number of nuclei cosmogenically induced for a particular isotope is N sat = f R p /λ. If the material is then exposed to cosmic rays at sea level for a time t, it is straightforward to derive the posterior evolution of the corresponding activity by solving the usual differential equation with simultaneous production and decay of a nuclear species for an initial number of nuclei N sat : The production rate R p of an isotope can be obtained from eq. (3.3) using its initial activity underground, obtained in section 3.1, and fixing the correction factor f and exposure time at seal level t. Since the geographic coordinates of Grand Junction are similar to those of Denver and Leadville, their correction factors f have been adjusted to Grand Junction altitude (4593 ft=1397 m) following the method proposed in [31] (see table 5); the average of both estimates has been taken hereafter as f value. Table 6 summarizes the production rates obtained in this way for f =3.6±0.1, t =30±5 days, λ from T 1/2 values in table 1 and the A 0 values presented in tables 2 (considering the isotope decay along the data taking) and 4. Reported uncertainties include contributions from those derived for A 0 and also from λ, f and t. The relative uncertainty quantified in the production rate is ≤5% (except for 126 I); assuming much larger uncertainties in f and t, as for instance f =3.6±0.5 and t =30±10 days, increases the relative uncertainties in the production rates to ∼15%. Contribution from uncertainty in half-lives is negligible.  Table 5. Estimate of the correction factor f for the cosmic neutron flux at Grand Junction from the factors at Denver and Leadville. The ratio between the flux at these locations and Grand Junction due to different altitude has been derived following [31]; the altitude of Grand Junction is 4593 ft (1397 m). For 121 Te, assuming an exposure to cosmic rays in Colorado until reaching saturation followed by a short exposure at sea level (as made for the other cosmogenic products), the estimate of the production rate from the measured initial activity underground is not so easy due to the combination of direct production and production through the metastable state. In the following, R p,g and R p,m correspond to the production rates at sea level of the ground and metastable states of 121 Te.
At an altitude with a correction factor f for the cosmic neutron flux, the number of nuclei N of each species is given by the system of coupled equations: The system in eq. (3.4) can be solved assuming that initially N g = N m = 0. Then, saturation values (for t → ∞) are: At sea level, the system of coupled equations giving the number of nuclei of ground and metastable isotopes is like that in eq. (3.4) with f =1 and can be solved assuming initially the saturation values from eq. (3.5). Then, the activity for the metastable state is given by eq. (3.3), whereas activity for the ground state is obtained as: λ g λ g − λ m (3.6) Table 6 presents the production rate at sea level R p,g , derived from eq. (3.6), using the corresponding initial activity underground deduced in table 3 from the 1-parameter fit, the values of λ m , λ g and P IT and the values of f and t already assumed for the other cosmogenic products. Uncertainty in the activity give the main contribution to the uncertainty in the production rate.

Calculation of production rates at sea level
As pointed out in section 1, cosmogenic activation of materials is an hazard for experiments looking for rare events and therefore it must be properly assessed in the background models of this type of experiments. Production rates are the main input to compute the induced activities of the relevant isotopes. Since direct estimates of the production rates, as those presented here for NaI, are not available for all materials and products, in many cases production rates must be evaluated using eq. (4.1) from the cosmic ray flux φ and the production cross sections σ, both dependent on the energy E of cosmic particles: This calculation method is strongly dependent on the excitation functions and assumed cosmic ray spectrum. In order to validate different approaches to such a calculation, production rates at sea level for all the cosmogenic products identified and quantified in ANAIS-25 crystals will be evaluated and then compared with the values obtained from the measurements presented in section 3. Although not observed in ANAIS-25 data mainly due to its short half-life (T 1/2 =4.1760±0.0003 days), 124 I has been considered in this study too. Production rate of 123 Te, having low energy emissions dangerous for dark matter searches, has been also evaluated. As proposed in [5], first, all the available information on the relevant production cross sections by neutrons and protons from both measurements and calculations will be presented; then, the best description of the excitation functions of products over the whole energy range will be chosen (section 4.1). Finally, production rates considering a particular cosmic ray spectrum (section 4.2) and the selected excitation functions will be calculated and compared with the experimental values (section 4.3).

Excitation functions
The excitation function for the production of a certain isotope by nucleons in a target over a wide range of energies (from some MeV up to several GeV) can be hardly obtained only from experimental data, since the measurements of production cross-sections with beams are demanding. The use of computational codes is therefore mandatory, although experimental data are essential to check the reliability of calculations. The difficulties for properly describing excitation functions are discussed for instance at [32].
Most of the considered cosmogenic isotopes are induced by interaction on 127 I; only 22 Na is produced on 23 Na. Both 127 I and 23 Na have a 100% isotopic abundance in natural iodine and sodium respectively. Figures 9-16 in the appendix show the information compiled for the excitation functions of the cosmogenic products induced on 127 I or 23 Na, taken from different sources: • The EXFOR database (CSISRS in US) [33] provided most of the experimental data.
• Libraries MENDL-2 and MENDL-2P (Medium Energy Nuclear Data Library) [37] gave results obtained by Monte Carlo simulation of the hadronic interaction between nucleons and nuclei using codes of the ALICE family, for neutrons up to 100 MeV and protons up to 200 MeV.
• At library TENDL-2013 (TALYS-based Evaluated Nuclear Data Library) [38], cross sections obtained with the TALYS nuclear model code system were found for neutrons and protons up to 200 MeV. An outstanding feature of this library is that production of metastable and ground states is evaluated separately.
• The YIELDX routine was used to compute production cross section at the whole energy range. This routine implements the most updated version of the popular Silberberg & Tsao semiempirical formulae [39]- [41] giving nucleon-nucleus cross sections for different reactions. They can be used for a wide range of targets and products at energies above 100 MeV; since they are based only on proton-induced reactions, cross section must be assumed to be the same for neutrons and protons.
• Library HEAD-2009 (High Energy Activation Data) [42] merges data for neutrons and protons up to 1 GeV. Here, data from HEPAD-2008 (High-Energy Proton Activation Data) sub-library have been considered, obtained using a selection of models and codes (CEM03.01, CASCADE/INPE, MCNPX 2.6.0,. . . ) dictated by an extensive comparison with EXFOR data.
At low energies (LE, below 200 MeV) excitation functions for protons and neutrons are clearly different, from both MENDL and TENDL-2013 data; therefore, it is very important to quantify cosmogenic yields considering cross sections specifically evaluated for neutrons. For I isotopes, experimental data on cross sections are abundant and are in general in reasonable agreement with the excitation functions provided by MENDL, TENDL-2013, YIELDX and HEAD-2009. However, for Te isotopes, experimental information is scarce. It is worth noting that some measurements distinguish production cross sections for ground (G) or metastable (M) states of Te isotopes, but this is not the case for the considered calculations except for TENDL-2013.
To quantify the agreement between the available experimental data and the different calculations helps to choose the best description of the excitation functions. For that, deviation factors F have been evaluated following a definition typically used in the literature: being N the number of pairs of experimental and calculated cross sections σ exp,i and σ cal,i at the same energies. Parameter d is the mean square logarithmic deviation.  found for 127m Te and HEAD-2009 calculation, predicting very low production cross sections well below available measurements for both metastable and ground states. Calculations at high energy show slightly better deviation factors than at low enegies; in the overall results not very significant differences are found between HEAD-2009 and YIELDX nor between MENLD-2 and TENDL-2013.
Different descriptions of the excitation functions have been analyzed for the calculation of production rates, based on the available data and calculations at the low and high energy ranges: • For iodine isotopes, thanks to the large amount of experimental data available, the description could be fixed for each isotope following the lowest deviation factors at low and high energies according to table 7.
• For tellurium isotopes, TENDL-2013 data for neutrons were always considered below 200 MeV, since they have been calculated separately for metastable states; HEAD-2009 data, having better deviation factors than YIELDX, were used above that energy for most of the isotopes. For 127m Te, it was preferred to extrapolate the TENDL-2013 value at 200 MeV to higher energies. The branching ratios for the production at ground or metastable state given by TENDL-2013 data at 200 MeV (∼0.7 for metastable production for all the considered isotopes) have been applied to data at higher energies, which do not distinguish between ground or metastable final states.
• Regarding 22 Na, the only description possible was using TENDL-2013 at low energies and YIELDX results at high energies.

Cosmic rays spectrum
Following eq. (4.1), the other ingredient required to evaluate production rates, together with the excitation function, is the cosmic ray spectrum. Different forms of the neutron spectrum at sea level have been used in cosmogenic activation studies, like the historical Hess spectrum Energy (MeV) Differential neutron flux (n/cm2/s/MeV) Figure 8. Energy spectrum of cosmic neutrons following the analytic expression in [44] (for conditions of New York City at sea level) which has been considered in this work to calculate production rates of cosmogenic isotopes. [43], that proposed by Lal & Peters [1] or the parametrization based on a compilation of measurements presented by Ziegler [31]. In [44], a set of measurements of cosmic neutrons on the ground across the US was accomplished using Bonner sphere spectrometers; an analytic expression fitting data for energies above 0.4 MeV was proposed and has been considered in this work. Figure 8 shows the energy spectrum given by this description of cosmic neutrons; applied to conditions of New York City at sea level, the integral flux from 10 MeV to 10 GeV is 3.6×10 −3 cm −2 s −1 . Table 8 presents the obtained production rates for the analyzed cosmogenic isotopes considering the selected excitation functions. Results are given per day and per kg of NaI. Contributions from the LE and HE descriptions of the cross sections are shown together with the total rate. For most of the isotopes, the dominant contribution comes from the low energy region; for instance, for iodine isotopes it accounts for ∼90% of the rate. Table 8 shows also the production rates deduced experimentally in section 3.2 for comparison. Ratios between the calculated and experimental rates are given too. Most of the calculated production rates are larger than the experimental ones. Good agreement has been found for I isotopes, having cross sections well validated not only at high energy with protons but also at low energies with neutrons. The largest deviation is obtained for 125m Te; the lack of experimental data in the whole energy range (see Fig. 13) makes difficult try to improve the selection of the excitation function. For 22 Na, there is an important step in the description of the excitation curve at LE and HE regions (see Fig. 16) and then contribution from HE could be overestimated and/or that from LE undervalued, justifying the difference between the calculated and experimental rates; in any case, they are compatible with the calculated maximum rate of ≃ 100 kg −1 d −1 at sea level quoted in [19]. The fact of observing a rate below the calculation could indicate that saturation was indeed not reached for this isotope. There is a similar mismatch in the description of the excitation function at LE and HE regions for 121 Te (see Fig. 15); but in this case re-scaling the curves would give production rates still compatible with the observed value. Given the discrepancies found between calculated and experimental production rates, it can be concluded that calculations of activation yields in sodium iodide for ultra-low background experiments are precise enough, provided the used cross sections are well validated along all the relevant energy range.

Conclusions
Cosmogenic activation in NaI(Tl) crystals has been evaluated from the measurements taken at LSC using two 12.5 kg NaI(Tl) detectors produced by Alpha Spectra, in Colorado (US), and then shipped to Spain. The prompt data taking starting, the low radioactive background conditions and the very good energy resolution of the detectors made possible a reliable quantification of production of isotopes with half-lives larger than ten days.
Signatures from 125 I, 127m Te, 125m Te, 123m Te and 121m Te were directly identified in the energy spectra registered by the detectors while production of 126 I, 121 Te and 22 Na was revealed by analyzing coincidences between both detectors. The evolution of the identifying signature of all these isotopes, except 22 Na, for a period of 210 days, following the expected exponential decay, allowed to confirm the hypothesized origin of the cosmogenic signals and to evaluate the activities at the start of the underground storage. Signal from 22 Na was obtained from spectra taken after fifteen months. Independent analysis performed for each detector or different data sets and following different methodology gave compatible results within uncertainties. Tables 2-4 summarize the measured initial activity underground for all the observed cosmogenic products. Production rates at sea level were estimated from these activities considering that saturation was reached during detector production in Colorado (assuming a cosmic neutron flux 3.6 times higher than at sea level) and that detectors were further exposed for 30 days during the trip by boat from US to Spain. The hypothesis of saturation is reasonable for the short-lived isotopes analyzed. Results are presented in table 6. The analysis of the evolution in time of 121 Te signal was not so easy because direct production and production from the metastable state, also cosmogenically induced with a longer life time, were combined. Values for the metastable isotope deduced in this analysis were compatible with those obtained directly from the 121m Te signature.
Calculations of production rates for all the cosmogenic isotopes observed were carried out considering the cosmic neutron spectrum given in [44] and selecting carefully the excitation functions. For each product, the available information on production cross sections by neutrons and protons was collected, considering measurements, results compiled in different libraries and calculations from semiempirical formulae made on purpose. Data from MENDL-2 or TENDL-2013 for neutrons were considered below 100 or 200 MeV while at higher energies YIELDX calculations or results from the HEAD-2009 library were used. For the analyzed isotopes, the overall deviation factors range from 1.4 to 1.8, being slightly better at high energies. Calculated production rates are summarized in table 8; agreement with the rates experimentally deduced is very good for some of the products, having excitation functions well validated against measurements. This comparison allows to quantify the systematic uncertainty when applying the followed methodology to the study of other cosmogenic products. Uncertainties in this kind of calculations come mainly from the difficulties encountered on the accurate description of cosmic ray spectra, on the precise evaluation of inclusive production cross-sections in all the relevant energy range and on the good knowledge of the material history; the low and medium energy regions below a few hundreds of MeV are the most problematic ones since neutron data are scarce and differences between neutron and proton cross sections may be important.
Some of the cosmogenic isotopes observed in this study generate emissions which could affect the very low energy region where dark matter signal is expected to appear in NaI(Tl) detectors. The effect of long-lived nuclides like 123 Te or 22 Na must be considered. Below 10 keV, the estimated contribution of 22 Na, based on GEANT4 simulation for the ANAIS-25 set-up, is of 3.8 kg −1 d −1 peaked around 0.9 keV for the activity deduced from data set III; as a reference, contribution from the measured level of 40 K in ANAIS-25 crystals is 5.0 kg −1 d −1 peaked at 3.2 keV. Concerning 123 Te, considering the production rates deduced in table 8 for both metastable and ground states, for an exposure lower than 10 10 y, present activity should be lower than 0.03 kg −1 d −1 . The other products have mean lives low enough to decay in reasonable times and they will not be a problem for dark matter detection provided crystals are left to cool underground before starting data taking.