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. \(\mathrm {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 \(\mathrm {{}^{90}Sr}\) (28.8 years) and \(\mathrm {{}^{137}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 \(\mathrm {Sr}\) that is isovalent with \(\mathrm {Ca}\) and has been shown to replace it in bones and teeth [3]. This ability to replace \(\mathrm {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 \(\mathrm {Sr}\) and \(\mathrm {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 \(\mathrm {{}^{90}Sr}\) substitution in cement. This species undergoes \(\beta \)-decay to \(\mathrm {{}^{90}Y}\) with a half-life of 28.8 years which, in turn, \(\beta \)-decays to stable \(\mathrm {{}^{90}Zr}\) with a half-life of 64 hours. There are several aspects of \(\beta \)-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 a bulk (solid) system. Secondly, the emitted \(\beta \)-electron typically travels a few \(\mathrm {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 \(\beta \)-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 \(\mathrm {Ca_{72}Si_{44}O_{235}H_{150}}\), with a C/S ratio of 1.64 and a density of \(2.33~\mathrm {g/cm^3}\), in very good agreement with experimental data. Figure 1 shows the silicate chains formed by corner sharing \(\mathrm {SiO_4^{-4}}\) tetrahedra. It can also be seen that the \(\mathrm {Ca^{2+}}\) ions appear in two distinct environments: intralayer and interlayer. Water, in the form of \(\mathrm {H_2O}\) molecules and \(\mathrm {OH^-}\) groups, is present in the interlayer space.

Fig. 1
figure 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 \(\mathrm {SiO_{4}^{4\mathrm{-}}}\) tetrahedra forming the silicate chains are also shown. The interlayer \(\mathrm {Ca}\) atom (marked as ‘Ca56’) is substituted by \(\mathrm {Sr}\), which later decays into \(\mathrm {Y}\) and subsequently into \(\mathrm {Zr}\). The marked atom (Ca56) is in the initial position of the Zr and Y atoms for the recoil simulations

Full size image

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.

Starting from the structure representative of cement shown in Fig. 1, we replace a \(\mathrm {Ca}\) with a \(\mathrm {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 \(\mathrm {Sr}\) replaced the original \(\mathrm {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 \(\mathrm {Sr}\) nucleus with an \(\mathrm {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 \(\mathrm {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 \(\mathrm {Ca_{71} Y Si_{44} O_{235} H_{150}}\) 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 calculation 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~\mathrm {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~\mathrm {Ry}\), and all the calculations were performed at constant volume with a single \(k\)-point (the \(\varGamma \)-point).

2.2.3 Excess holes and electrons

The passage of the \(\beta \)-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 and excess hole in the sample shown in Fig. 1, i.e., without \(\mathrm {Sr}\) substitution. This is a valid procedure because the concentration of \(\mathrm {Sr}\) is relative small (maximum 4%, given by the solubility limit [6]) and the chances that charge defects localize next to another \(\mathrm {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 \(\varGamma \)-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 \(\beta \)-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 \(\beta \)-electron emission.

3.1 Sudden increase of the nuclear charge

3.1.1 Probability of \(\beta \)-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 \(\beta \)-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 \(\mathrm {Sr^{2+}}\) and \(\mathrm {Y^{3+}}\) are shown in Fig. 2 together with the changes upon transmutation. Also reported are the overlap integrals of the radial wavefunctions,

$$\begin{aligned} w_{n\ell } = \int _0^\infty {{\tilde{P}}}_{n\ell }(r) P_{n\ell }(r) \mathrm {d}r, \end{aligned}$$
(1)

where \({{\tilde{P}}}_{n\ell }(r)\) is the radial wavefunction of \(\mathrm {Sr^{2+}}\), and \(P_{nl}(r)\) is that of \(\mathrm {Y^{3+}}\).

Table 1 Shake-off ionization probabilities (in percent) for the different electronic shells of \(\mathrm {Y^{3+}}\) after transmutation from \(\mathrm {Sr^{2+}}\)
Full size table
Fig. 2
figure 2

The radial wavefunctions \({{\tilde{P}}}_{n\ell }(r)\) of \(\mathrm {Sr^{2+}}\) (solid red) and \(P_{n\ell }(r)\) of \(\mathrm {Y}^{3+}\) (dashed blue) versus r in [Å]. The numerical value on each panel is the overlap integral \(w_{n\ell }\) between the radial wavefunctions. Solid green curves show \(({\tilde{P}}_{n\ell }(r) - P_{n\ell }(r)) P_{n\ell }(r)\), whose integral from \(0\) to \(\infty \) gives the difference between \(w_{n\ell }\) and \(1\)

Full size image

The probability of ionization for the \(n\ell \) subshell is given by:

$$\begin{aligned} p_{nl} = 1 - \left( w_{n\ell }^2\right) ^{N_{n\ell }} \simeq N_{n\ell } (1 - w_{n\ell }^2) \end{aligned}$$
(2)

where \(N_{n\ell } = 2(2\ell +1)\) is the number of electrons in the \(n\ell \) subshell. Note that Eq. (2) is valid only when \(p_{n\ell } \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.

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 \(\mathrm {Ca}\) with a \(\mathrm {Sr}\) atom and relaxed the structure using ab initio DFT-GGA electronic structure calculations. We next replaced the \(\mathrm {Sr}\) with a \(\mathrm {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 \rightarrow 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 \(\mathrm {Ca_{71} Y Si_{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 \(\mathrm {Y}\) atom. After \(0.5~\mathrm {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 \(\mathrm {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).

Fig. 3
figure 3

Valence charge around the \(\mathrm {Sr/Y}\) site during the first femtosecond after the \(\beta \)-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 \(\mathrm {Sr}\) and final \(\mathrm {Y}\) to the corresponding ground state values (for fixed ion positions)

Full size image

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~\mathrm {fs}\) after transmutation.

Fig. 4
figure 4

Change in electronic charge distribution \(0.3~\mathrm {fs}\) after the \(\beta \)-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

Full size image

The presence of neighboring atoms chemically bound to the \(\mathrm {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 replacement of \(\mathrm {Sr}\) by \(\mathrm {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 around the \(\mathrm {Y}\)-atom in a longer time scale. They are in addition to the sudden forces induced by the recoil of the \(\mathrm {Y}\) nucleus at the moment of the \(\beta \) ejection.

Fig. 5
figure 5

Upper panel: structure of first neighbors around \(\mathrm {Y}\) (before \(\beta \)-decay). Lower panel: forces on the neighboring \(\mathrm {O}\)-atoms in the direction of the Y–O bonds, as a function of time after \(\beta \)-decay

Full size image

We will see in the next section that the \(\mathrm {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 \(\beta \)-decay from \(\mathrm {Y}\) to \(\mathrm {Zr}\) occurs (the lifetime of the \(\mathrm {Y}\) nucleus is 64 hours), the system is in its new equilibrium state. This second \(\beta \)-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 describ- ed 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 \(\mathrm {SrTiO_3}\) and \(\mathrm {SrH_2}\) [19]. In our case, the maximum kinetic energy of the recoiling atoms is of the order of \(5~\mathrm {eV}\) for the first decay, \({}^{90}\mathrm {Sr}\) to \({}^{90}\mathrm {Y}\), and \(20~\mathrm {eV}\) for the second one, \({}^{90}\mathrm {Y}\) to \({}^{90}\mathrm {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 \(\mathrm {eV}\).

Fig. 6
figure 6

Molecular dynamics simulation of the recoil excerpted on the \(\mathrm {Zr}\) daughter atom. The calculations have been performed using the Quantum Espresso code [10] with a time step of \(1~\mathrm {fs}\). a Evolution of the kinetic and potential energy of the simulation box. b Relative displacement of the \(\mathrm {Zr}\) atom with respect to its time equal zero position. c Initial atomic configuration. d Final atomic configuration

Full size image

Since the outcome of the recoil process is not obvious, we have conducted adiabatic first-principles molecular dynamics simulations where first the \(\mathrm {Y}\), and then the \(\mathrm {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 \(\mathrm {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 \(\mathrm {Zr}\).

The results corresponding to the \(\mathrm {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~\mathrm {ps}\) at \(300~\mathrm {K}\). Then the \(\mathrm {Zr}\) kinetic energy has been suddenly increased to \(20~\mathrm {eV}\) by giving the \(\mathrm {Zr}\)-atom a suitable velocity along the desired direction. This can be seen as a jump in the the total kinetic energy at  \(t = 0.2~\mathrm {ps}\). The simulation then continued during \(1~\mathrm {ps}\) where after a short transient of about \(100~\mathrm {fs}\) the excess kinetic energy is spread among all the atoms in the simulation box reaching a final temperature of about \(450~\mathrm {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 \(\mathrm {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 \(\beta \)-electron in its tracks

3.3.1 Radiation cascade

The \(\beta \)-electron resulting from the decay of \(^{90}\)Sr to \(^{90}\)Y is emitted with an energy of \(0.5~\mathrm {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~\mathrm {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,21,22]. This approach uses the energy distribution of the \(\beta \)-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 \(\mathrm {H}_2\) from water dissociation [23]. The accumulation of \(\mathrm {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 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.

Fig. 7
figure 7

Spin density (yellow surface) representing the excess electron (left) and hole (right) in cement, generated by irradiation

Full size image

The excess electron in C–S–H is observed to localize in the calcium and water-rich interlayer regions across the water molecules, \(\mathrm {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 \(\text {H}^{\bullet }\) and \(\mathrm {H_2}\) [24]. In fact, rather than reducing the \(\mathrm {Ca}\) ions, the excess electron prefers to be solvated by the water molecules [25]. Such solvated electrons are well-known precursors of \({\mathrm {H}}_2\). A similar situation was observed for tobermorite-11 in Ref. [23]. While the \(\mathrm {Ca/Si}\) (\(\mathrm {C/S}\)) ratio of this material (\(\mathrm {C/S} = 0.7-0.8\)) is significantly lower than experimental values determined for typical C–S–H gels (\(\mathrm {C/S}=1.8\)) [26] and than the one obtained for the present model (\(\mathrm {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 \(\mathrm {SiO_4}\) tetrahedron, a water molecule or an \(\mathrm {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 \(\text {OH}^{\bullet }\). Alternatively, the electron can be removed from a water molecule producing the \(\text {H}_2\text {O}^{\bullet +}\) radical cation, which then dissociates into \(\text {OH}^{\bullet }\) and \(\mathrm {H^+}\). When reduced by a solvated electron, \(\mathrm {H^+}\) forms a hydrogen atom, which is a precursor of \(\mathrm {H_2}\). 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 \(\mathrm {SiO_4}\) 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

\(\mathrm {Sr}\) is known to replace \(\mathrm {Ca}\) due to their isovalent properties and, indeed, \(\mathrm {{}^{90}Sr}\) has been detected in the concrete walls of buildings nearby nuclear accident locations like Fukushima and Chernobyl [27]. \(\mathrm {Sr}\) is also known to replace \(\mathrm {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 \(\beta \)-decay.

In this paper, we have considered the chain of \(\beta \)-decay events \(\mathrm {{}^{90}Sr \rightarrow {}^{90}Y \rightarrow {}^{90}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 \(\beta \)-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 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, \(\mathrm {{}^{90}Sr}\) and \(\mathrm {{}^{90}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 \(\mathrm {{}^{90}Y}\) atom is smaller than \(5~\mathrm {eV}\), which is below the displacement threshold energy, and it is rapidly dissipated into the environment. Moreover, since \(\mathrm {Y}\) and \(\mathrm {Sr}\) are both stable in the \(\mathrm {Ca}\) site [6], the first decay \(\mathrm {{}^{90}Sr \rightarrow {}^{90}Y}\) does not lead to any significant local effects. The second decay \(\mathrm {{}^{90}Y \rightarrow {}^{90}Zr}\) is however a different matter. There the larger energy of the \(\beta \)-electron leads to a recoil energy of up to \(20~\mathrm {eV}\), which is sufficient to provoke the translocation of the \(\mathrm {Zr}\) atom. This, together with the observation that the \(\mathrm {Zr}\) atom is not stable in the \(\mathrm {Ca}\) location, preferring the \(\mathrm {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 \(\beta \)-electron. Instead of simulating the radiation cascade, we assumed that the main effect was to ionize the material via electron 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 \(\mathrm {Ca}\) ions, water molecules and \(\mathrm {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 \(\text {H}^{\bullet }\) and hence \(\mathrm {H_2}\). 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 cascade due to the \(\beta \)-electron, i.e., which species are created 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 \(\mathrm {{}^{137}Cs}\) in cement, which has a great propensity to leach [29]. Also of major interest is to study the effect of the decay of \(\mathrm {{}^{90}Sr}\) and \(\mathrm {{}^{137}Cs}\) in stainless steel, which is the other structural material used as a barrier to contain radioactive waste [30].