# Radioactive decay of <sup>90</sup>Sr in cement: a non-equilibrium first-principles investigation

Jorge Kohanoff, Alfredo Correa, Gleb Gribakin, Conrad Johnston, Andrés Saúl

► To cite this version:

Jorge Kohanoff, Alfredo Correa, Gleb Gribakin, Conrad Johnston, Andrés Saúl. Radioactive decay of <sup>90</sup>Sr in cement: a non-equilibrium first-principles investigation. The European Physical Journal D: Atomic, molecular, optical and plasma physics, 2021, 75 (9), 10.1140/epjd/s10053-021-00202-8. hal-03358070

HAL Id: hal-03358070

https://amu.hal.science/hal-03358070

Submitted on 4 Oct 2021

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Radioactive decay of <sup>90</sup>Sr in cement: a non-equilibrium first-principles investigation

Jorge Kohanoff<fnref n="1" />,<fnref n="2" />,ª, Alfredo A. Correa<fnref n="3" />, Gleb Gribakin<fnref n="4" />, Conrad Johnston<fnref n="1" />, and Andrés Saúl<fnref n="5" />,ᵇ

Received 31 March 2021 / Accepted 14 June 2021
© The Author(s) 2021

**Abstract.** Cement is an inexpensive and relatively easily manageable material that is used as a last barrier for nuclear waste disposal. Under these conditions, the cement is in contact with low radiation doses, but there is a distinct possibility of being contaminated with radioactive products. Of particular concern is the medium lived half-life product <sup>90</sup>Sr (28.8 years) due to its ability to replace Ca. <sup>90</sup>Sr undergoes β-decay to <sup>90</sup>Y which, in turn, β-decays to stable <sup>90</sup>Zr. In this work, we discuss systematically the chain of non-equilibrium processes that result as a consequence of β-decay events in cement. We first use density functional-based methods to study the consequences of the sudden increase of the nuclear charge from Z to Z + 1, a possible induced ionization and the perturbation of the surrounding electronic charge. Secondly, we use molecular dynamics simulations to study the recoil of the daughter nucleus. Finally, we discuss the damage caused by the ionization cascade produced during the propagation of the β-electron and the resulting chemical and structural perturbation.

## 1 Introduction

Nuclear fission processes generate a large variety of radioactive products of different lifetimes. Those with a short half-life, e.g., months to a few years, are relevant during the period of operation of the power plant and, after decommissioning, are managed simply by storing the used *hot* fuel for a few years until it can be recycled or buried away. Pu and other very long-lived transuranic elements—over tens of thousands of years—are a major concern for long-term nuclear waste technologies. Some countries have chosen the road of sealing the used fuel components in steel canisters and burying them in geologically stable areas [1]. Others are looking for stable mineral forms that can hold and stand the decay of radioactive atoms, after separating them via chemical processes [2].

A particularly challenging problem is posed by intermediate half-life fission products like <sup>90</sup>Sr (28.8 years) and <sup>137</sup>Cs (30 years). Under accidental conditions, these

isotopes can escape into the environment and cause health problems in the population by entering the food chain. This is particularly relevant in the case of Sr that is isovalent with Ca and has been shown to replace it in bones and teeth [3]. This ability to replace Ca manifests itself also in geological and building materials like clays and cements, which are almost invariably used as the last barrier for nuclear waste disposal, and are in contact with the terrain and exposed to ground waters and radiation. In the case of accidents like Chernobyl, radioactive Sr and Cs have been found in the concrete walls of surrounding buildings. It is then of great interest to understand what are the effects of the radioactive decay of these species in cement.

In this work, we will focus on the behavior of the <sup>90</sup>Sr substitution in cement. This species undergoes β-decay to <sup>90</sup>Y with a half-life of 28.8 years which, in turn, β-decays to stable <sup>90</sup>Zr with a half-life of 64 hours. There are several aspects of β-decay that are worth considering. First, there is the transmutation in atomic number from Z to Z + 1, which induces electronic and structural rearrangements. In the isolated atom, this electronic relaxation can lead to a sequence of ionization events, known as shake-off, a process that has been studied by Carlson et al. for each of the atomic shells for all the elements from Z = 2 to Z = 92 [4]. Contrary to all previous calculations of electron shake-off, the one presented here takes into account that the process occurs in

<footnote n="1">¹ Atomistic Simulation Centre, School of Mathematics and Physics, Queen's University Belfast, Belfast BT7 1NN, UK</footnote>

<footnote n="2">² Instituto de Fusión Nuclear “Guillermo Velarde”, Universidad Politécnica de Madrid, c/José Gutierrez Abascal 2, 28006 Madrid, Spain</footnote>

<footnote n="3">³ Lawrence Livermore National Laboratory, Livermore, CA 94588, USA</footnote>

<footnote n="4">⁴ Center for Theoretical Atomic and Optical Physics, Queen's University Belfast, University Road, Belfast BT9 1NN, UK</footnote>

<footnote n="5">⁵ Centre Interdisciplinaire de Nanoscience de Marseille-CNRS (UMR 7325), Aix-Marseille University, Marseille, France</footnote>

<footnote>Advances in Multi-Scale Modelling of Intense Electronic Excitation Processes Guest editors: Jorge Kohanoff, Antonio Rivera, Eduardo Oliva Gonzalo, Andrey V. Solov'yov and Tzveta Apostolova.</footnote>

<footnote>ªe-mail: j.kohanoff@upm.es</footnote>

<footnote>ᵇe-mail: saul@cinam.univ-mrs.fr (corresponding author)</footnote>
a bulk (solid) system. Secondly, the emitted β-electron typically travels a few mm producing a complex cascade in which it generates secondary electrons, X-rays, and Auger electrons. By momentum conservation, the daughter nucleus experiences a recoil that, in principle, can also result in structural consequences.

In this article, we will study the consequences of each of these effects in cement using adiabatic and non-adiabatic first-principles simulations. In Sect. 2, we present the atomic structure of the model cement paste and the details of the different density functional theory (DFT)-based calculation methods used in this work. In Sect. 3.1, we discuss the consequences of the sudden increase of the nuclear charge (Z) due to β-decay: a possible induced ionization (Sect. 3.1.1) and the perturbation of the surrounding electronic charge (Sect. 3.1.2). In Sect. 3.2 we present the results of molecular dynamics simulations (MD) of the recoil of the daughter nucleus. In Sect. 3.3 we discuss the damage caused by the ionization cascade produced during the propagation of the emitted electron (Sect. 3.3.1) and the resulting chemical and structural perturbations (Sect. 3.3.2). Section 4 is devoted to the discussion and conclusions.

# 2 Calculation details

## 2.1 The structure of the model cement system

C-S-H (calcium-silicate-hydrate) is the principal binding product of cement hydration. The atomic structure for C-S-H used in this work is presented in Fig. 1. It was generated following the procedure reported in Ref. [5] and used in a previous paper [6]. The cell contains 501 atoms and has the chemical composition Ca<sub>72</sub>Si<sub>44</sub>O<sub>235</sub>H<sub>150</sub>, with a C/S ratio of 1.64 and a density of 2.33 g/cm<sup>3</sup>, in very good agreement with experimental data. Figure 1 shows the silicate chains formed by corner sharing SiO<sub>4</sub><sup>-</sup> tetrahedra. It can also be seen that the Ca<sup>2+</sup> ions appear in two distinct environments: intralayer and interlayer. Water, in the form of H<sub>2</sub>O molecules and OH<sup>-</sup> groups, is present in the interlayer space.

## 2.2 Density functional theory-based methods

To simulate the various aspects of the problem, we have used a variety of methods and codes, although uniformly based in DFT. The choice of code has been dictated by their efficiency for the specific task. They are described in this section. Here, it is important to remark that after more than 20 years of development, testing, and collaboration, the majority of the DFT codes available nowadays, when used correctly, give consistent results [7].

In order to study the successive evolution, we need to follow the electronic dynamics (and later atomic dynamics) starting from meaningful initial conditions.

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Schematic representation of the C-S-H structure. Hydrogen atoms are in white, oxygen atoms in red, calcium in yellow, and Si in blue. The corner sharing SiO<sub>4</sub><sup>-</sup> tetrahedra forming the silicate chains are also shown. The interlayer Ca atom (marked as ‘Ca56’) is substituted by Sr, which later decays into Y and subsequently into Zr. The marked atom (Ca56) is in the initial position of the Zr and Y atoms for the recoil simulations</figcaption></figure>

Starting from the structure representative of cement shown in Fig. 1, we replace a Ca with a Sr atom and relax the atomic structure using the same DFT-GGA approach as for the pure cement. Except serving as an initial condition for the rest of the work, we are not interested in how an individual Sr replaced the original Ca, which can take years as a result of diffusive processes.

### 2.2.1 Time dependent density functional theory calculations

To represent an instantaneous transmutation, at some arbitrary point in time (that we define as $t = 0$), we replace the Sr nucleus with an Y without modifying the electronic state (and electronic density). This is an application of the sudden approximation in quantum mechanics [8], where the electronic state is continuous in spite of the time discontinuity of the external potential. Therefore, the electronic state finds itself away from the ground state of the new nuclear potential (with the Y nucleus). It is the dynamics of the electronic state (and electronic density) on the femtosecond time scale that concerns us in this section.

The Ca<sub>71</sub>YSi<sub>44</sub>O<sub>235</sub>H<sub>150</sub> system, which is in an excited electronic state, is then allowed to evolve via real-time time-dependent density functional (rt-TDDFT) description of the electron dynamics combined with the Ehrenfest approximation for the calcula-
tion of the forces. This has been done using the QB@LL code [9], which describes the electronic Kohn–Sham orbitals and electronic density in a plane wave basis set, and replaces the core electrons with norm-conserving pseudopotentials. The time-dependent Kohn–Sham equations were integrated using a Runge–Kutta integrator with a time step of 1 attosecond. The simulations ran for 1 fs; therefore, the nuclei did not have sufficient time to react and move.

### 2.2.2 First principles molecular dynamics simulations

In order to investigate the structural relaxation induced by the recoil of the daughter nuclei, we performed first-principles molecular dynamics simulations within the DFT framework using the Quantum-espresso plane-wave code [10]. We used the Perdew–Burke–Ernzerhof generalized gradient approximation (PBE-GGA) [11] and ultrasoft pseudopotentials with explicit semi-core electrons. The plane-wave cutoff for the Kohn–Sham orbitals was set to 80 Ry, and all the calculations were performed at constant volume with a single k-point (the Γ-point).

### 2.2.3 Excess holes and electrons

The passage of the β-electron through the cement produces electronic excitation and ionization of the target, leaving holes in certain regions of the sample, and excess electrons in others. To find out where are these localized, we conducted ground state electronic structure calculations with an excess electron and an excess hole in the sample shown in Fig. 1, i.e., without Sr substitution. This is a valid procedure because the concentration of Sr is relative small (maximum 4%, given by the solubility limit [6]) and the chances that charge defects localize next to another Sr defect are meagre. The calculations were carried out using the GPW method embodied into the CP2K package [12]. In this method, the Kohn–Sham orbitals are expanded in a Gaussian basis set while the electronic density is represented by plane waves. The core electrons are replaced by Goedecker–Teter–Hutter (GTH) pseudopotentials [13], and the basis set used was TZV2P-MOLOPT-GTH for all atoms except Ca, for which a DZVP-MOLOPT-GTH was adopted [14]. These are a triple-zeta basis with two sets of polarization functions, and a double-zeta polarized basis set finely tuned for molecular systems. The charge of the electron or hole was compensated with a uniform background of opposite charge, so that the unit cell remained neutral. While previous work used the same GGA-PBE functional but different basis sets and pseudopotentials [6], to be consistent, we decided to maintain the same atomic structure without re-optimizing the geometry. We nevertheless checked that this choice did not lead to any consequences for the location of excess electrons and holes. Since GGAs are known to suffer from the self-interaction problem, which is very important for unpaired electrons and holes, we re-calculated the electronic ground state for the same

atomic structure, using the hybrid functional M06-2X, which contains 54% of Hartree-Fock exchange [15]. The electronic Brillouin zone was sampled at the Γ-point.

## 3 Results

As mentioned above, the transmutation nuclear reaction is composed of three different phenomena: (1) the increase of the nuclear charge from Z to Z + 1, (2) the emission of a β-electron, and (3) the accompanying recoil of the nucleus. In what follows, these three aspects will be addressed individually. We first study the relaxation following nuclear transmutation. The two aspects, electronic and structural, occur on very different time scales due to the mass difference between electrons and nuclei. Therefore, in a first instance, we follow the electron dynamics at fixed nuclear geometry, and then we study the consequences of the recoil. Finally, we study the effects of the β-electron emission.

### 3.1 Sudden increase of the nuclear charge

#### 3.1.1 Probability of β-induced ionization

The main consequence of the nuclear charge increase is to perturb the electronic subsystem, bringing it into a non-stationary state. This shake-up phenomenon is followed by the relaxation of the electronic component in ways that depend on the specific atomic environment.

In the gas phase, an important relaxation channel following β-decay is the emission of a secondary electron, in a phenomenon usually called shake-off, which has been first studied theoretically in the 1950s [16,17]. In order to understand which atomic shells are more prone to be ionized via shake-off, we performed a Hartree-Fock atomic calculation including core levels. The radial parts of the electronic orbitals of Sr²⁺ and Y³⁺ are shown in Fig. 2 together with the changes upon transmutation. Also reported are the overlap integrals of the radial wavefunctions,

$$w_{nl} = \int_{0}^{\infty} \tilde{P}_{nl}(r) P_{nl}(r) dr, \quad (1)$$

where $\tilde{P}_{nl}(r)$ is the radial wavefunction of Sr²⁺, and $P_{nl}(r)$ is that of Y³⁺.

The probability of ionization for the nℓ subshell is given by:

$$p_{nl} = 1 - (w_{nl}^2)^{N_{nl}} \approx N_{nl}(1 - w_{nl}^2) \quad (2)$$

where $N_{nl} = 2(2\ell + 1)$ is the number of electrons in the nℓ subshell. Note that Eq. (2) is valid only when $p_{nl} \ll 1$, i.e., when the overlaps $w_{nl}$ are close to unity, which is the case in the present calculation. Using the overlaps from Fig. 2, we obtain the probabilities reported in Table 1.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** The radial wavefunctions $\tilde{P}_{nl}(r)$ of Sr$^{2+}$ (solid red) and $P_{nl}(r)$ of Y$^{3+}$ (dashed blue) versus $r$ in [Å]. The numerical value on each panel is the overlap integral $w_{nl}$ between the radial wavefunctions. Solid green curves show $(\tilde{P}_{nl}(r) - P_{nl}(r))P_{nl}(r)$, whose integral from 0 to $\infty$ gives the difference between $w_{nl}$ and 1</figcaption></figure>

The results in Table 1 show that the total probability that transmutation is accompanied by the ejection of a core (subshells 1s to 3d) or a semicore (subshells 4s and 4p), electron is 6.5% and 11.9%, respectively.

### 3.1.2 Transmutation in the condensed phase

Starting from the structure representative of cement shown in Fig. 1, we replaced Ca with a Sr atom and relaxed the structure using *ab initio* DFT-GGA electronic structure calculations. We next replaced the Sr with a Y pseudopotential, *without* modifying the electronic state (and electronic density), to mimic an instantaneous transmutation as discussed above. Both

pseudopotentials include the semicore 4s and 4p electrons in the valence.

Since most of the transition probability is at valence (5s) and semi-core (4s, and 4p), an all-electron calculation is unnecessary and a frozen core approximation is justified even under these conditions, provided that the semi-core electrons are considered explicitly, as these participate in the shake-off. In summary, the use of pseudopotentials is a reasonable approximation when there is a sudden nuclear change from $Z \to Z + 1$. In other words, for core electrons the difference between $Z$ and $Z + 1$ is a relatively small perturbation for $Z = 38$. Within this framework, the transmutation event is represented by the instantaneous replacement of one pseudopotential with another, without relaxing the electronic structure.

The Ca$_{71}$YSi$_{44}$O$_{235}$H$_{150}$ system is initially in an excited electronic state. At this point we start an electron dynamics simulation using rt-TDDFT combined with Ehrenfest dynamics for the calculation of the forces. In Fig. 3, we show the time evolution of the charge in the vicinity of the Sr/Y atoms, calculated using the Bader and Voronoi constructions [18]. It can be seen that initially, in times of the order of 0.2 fs, almost one unit of charge (valence charge from 8.2 to 9) surrounds the Y atom. After 0.5 fs the charge levels off near the ground state value, after the initial overshooting and several oscillations. The reason for this is clear: After transmutation, there is a charge imbalance in the vicinity of the Y atom, which would attract an additional neutralizing electron, should there be one available. While this is not possible in the gas phase, and the energy is partly dissipated in ejecting another electron, in the condensed phase there are plenty of available electrons in the host and there are other channels available for dissipation. This is what is observed in Fig. 3. Note that the core electron charge (not simulated by rt-TDDFT) instead should remain practically unchanged according to the atomic Hartree-Fock calculation (Sect. 3.1.1).

A spatially resolved view of the electron capture can be seen in Fig. 4, which shows a snapshot of the change in the electronic charge distributions relative to the initial one, taken at a time 0.3 fs after transmutation.

The presence of neighboring atoms chemically bound to the Y atom opens a dissipation channel that is inaccessible to the prototypical shake-off of an isolated atom. The charge redistribution generates forces in these atoms that induce their motion. We calculated the time evolution of the force on the neighboring atoms during the electronic evolution (Fig. 5). These are non-adiabatic Ehrenfest forces that include, however, a nonzero adiabatic component due to the replace-

<table><caption>**Table 1** Shake-off ionization probabilities (in percent) for the different electronic shells of Y$^{3+}$ after transmutation from Sr$^{2+}$</caption><thead><tr><th></th><th>1s</th><th>2s</th><th>2p</th><th>3s</th><th>3p</th><th>3d</th><th>4s</th><th>4p</th></tr></thead><tbody><tr><th>p<sub>nl</sub></th><td>0.10</td><td>0.31</td><td>0.57</td><td>0.70</td><td>1.98</td><td>2.84</td><td>2.36</td><td>9.54</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Valence charge around the Sr/Y site during the first femtosecond after the β-decay (transmutation) event. Note the logarithmic scale in time. Bader and Voronoi criterion are used to assign charge to the site [18]. Horizontal lines show initial Sr and final Y to the corresponding ground state values (for fixed ion positions)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Change in electronic charge distribution 0.3 fs after the β-decay (transmutation) event. Blue (red) indicates excess (defect) of electrons relative to the initial state. The neighboring atoms to the Y are indicated in red (O), white (H), and brown (Si). The excess and defect of electrons follow the deformed outer electronic shells, the surrounding O atoms lose a fraction of electronic charge each or polarize</figcaption></figure>

ment of Sr by Y. Initially there is large force due to the unscreened ion-ion repulsion originated in the increased nuclear charge after transmutation. As the charge distribution settles into its new ground state, the forces decrease in magnitude and stabilize at a value of the order of a few eV/Å. These remnant forces are later responsible for the relaxation of the local structure

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Upper panel: structure of first neighbors around Y (before β-decay). Lower panel: forces on the neighboring O-atoms in the direction of the Y–O bonds, as a function of time after β-decay</figcaption></figure>

around the Y-atom in a longer time scale. They are in addition to the sudden forces induced by the recoil of the Y nucleus at the moment of the β ejection.

We will see in the next section that the Y nucleus recoil does not carry any relevant structural effects. Electron-phonon thermalization occurs in the ps–ns timescale, so that by the time the second β-decay from Y to Zr occurs (the lifetime of the Y nucleus is 64 hours), the system is in its new equilibrium state. This second β-decay can be analyzed by the same type of non-adiabatic simulation and should follow an analogous behavior for the density and the forces.

## 3.2 Recoil

Concomitant with the electronic redistribution described in the previous section, the transmuted nucleus experiences a recoil due to momentum conservation. Recoil events have been studied by Marks et al. in the specific cases of SrTiO₃ and SrH₂ [19]. In our case, the maxi-
mum kinetic energy of the recoiling atoms is of the order of 5 eV for the first decay, $^{90}$Sr to $^{90}$Y, and 20 eV for the second one, $^{90}$Y to $^{90}$Zr. These energies have been shown to be enough to produce some structural changes in molecules. In the solid state, the presence of the environment poses a more stringent challenge, with typical displacement thresholds of the order of tens of eV.

Since the outcome of the recoil process is not obvious, we have conducted adiabatic first-principles molecular dynamics simulations where first the Y, and then the Zr, were given the maximum allowed recoil energy in a few randomly selected directions. This was done as explained in Sect. 2.2.2. For the Y recoil, we have not observed any significant structural rearrangements, as the recoil energy is not sufficiently large. As we shall see below, this is not the case of Zr.

The results corresponding to the Zr recoil along the [111] direction are shown in Fig. 6. The time evolution of the kinetic and potential energy of the simulation box is shown in Fig. 6a. The system has been first thermalized during 0.2 ps at 300 K. Then the Zr kinetic

energy has been suddenly increased to 20 eV by giving the Zr-atom a suitable velocity along the desired direction. This can be seen as a jump in the total kinetic energy at $t = 0.2$ ps. The simulation then continued during 1 ps where after a short transient of about 100 fs the excess kinetic energy is spread among all the atoms in the simulation box reaching a final temperature of about 450 K, as seen in Fig. 6a. It is to be noticed that also the potential energy increases in a similar amount, following equipartition. The relative displacement of the Zr-atom is shown in Fig. 6b. As can be seen, only the *x* and *y* coordinates change significantly after the kick, indicating that the interlayer spacing acts like a diffusion channel for the contaminating species. Other simulations with initial velocity having a component along the *z*-direction do not result in the impurity crossing from one channel to another.

## 3.3 Following the β-electron in its tracks

### 3.3.1 Radiation cascade

The β-electron resulting from the decay of $^{90}$Sr to $^{90}$Y is emitted with an energy of 0.5 MeV [19] and propagates through the medium generating a rich cascade of events. Typically, the penetration depth of such an electron is about 0.5 mm, and the structure of the track is not a straight line but a branched, tree-like structure. The final result of the cascade is a distribution of defects, mainly in the form of electrons and holes in displaced positions. Recombination here is different that in the case of optical excitations, because secondary electrons are produced with a kinetic energy that takes them far enough from their respective holes. Moreover, cement is an electronic insulator. Therefore, unlike the case of metals, these defects are relatively stable. The track structure, i.e., the distribution and density of defects, can be calculated via Monte Carlo transport simulations [20–22]. This approach uses the energy distribution of the β-electron, which is obtained from nuclear physics considerations, to generate initial conditions that are then sampled.

### 3.3.2 Chemical effects: electrons and holes

Knowing where exactly in the cement matrix will electrons and holes localize is important to assess the chemical and structural stability of the cement upon irradiation. These species participate in chemical reactions leading, e.g., to the production of $H_2$ from water dissociation [23]. The accumulation of $H_2$ is undesirable, especially when cement is used as a container for nuclear waste decommissioning, as it can build up pressure and lead to fracture and hence leakage of radioactive material into the soil and ground waters. Here we focus on understanding where do electrons and holes localize, which new chemical species are formed, and how these can participate in chemical reactions.

To answer these questions, we carried out electronic structure calculations of excess electrons and holes in

<figure><img src="image_7.png" /><figcaption>Fig. 6 Molecular dynamics simulation of the recoil excerpted on the Zr daughter atom. The calculations have been performed using the Quantum Espresso code [10] with a time step of 1 fs. a Evolution of the kinetic and potential energy of the simulation box. b Relative displacement of the Zr atom with respect to its time equal zero position. c Initial atomic configuration. d Final atomic configuration</figcaption></figure>
the structure shown in Fig. 1, according to the procedure described in Sect. 2.2.3. The addition or removal of an elementary charge from the simulation cell, which is otherwise neutral, leads to a non zero spin density. This is so because the neutral system is spin paired and exhibits an energy gap, and the additional electron occupies the lowest unoccupied molecular orbital in one of the spin projections. The spin density thus corresponds to the electronic density of the excess electron. A similar argument applies to the hole. This is shown by yellow isosurfaces in Fig. 7.

The excess electron in C–S–H is observed to localize in the calcium and water-rich interlayer regions across the water molecules, OH⁻ groups and any nearby interlayer calcium ions. These observations suggest that the structure exerts a minimal influence on the location of the electron. Moreover, first-principles MD simulations show the excess electron moving easily within and between the interlayer spacings along the dynamics. The fact that excess electrons localize in these water-rich regions is important as these species may become mobile and may react with water to form a variety of radicals and other species such as H• and H₂ [24]. In fact, rather than reducing the Ca ions, the excess electron prefers to be solvated by the water molecules [25]. Such solvated electrons are well-known precursors of H₂. A similar situation was observed for tobermorite-11 in Ref. [23]. While the Ca/Si (C/S) ratio of this material (C/S = 0.7 – 0.8) is significantly lower than experimental values determined for typical C–S–H gels (C/S = 1.8) [26] and than the one obtained for the present model (C/S = 1.64) [5], it appears that the results for electron localization are remarkably independent on the C/S ratio.

The excess hole appears to localize mainly on a single oxygen which forms either part of a SiO₄ tetrahedron, a water molecule or an OH⁻ group in the interlayer spacing. This is to be expected, as electrons are removed from the highest occupied molecular orbitals which, over a wide range of oxides, are located on the oxygen atoms. The peculiar characteristic here is

that there are several different nonequivalent oxygen species in cement. When an electron is removed from the hydroxyl anion, the latter becomes the highly reactive hydroxyl radical OH•. Alternatively, the electron can be removed from a water molecule producing the H₂O•⁺ radical cation, which then dissociates into OH• and H⁺. When reduced by a solvated electron, H⁺ forms a hydrogen atom, which is a precursor of H₂. For the hole localization, there are significant differences between models. For example, considering two variants of tobermorite-11, it has been observed that holes localize on the intralayer or interlayer SiO₄ oxygens depending on the C/S ratio [23]. Holes localized in the intralayer space are unlikely to be involved in further reactions. Studying chemical reactions that follow the generation of excess electrons and holes will be subject of future investigations.

## 4 Discussion and conclusions

Sr is known to replace Ca due to their isovalent properties and, indeed, ⁹⁰Sr has been detected in the concrete walls of buildings nearby nuclear accident locations like Fukushima and Chernobyl [27]. Sr is also known to replace Ca in bone [3], potentially leading to diseases like leukemia due to its proximity to bone marrow, where stromal cells are originated and stored [28]. It is therefore of great interest to understand what are the consequences of β-decay.

In this paper, we have considered the chain of β-decay events ⁹⁰Sr → ⁹⁰Y → ⁹⁰Zr in a recently proposed model for cement with density and C/S ratio consistent with experimental data [5]. We focused on the three aspects of the decay event: nuclear transmutation, recoil, and β-electron track.

We showed that, contrary to the gas-phase phenomenon of shake-off by which, after transmutation, an important channel for dissipating the excess energy is the ejection of a second electron, in the condensed phase there is a quick contraction of the electronic density followed by the capture of an electron from the surrounding medium for neutralization purposes, which is completed in less than a femtosecond, well before any atom can move. Therefore, this aspect of the decay does not lead to any significant direct effects on the cement matrix. This observation stands for both decay events, ⁹⁰Sr and ⁹⁰Y.

Due to momentum conservation and the very small mass ratio between the ejected electron and the recoil nucleus, the kinetic energy of the daughter nucleus in the recoil event is rather small. While it has been shown that it can lead to chemical transformations in molecular systems [19], in the cement matrix the effects are quite small. The recoil energy for the ⁹⁰Y atom is smaller than 5 eV, which is below the displacement threshold energy, and it is rapidly dissipated into the environment. Moreover, since Y and

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Spin density (yellow surface) representing the excess electron (left) and hole (right) in cement, generated by irradiation</figcaption></figure>
Sr are both stable in the Ca site [6], the first decay
$^{90}\text{Sr} \rightarrow ^{90}\text{Y}$ does not lead to any significant local effects.
The second decay $^{90}\text{Y} \rightarrow ^{90}\text{Zr}$ is however a different
matter. There the larger energy of the $\beta$-electron leads
to a recoil energy of up to 20 eV, which is sufficient
to provoke the translocation of the Zr atom. This,
together with the observation that the Zr atom is
not stable in the Ca location, preferring the Si site
in the C-S-H matrix [6], suggest that the second $\beta$-
decay has a small, but observable effect on the local
structure.

Finally, we have studied the consequences of the
irradiation of the cement by the β-electron. Instead
of simulating the radiation cascade, we assumed that
the main effect was to ionize the material via elec-
tron impact, thus leaving holes and excess electrons
in its tracks. By means of first-principles calculations
we showed that secondary electrons tend to localize
in the interlayer region which is populated with Ca
ions, water molecules and OH⁻ ions. First-principles
MD simulations show that the electron is very mobile
in this region, eventually becoming solvated by water
and reacting with protons to form H• and hence H₂. A
similar situation occurs with the holes localized in the
interlayer region but holes are stable in the intralayer
region.

An important aspect that was not addressed here,
and which deserves investigation, is the radiation cas-
cade due to the β-electron, i.e., which species are cre-
ated and where. In addition to ionization, there is also
vibrational and electronic excitation, which can also
lead to chemical reactions and the formation of other
species. Once it is known what species are present, then
it would be interesting to address the fate of these
species, their solubility in water, transport properties,
segregation, and the possibility that they compromise
the structural stability of the cement producing cracks
that allow radioactive material to leach into the soil
and ground waters. Within this context, it is also of
great importance to study the fate of ¹³⁷Cs in cement,
which has a great propensity to leach [29]. Also of major
interest is to study the effect of the decay of ⁹⁰Sr and
¹³⁷Cs in stainless steel, which is the other structural
material used as a barrier to contain radioactive waste
[30].

**Acknowledgements** A.A.C. work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. A.A.C. was supported by the U.S. Department of Energy, Office of Science, Materials Sciences and Engineering Division. Computing support for this work came from the Lawrence Livermore National Laboratory Institutional Computing Grand Challenge program.

**Data Availability Statement** This manuscript has no associated data or the data will not be deposited. [Authors' comment: All the results presented in the paper can be reproduced with the information provided in the manuscript. The authors can be contacted for further information.]

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References

1. A.B. Veksler, V.F. Fisenko, Hydrotech. Constr. **31**, 106 (1997). https://doi.org/10.1007/bf02766868

2. L.B. Silverio, Lamas W. de Queiroz, Energy Policy **39**, 281 (2011). https://doi.org/10.1016/j.enpol.2010.09.040

3. E. Storey, Arch. Biochem. Biophys. **124**, 575 (1968).
https://doi.org/10.1016/0003-9861(68)90368-8

4. T.A. Carlson, C.W. Nestor, T.C. Tucker, F.B. Malik, Phys. Rev. **169**, 27 (1968). https://doi.org/10.1103/physrev.169.27

5. M.J. Abdolhosseini Qomi, K.J. Krakowiak, M. Bauchy, K.L. Stewart, R. Shahsavari, D. Jagannathan, D.B. Brommer, a Baronnet, M.J. Buehler, S. Yip, F.-J. Ulm, K.J. Van Vliet, R.J.-M. Pellenq, Nat. Commun. **5**, 4960 (2014). https://doi.org/10.1038/ncomms5960

6. L. Dezerald, J.J. Kohanoff, A.A. Correa, A. Caro, R.J.-M. Pellenq, F.J. Ulm, A. Saúl, Environ. Sci. Technol. **49**, 13676 (2015). https://doi.org/10.1021/acs.est.5b02609

7. K. Lejaeghere, G. Bihlmayer, T. Bjorkman, P. Blaha, S. Blugel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. Dal Corso, S. de Gironcoli, T. Deutsch, J. K. Dewhurst, I. Di Marco, C. Draxl, M. Du ak, O. Eriksson, J. A. Flores-Livas, K. F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Granas, E. K. U. Gross, A. Gulans, F. Gygi, D. R. Hamann, P. J. Hasnip, N. A. W. Holzwarth, D. Iu an, D. B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Kucukbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordstrom, T. Ozaki, L. Paulatto, C. J. Pickard, W. Poelmans, M. I. J. Probert, K. Refson, M. Richter, G.-M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunstrom, A. Tkatchenko, M. Torrent, D. Vanderbilt, M. J. van Setten, V. Van Speybroeck, J. M. Wills, J. R. Yates, G.-X. Zhang, S. Cottenier, Science **351**, aad3000 (2016) https://doi.org/10.1126/science.aad3000

8. A. Messiah, *Title Quantum Mechanics*, series Dover books on physics No. number v. 2 (Dover Publications, 1999)

9. E.W. Draeger, X. Andrade, J.A. Gunnels, A. Bhatele, A. Schleife, A.A. Correa, J. Parallel Distrib. Comput. **106**, 205 (2017). https://doi.org/10.1016/j.jpdc.2017.02.005

10. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococ-
cioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A.P. Seitsonen, A. Smogunov, P. Umari, R.M. Wentzcovitch, J. Phys. Condens. Matter **21**, 395502 (2009). https://doi.org/10.1088/0953-8984/21/39/395502

11. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. **77**, 3865 (1996). https://doi.org/10.1103/PhysRevLett.77.3865

12. T.D. Kühne, M. Iannuzzi, M. Del Ben, V.V. Rybkin, P. Seewald, F. Stein, T. Laino, R.Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M.H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A.S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G.K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C.J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, J. Hutter, J. Chem. Phys. **152**, 191403 (2020). https://doi.org/10.1063/5.0007045

13. S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B **54**, 1703 (1996). https://doi.org/10.1103/PhysRevB.54.1703

14. J. VandeVondele, J. Hutter, J. Chem. Phys. **127**, 114105 (2007). https://doi.org/10.1063/1.2770708

15. Y. Zhao, D. Truhlar, Theor. Chem. Accounts **120**, 215 (2008). https://doi.org/10.1007/s00214-007-0310-x

16. R. Serber, H.S. Snyder, Phys. Rev. **87**, 152 (1952). https://doi.org/10.1103/PhysRev.87.152

17. J.S. Levinger, Phys. Rev. **90**, 11 (1953). https://doi.org/10.1103/PhysRev.90.11

18. W. Tang, E. Sanville, G. Henkelman, J. Phys. Condens. Matter **21**, 084204 (2009). https://doi.org/10.1088/0953-8984/21/8/084204

19. N. Marks, D.J. Carter, M. Sassi, A.L. Rohl, K.E. Sickafus, B.P. Uberuaga, C.R. Stanek, J. Phys. Condens. Matter **25**, 065504 (2013). https://doi.org/10.1088/0953-8984/25/6/065504

20. F. Blanco, A. Muñoz, D. Almeida, F. Ferreira da Silva, P. Limão-Vieira, M.C. Fuss, A.G. Sanz, G. García, Eur. Phys. J. D **67**, 199 (2013). https://doi.org/10.1140/epjd/e2013-40276-1

21. V. Lipp, N. Medvedev, B. Ziaja, In *Damage to VUV, EUV, and X-ray Optics VI*, vol. 10236, edited by L. Juha, S. Bajt, and R. Soufli, International Society for Optics and Photonics (SPIE, 2017) pp. 40–48. https://doi.org/10.1117/12.2267939

22. E. García-Toraño, V. Peyres, F. Salvat, Comput. Phys. Commun. **245**, 106849 (2019). https://doi.org/10.1016/j.cpc.2019.08.002

23. S. Le Caër, L. Dezerald, K. Boukari, M. Lainé, S. Taupin, R.M. Kavanagh, C.S.N. Johnston, E. Foy, T. Charpentier, K.J. Krakowiak, R.J.-M. Pellenq, F.J. Ulm, G.A. Tribello, J. Kohanoff, A. Saúl, Cement Concr. Res. **100**, 110 (2017). https://doi.org/10.1016/j.cemconres.2017.05.022

24. S. Le Caër, Water **3**, 235 (2011). https://doi.org/10.3390/w3010235

25. C. Fourdrin, H. Aarrachi, C. Latrille, S. Esnouf, F. Bergaya, S. Le Caër, Environ. Sci. Technol. **47**, 9530 (2013). https://doi.org/10.1021/es401490t

26. A.J. Allen, J.J. Thomas, H.M. Jennings, Nat. Mater. **6**, 311 (2007). https://doi.org/10.1038/nmat1871

27. E.B. Farfán, S.P. Gaschak, A.M. Maksymenko, E.H. Donnelly, M.D. Bondarkov, G.T. Jannik, J.C. Marra, Health Phys. **101**, 311 (2011). https://doi.org/10.1097/HP.0b013e3182103242

28. S. Musilli, N. Nicolas, Z. El Ali, P. Orellana-Moreno, C. Grand, K. Tack, S. Kerdine-Römer, J.M. Bertho, Sci. Rep. **7**, 41580 (2017). https://doi.org/10.1038/srep41580

29. K.G. Papadokostaki, A. Savidou, Hazard Mater. **171**, 1024 (2009). https://doi.org/10.1016/j.jhazmat.2009.06.118

30. A.R. Lang, D.L. Engelberg, C. Walther, M. Weiss, H. Bosco, A. Jenkins, F.R. Livens, G.T.W. Law, ACS Omega **4**, 14420 (2019). https://doi.org/10.1021/acsomega.9b01311