Clémentine Panetier

Antonia Ruiz Moreno

François Rossi

Théo Roubille

Gašper Žerovnik

Arjan Plompen

Nathalie Moncoffre

Yves Pipon
email: yves.pipon@univ-lyon1.fr

A Ruiz-Moreno

G Zerovnik

Molecular dynamics simulations of Mo nanoparticles sputtering under irradiation

Keywords: Molecular Dynamic, Nanoparticle, Kinematic recoil method, Irradiation, Sputtering

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The use of radioisotopes has become an essential part of the innovative medical treatments either for diagnostics, therapy or the combination of both (theranostics). 99m Tc- 99 Mo generators are needed to respond to the worldwide demand of 99m Tc radionuclide used in nuclear medicine as tracer for diagnostic imaging. Up to now, about 90% of 99 Mo, the parent radionuclide of 99m Tc, is produced from fission occurring in highly enriched uranium (HEU) targets. In the world, only five nuclear reactors work with HEU. They were built in the 60's and are therefore reaching the end of their designed lifetime. The risk of penury induced by this situation, added to the proliferation issue due to the HEU use, show the importance of diversifying 99 Mo production. In order to achieve this nonproliferation goal, all major 99 Mo-producing countries have agreed to use low-enriched uranium (LEU) targets [1]. This decision implies to focus researches on the use of small particles accelerators. That is the reason why several accelerator-based reaction channels and methods were investigated over the years.

In their paper, Starovoitova et al. [2] discussed the possibility of photonuclear production of Mo radioisotope using linear electron accelerators. They have shown both by Monte Carlo calculations and experimental results that the specific activity was not enough for 99 Mo to be separated from 100 Mo using standard alumina columns. However, they proposed a kinematic recoil method for the separation of the isotopes that could make the photoneutron production in accelerators feasible. It involves a target surrounded by a so-called "catcher" material. Subsequently after the target irradiation, recoil nuclei with sufficient kinetic energy escape the target and are trapped in the catcher material, from which they will be further separated. To be efficient, this method requires the use of small targets such as thin foils or suspension of nanoparticles (NPs). In the case of NPs, the drawback of this method is the generation of unwanted nuclei mostly from the NP surface sputtering. Dikiy et al. studied MoO3 NPs suspended in isopropyl alcohol and irradiated by Bremsstrahlung with Emax = 12.5 MeV [3]. Being well aware of the sputtering problem, each NP was covered with a 0.5 nm layer of ethylene glycol. Because of this additional layer, NP size should be the smallest possible to allow the escape of the 99 Mo nuclei which complicates the NP elaboration part. Another possible nuclear reaction is the 100 Mo(n, 2n) 99 Mo reaction. Nagai et al. used a 40 MeV deuteron ion beam to irradiate carbon foils which produce neutrons of around 14 MeV [4]. These neutrons were guided towards a MoO3 sample and a sublimation separation of 99m Tc was operated from the 100 MoO3 sample by using a dedicated furnace.

It should be mentioned that 99 Mo activity produced by fission reaction is around 10 3 -10 4 times higher than any alternative production method [5]. It is therefore really important to optimize each method in order to maximize each production yield. Several private-sector companies in the United States are planning to produce 99 Mo for medical use such as SHINE Medical Technologies which plan to use deuterium/tritium accelerator technology to induce subcritical fissioning of 235 U in an LEU uranyl sulfate solution. We could also mention the TRIUMF's project (Canada) [6], which aims to produce 99m Tc via the 100 Mo(p,2n) 99m Tc reaction using cyclotrons.

In the present paper, we propose to explore more in-depth the use of "small" ion beam accelerators to produce 99 Mo with the objective to reduce at maximum the complexity of the process (for instance avoiding any NP coating). The initial setup would consist in the immersion of 100 Mo NPs in liquid followed by irradiation in the tandem accelerator MONNET or the linear accelerator GELINA (located in Geel, Belgium) [7] to enforce the two nuclear reactions 100 Mo(n,2n) 99 Mo or 100 Mo(γ,n) 99 Mo respectively. After irradiation, the 99 Mo isotopes could be extracted by solvent extraction or mass separation techniques such as centrifugation. However, performing an experimental parametric study, e.g. by varying NP size and concentration, to maximize the production yield is rather resource intensive. In this regard, atomistic simulations can be very useful to theoretically evaluate the outcome of the reactions and guide the experimental efforts. Most of the time, Monte Carlo calculations with codes such as GEANT4 [8] or SRIM [START_REF] Ziegler | Stopping and Ranges of Ions in Matter[END_REF] are performed in order to evaluate the sputtering yield and the projected range of the recoil atoms. However, the geometry of the target can impact the sputtering and these codes are rather limited for that (especially SRIM). Moreover, the irradiation physics behind these Monte Carlo codes are based on binary collision approximation (BCA), which is not accurate in order to simulate recoil atoms dynamics with low kinetic energy. Indeed multiple collision sequences and thermal spikes are the dominant effects because of the relative NP and collision cascade sizes. In this paper, we propose to use Molecular Dynamics (MD) simulations with semi-empirical potentials to study the behaviour of Mo NPs of two different sizes: 5 nm, chosen because a stable solution of smaller particles would be very difficult to achieve experimentally, and 15 nm which would be a good compromise to compare bigger NPs and save calculation time. The NPs were placed either in vacuum or in water for three different irradiation conditions. Indeed, we used three different energies for the primary knock-on atom (PKA): 19 keV for the photonuclear reaction on one hand and 100 and 150 keV for simulating the effects of incoming neutrons between 12 and 16 MeV on the other hand. These energy values have been determined mostly by Monte Carlo calculations as it is explained in the following paragraph.

Method and simulation parameters

Calculation of the recoil energies to be used in MD simulations

The primary recoil energy spectra resulting from 100 Mo(γ,n) 99 Mo reaction were determined for a generic irradiation setup with a monoenergetic 75 MeV electron beam. This requires combining three calculations. The first one determines the photon energy spectrum in the sample, the second one the neutron spectra of the (,n) reaction and the third one uses the neutron spectra to calculate the recoil spectra. The setup consisted of an electron beam of 1 cm diameter surrounded by lead with a 5 mm thick tantalum radiator at the end to produce high energy photons by Bremsstrahlung. As sample, a 5 mm thick, 1 cm diameter disc of water was placed at a distance of 5 mm from the tantalum radiator. The photon flux spectrum in the water sample was calculated using MCNP code [START_REF] Pelowitz | MCNPX Users Manual Version 2[END_REF] and assuming a normal incidence of the electron beam. Then, the neutron spectra in the center of momentum system were calculated with the nuclear reaction code TALYS [START_REF] Koning | Talys: a new code for calculating nuclear data[END_REF]. Such neutron spectra determine the recoil energy in the center of mass system due to the fact that the momenta of the neutron and the recoil in that system are equal in size but directed oppositely. It was observed that the spectrum was dominated by neutrons generated from a photon energy of 14 MeV, which exhibited a neutron energy peak at 0.4 MeV. Finally, the 99 Mo recoil energy was calculated using the simple kinematics of binary reactions where there is conservation of momentum, assuming that the outgoing recoil angular distribution is isotropic and neglecting the momentum of the outgoing photons. For each photon energy, the differential cross section from the TALYS center of mass neutron spectra was multiplied by the photon flux determined with MCNP. The recoil energy spectrum presented a peak at 4.9 keV and a mean energy of 19 keV in agreement with literature [START_REF] Np | [END_REF]. The mean energy of the recoil spectrum was used as PKA recoil energy for MD simulation.

Regarding the 100 Mo(n,2n) 99 Mo reaction, the recoil energy was estimated applying the momentum conservation principle and assuming incoming neutrons of 12 -16 MeV, which give a cross section of about 1.4 barns, and two outgoing neutrons of about 1 MeV. This results in 99 Mo recoil energies above 100 keV. For MD simulations, PKA energies of 100 and 150 keV were selected.

Description of the potentials used in the MD simulations

All MD simulations were performed with the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [13] code and visualization of the systems were done with the OVITO 3.0 visualization software [14].

To describe Mo pair interactions, we used the Ackland and Thetford (A&T) potential [15], which is derived from the Finnis-Sinclair potential [16]. The many-body A&T potential was shown to be well suited to simulate the production of defects (i.e. Frenkel pairs) after Mo irradiation with low energy PKA [START_REF] Np | [END_REF][13]. However, for PKA energies higher than 5 keV, it was shown by Selby et al. [17] that another term should be added to better describe the interactions at short distances (< 2 Å). This is why a purely repulsive Ziegler-Biersack-Littmark (ZBL) potential term [START_REF] Ziegler | Stopping and Ranges of Ions in Matter[END_REF] was added to account for the screened nuclear repulsion, important for better describing "high"-energy collisions between atoms. In addition, for PKA energies higher than 100 keV, a specific LAMMPS command (fix electron/stopping) was used to take into account inelastic energy loss via electronic interactions. It uses electronic stopping power values of Mo calculated with SRIM for PKA energies ranging from 10 eV to 150 keV.

The description of a water molecule can be done in various ways as many potentials exist that can take into account polarizability effects, flexibility of the water molecule geometry and other physico-chemical properties. In our study, we chose the TIP4P/2005 water model [18] which was shown to accurately simulate the water phase diagram, which is mandatory due to the expected temperature increase provoked by the recoil atoms losing their energy in water. This model exhibits four possible interaction sites: one oxygen atom, two hydrogen atoms and a pseudo atom "M" just below the O atom holding a negative charge. This allows a better distribution of the global charge of the water molecule. The O-H bond length and the H-O-H angle are fixed during the whole simulation thanks to the SHAKE algorithm [19] with a precision better than 10 -4 .

To simulate the interactions between Mo atoms and O atoms from water molecules (no Mo-H interactions were considered), we used the Lennard-Jones (LJ) potential which is described below:

V(Mo -O) = 4 [(  r ) 12 -(  r ) 6 ]
This method was previously applied with success by Köhler et al. in order to simulate MoS2 nanomembranes in water [20]. We chose to work with two mixing rules in order to calculate the LJ potentials parameters: (i) the Lorenz-Berthelot (L-B) mixing rule [21] and, (ii) the Kong rule [22] which was used once to see the impact of different mixing rules. The parameter values are given in Table 1. A cut off value was set to 10 Å as done by Heyhat et al. [23] for a {Ag -water} system. The interactions between Mo and O are therefore not calculated beyond this distance.

NP modelling in vacuum or in water

In order to obtain a stable polymorph of Mo, we used the WulffPack Python package based on the Wulff theory [24] which was developed to find the most stable nanoparticle shape. The library requires the energy of the main surfaces obtained by cutting a BCC Mo crystal with different slabs. The surface energy was therefore calculated as follows:

 = 1 2𝐴 (𝐸 𝑠𝑙𝑎𝑏 -𝑛. 𝐸 𝑏𝑢𝑙𝑘 )
With: A, being the plane surface ; Eslab, being the total surface energy ; Ebulk, being the bulk energy and n, being the number of atoms in the considered plane.

Table 2 sums up the surface energy for the main surfaces calculated in the present work. The nanoparticle shape resulting from the code and the A&T surface energies consists in a polyhedral form with only two different planes. The main one is the (110) plane with 89% of the total surface and the (100) plane with 11% of the total surface. NPs containing 4,461 and 110,485 atoms were afterwards constructed. Approximating the NP shape by a perfect sphere, it would correspond to diameters of 5 and 15 nm respectively.

First, as displayed in Figure 1, the constructed NPs were placed at the centre of an empty cube with a 50 nm length for a 15 nm NP or a 20 nm length for a 5 nm NP. The empty cube simulates a vacuum environment. Second, the NPs were placed in a big box filled with water molecules. In order to save some calculation time, the size of the cube was decreased. Then, water molecules were randomly added around the NP. The number of water molecules "N" was fixed in relation with the cube volume "V" and with the liquid water density "" according to the relation:

 = 𝑁. 𝑀 𝑁 𝐴 . 𝑉
With: M, the molecular mass of water (18 g/mol) and NA, the Avogadro constant.

It means that a 5 nm Mo NP was associated with 28,758 water molecules in a 10 nm cube (figure 2) and that a 15 nm Mo NP was associated with 208,768 water molecules in a 20 nm cube. Our approach was validated by checking the water density and the oxygen diffusion coefficient at 298 K and 1 bar pressure. The mean water density was found to be 0.997 g/cm 3 in excellent agreement with the experimental reference value of 0.9971 g/cm 3 taken from [25]. The oxygen diffusion coefficient "D" was calculated thanks to the Einstein's diffusion equation:

〈(𝑟(𝑡) -𝑟(𝑡 0 )) 2 〉 = 6𝐷(𝑡 -𝑡 0 )
We have recorded the mean square displacement (MSD = 〈(𝑟(𝑡) -𝑟(𝑡 0 )) 2 〉) evolution during 100 ps at 298 K. A diffusion coefficient mean value of 2.3×10 -9 m 2 /s was found in agreement with literature [26].

It was also important to consider the thermal exchanges between the Mo NP and water. For that, two Langevin thermostats were used: one for the NP and one for water. We first heated water and Mo at 300 K in the NPH ensemble. Then we raised the Mo NP temperature to 400 K during 100 ps.

Next, the Mo NP thermostat was disabled solely and the Mo temperature decrease was recorded during 100 ps. The timestep used is 1 fs. Figure 3 shows the Mo temperature evolution during the last 200 ps of the simulation. This procedure was repeated with the Kong mixing rule in order to test the impact of different Mo-O potentials on the thermal exchanges. The temperature decrease is rather low as only 25 K are lost after 100 ps for the two mixing rules. For comparison, a similar study was done by A. Rajabpour with a silver NP in water [27] and 100 K were lost after only 15 ps. Even if the thermal conductivity of silver is more than three times higher than that of molybdenum, the potentials used in the present study display a poor capacity for the heat exchanges. This means that the temperature values can be overestimated in our simulations. Moreover, the total number of water molecules in our different systems is rather small and can lead to an overheating of water. If we add the fact that the NP cooling is slow, the thermal equilibrium cannot be achieved in a reasonable simulation time. Therefore, we have chosen to force the water thermalization with three different velocities taken among all those tested for a Mo recoil entering the water medium (Figure 4). The first one is the strongest thermalization as only 100 steps (100×dt) are needed to reach the final temperature (298 K) preventing any temperature increase. The second one corresponds to 10,000×dt (~1ps to reach the final temperature). The last one is the weakest as 100,000 steps are needed to thermalize the system.

Irradiation parameters and data processing

We have first simulated the irradiation of a Mo bulk system for different PKA energies. Periodic boundary conditions (PBC) were therefore applied to simulate a "infinite" crystal. It is important to note that the simulation box should have a sufficient size to avoid that any atom displacement through a box border occurs and that an atom emerges again at another place of the box. Different box sizes were chosen according to the PKA energy. The maximum system size was obtained by replicating the primitive lattice 120 times in the three-space dimension corresponding to 3,456,000 atoms. This procedure allows a maximum PKA energy of 50 keV which is sufficient in order to compare with the calculations of Selby et al. [15]. Irradiation simulations were performed in the NVE ensemble during 100 ps. An initial velocity is assigned to a Mo atom; this atom is then referred as the PKA. To obtain a good statistic, at least 40 simulations were run for each energy with a PKA initial position and velocity direction chosen randomly. The direction was randomly chosen by determining the spherical coordinates of the velocity vector satisfying a chosen velocity norm (for instance: 194,8629182 nm/ps for 19 keV) for different azimuthal angles, and polar angles. An adaptive timestep is used during the collision events with a rescaling method that allows to reset the timestep to a lower value (10 -5 ps) if an atom moves more than 0.0003×a0 (with a0 the lattice parameter of Mo BCC system) within a single step. This procedure allows to decrease the simulation duration. To analyse the quantity and nature of defects induced by each collision cascade, the Wigner-Seitz (WS) method implemented in the OVITO 3.0 software was used. In this method, Voronoi cells are constructed in a perfect BCC Mo lattice. Before irradiation, each cell contains only one atom in its centre. After irradiation and the possible displacements, a cell can host zero, one, two or more atoms. The WS analysis consists in counting the atoms in each cell allowing to identify Frenkel pair defects and clusters of these defects. Figure 5 shows the mean number of defects (averaged over all simulations) with the associated standard deviation generated for each PKA energy. Results obtained by Selby et al. were compared to three sets of potentials: A&T potential, A&T potential with ZBL and finally the two last potentials with the Se term. Figure 5 clearly demonstrates the necessity to stiffen the potential interaction at short distances as the A&T potential induces a decrease of the defect number for a PKA energy higher than 5 keV. With the addition of the ZBL potential, a linear relationship is observed between the number of defects and the PKA energy, in agreement with the Norgett-Robinson-Torrens (NRT) model [28]. Finally, results display no contribution of the electronic stopping power for PKA energies lower than 50 keV. As already explained, no simulation was performed above this energy because of the size of the simulation box which would involve too many Mo atoms.

The other irradiations involve the two following finite systems (free boundary conditions): (i) Mo NP surrounded by vacuum and, (ii) Mo NP surrounded by water. Before any MD calculation, systems were optimized with the conjugate gradient (CG) and the FIRE algorithms [29] to a force norm of 10 -8 eV/Å, as suggested in [30]. For the {vacuum + NP} system, the Nosé/Hoover thermostat in NVT ensemble is used to work at 300 K. The system was equilibrated during 50 ps to ensure a good minimization as shown in [31]. The pressure remained close to 0 GPa during the simulations. For the {water + NP} system, as previously written, two Langevin thermostats were used and the pressure was maintained at 1 bar. The NP temperature was fixed at 300 K, while for the water the temperature was initiated at 1K before being increased up to 300 K during 100 ps with a timestep of 1 fs. This procedure allows to maintain an equilibrated system with normal pressure and temperature. It has to be noted that Mo atoms are stopped in the last 0.2 angstroms before crossing the border of the box in order to avoid loss of atoms.

For each system, irradiation simulations were then performed in the NVE ensemble during 100 ps as explained for the Mo bulk irradiation. Around 100 runs for each energy and NP size were performed in order to have sufficient statistics (random PKA initial location and direction). It should be mentioned that the WS method cannot be applied as the NP rotates and moves slightly after the recoils. Therefore, we have only counted the atoms outside the NP after the collision cascades without considering the defects inside the NP.

NP irradiation results

Several possibilities were identified depending on the PKA initial position and direction. The nanoparticle was divided into three regions: surface (the first atomic layer), edge (3 atomic layers after the so-called "surface") and centre (others atoms). Hence, results will vary according to the initial PKA position into one of these regions and its direction. We spotted several possibilities, depicted in The first possibility corresponds to a PKA at the surface directed outwards the NP. In this case, it escapes directly without disturbing the lattice/atoms around. There is no cascade and thus no sputtering. For a PKA located at the edge or in the centre of the NP, the displaced atoms have to go through the NP to escape, thus provoking collision cascades except in the case of a channelling process where the PKA goes through atomic planes without interaction (with a very low probability). For the other cases, the collision cascades induce sputtering, which intensity varies with the energy transfer between each collision cascade and, of course, the direction of the energetic displaced atoms. Generally, the sub cascades are the key parameter that controls the final number of sputtered atoms. During each irradiation simulation, three consecutive stages were shown to occur as follows:

1) First, the PKA of high kinetic energy, exits the NP in about 10 fs of time.

2) The energy left by the PKA in the NP leads to primary collision cascades. Some SKA (secondary knock-on atoms) with a few keV kinetic energy escape the NP within the first 50 fs.

3) Finally, atoms located at the surface are released from the NP (sputtering process). The Mo surface binding energy, which is roughly the sublimation heat, is around 6.5 eV. This value indicates the minimum energy to be transferred via the collision cascades for a Mo atom to leave the surface of the NP. After 10 ps, no more atom is displaced and the annealing phase (defect recombination) occurs.

In the following, we have divided the results into two parts to distinguish the simulations with NPs in vacuum from those in water. Each part shows the number of sputtered atoms for the three PKA energies (19, 100 and 150 keV) and for the two NP sizes (5 and 15 nm). The analysis will be now focused on the number of sputtered 100 Mo atoms and whether or not the PKA (simulating the 99 Mo atom) exits the NP.

NP irradiation in vacuum

Figure 8 shows the number of sputtered atoms for the different cases. The column "zero" means that only the PKA exits the NP (stage 1). The first vertical dashed line marks the limit between stage 2 and stage 3. For some runs (only for the {vacuum + NP} system), another stage was observed. It is visible on Figure 8 with the second vertical dashed line. Beyond this line is represented the case where the number of sputtered atoms is greater than 80 atoms. This fourth stage corresponds to a partial destruction of the nanoparticle. Figure 9 shows an example of this stage where 463 atoms exit the NP leading to initial particle fragmentation and subsequent NP reconstruction but without recovering its original shape. For the 5 nm NP, the higher the PKA energy, the lower probability fourth stage because the PKA exits the NP without much collisions. This is exactly the opposite for the 15 nm NP.

Table 3 gathers the results on the percentage of runs where the PKA escapes the NP and where no additional atoms are released (PKA escaping the NP solely). For the 5 nm NP, it should be noted that the PKA is almost always escaping the NP ("only" 92 % for the 19 keV PKA energy). Moreover, almost one third of the calculations simulating neutron irradiation (PKA of 100 and 150 keV initial recoil energy) show that the PKA escapes solely the NP. This corresponds to the best case, as it means that there would be no contamination of unwanted isotopes in the trapping material. It can be observed in Figure 8 that in the majority of the runs, only a few atoms (less than 10) are sputtered.

On the contrary, for the 15 nm NP, only two third of our simulations lead to the exit of the PKA and only 5% of the PKA exit solely the NP for the calculations simulating the photonuclear reaction (PKA of 19 keV initial recoil energy).

For the calculations simulating the neutron irradiation, the PKA is quasi always exiting the NP but it is almost always (~90% of the time) with a large number of sputtered atoms as shown in figure 8.

NP irradiation in water

The simulations are time-consuming and this is why the 150 keV PKA energy simulation was not done and the 15 nm NP was only simulated with a single thermalization velocity. It should be also noted that the fourth stage does not exist as the water molecules around the NP prevent any fragmentation. Indeed, the thermalization velocity significantly decreases the number of sputtered atoms. This effect is shown in Figure 10 which presents the number of sputtered atoms for a 19 keV PKA. For the strongest water thermalization (100 dt), the majority of the sputtered atoms (> 70 %) in the different runs are less than 5. The two other thermalizations (10,000 and 100,000 dt) are quite similar and less than 50 % of the simulations show 5 or less sputtered atoms. Table 4 sums up the main information for the calculations simulating the photonuclear reaction. The thermalization velocity does not affect much the exit of the PKA as expected. However, the thermalization velocity impacts the PKA exiting the NP sole or with other Mo atoms.

Of course, the strongest thermalization velocity corresponds to the most interesting case with a higher chance for the PKA to exit alone.

Figure 11 and Table 5 show the results for the 100 keV PKA energy. Because of its high energy, the PKA almost always exits the NP. The thermalization velocity affects similarly the number of sputtered atoms as the previous case. We can also note that the number of sputtered atoms is lower than for the 19 keV PKA energy.

Discussion

Comparison of the different simulation parameters

The first point to consider is the environment surrounding the Mo NP, water or vacuum. The main difference is the ability of the Mo NP to break apart (4 th stage) and, of course, an increased number of sputtered atoms in vacuum. This is explained by the fact that atoms ejected from the NP with very low kinetic energy return, most of the time, to their initial site in presence of water. This effect is enhanced as the thermalization velocity increases. With the 100dt thermalization velocity, the probability to obtain less than 10 sputtered atoms is higher than 85 % (whatever the PKA energy or NP size) which is a promising result, as this indicates that the product would have a good purity and specific activity. Outside this effect, the same trends have been observed whatever the Mo environment. It is clearly shown that the probability is higher when the NP is surrounded by water.

Let us now compare the results depending on the PKA initial energy and on the NP size. A 19 keV initial energy PKA exits most of the time (~90% of the runs) from a 5 nm NP but only at ~65% from a 15 nm Mo NP. This was expected as the projected range of a 99 Mo in a Mo matrix is of 7 nm (from SRIM calculations). Clearly, in order to maximize the probability of 99 Mo to exit the NP, the NP size should be as small as possible. This is, of course, different for a PKA energy of 100 keV as the PKA quasi always exits the NP whatever its size (at least 96 % of the runs). Again, this was expected as the projected range of 99 Mo in a Mo matrix is of 24 nm (from SRIM calculations). This indicates that this nuclear reaction would allow the use of larger NP targets than for the photonuclear reaction. sputtered atoms (over all runs) divided by the number of times the PKA exits the NP (over all runs). Without surprise, the ratio is very important for the Mo NP in vacuum which, probably, does not reflect the reality. The water presence significantly lowers the number of sputtered atoms and therefore the 100 Mo/ 99 Mo ratio. As observed, the ratio for the 100 dt thermalization velocity is ~3 times lower than the ratios of the other velocities which are roughly similar. Clearly, based on Table 6 values, the 5 nm NP irradiated with neutrons (100 keV PKA) corresponds to the best condition to avoid unwanted nuclei in water. For NP with increasing size, the ratio is lower for the photonuclear reaction but the values are somewhat close.

Limitations

Even though the results obtained are encouraging, some parameters should be considered in order to be closer to the real experimental conditions.

The first parameter to consider is the irradiation fluence of the NPs. We have performed a limited number of runs (100) in order to have sufficient statistics because the PKA initial site and direction can lead to different results depending on the different collision cascades that occur during irradiation. Even if the MD approach is mandatory to be closer to the NP irradiation process, statistics are also mandatory to consider all cases and have a mean picture.

We have not considered the case of a NP impacted several times by the irradiation. However, we have observed that the NP annealed itself rather quickly (~10ps) and the probability to have another reaction in this short duration is really low.

It should also be considered that the PKA exiting the NP could potentially impact another NP provoking additional sputtering. To check this effect, we have determined the energy distribution of the sputtered atoms in water. Then, the energy corresponding to the 95 th percentile of this distribution was calculated, which allows us to estimate with SRIM, the minimum projected range of most atoms in water. Results are gathered in table 7 for a 100 keV PKA (the 19 keV PKA case being less probable). Even with a low thermalization energy, these results show that most atoms keep at least 3 keV in water leading to a projected range higher than 10 nm. This result can be very problematic for the neutron irradiation as it could increase the number of unwanted 100 Mo atoms in water especially if there are some agglomerates. It was decided to simulate a small aggregate of 5 nm Mo NPs as shown in Figure 12. In this example, the PKA direction was chosen in order to have an impact on the NP close to the one hosting the PKA. It is clear that the Mo NP located above the one being irradiated induced a few more sputtered atoms in water. As it can been seen, if the agglomerate contains a lot of NPs, multiple NPs will be impacted and the effect will be enhanced but difficult to determine by calculations.

Conclusion

This paper shows that atomistic simulations (MD with semiempirical potentials) can bring useful data to guide experimental efforts on the development of irradiation setups. Two different irradiation types (photonuclear and neutron) were simulated to explore the possibility of a new route to produce 99 Mo isotopes from Mo nanoparticles. From our results, the neutron irradiation in 5 nm NPs seems to be the best choice to obtain a maximum of 99 Mo while limiting sputtering whereas larger NP sizes lead to a decrease of the process efficiency.

It must be emphasized that experimentally, monodisperse NPs are difficult to produce either by chemical or physical (cluster beam) synthesis. Polydisperse NPs are normally obtained leading to agglomerates. Agglomerates can increase the unwanted sputtered atoms as the PKA can escape from one single NP and interact with other NPs from the same agglomerate. This is expected to decrease the efficiency of the whole process and should be checked in further experimental works.

Other engineering challenges related to this new production would be to consider the incident neutron self-absorption which lowers the flux, the reactivity of the Mo NPs in suspension and separation of the produced isotopes from the treated NP. On the other hand, recirculating the NPs within the beam would allow an increase of the production yield which could counterbalance the drawbacks.

Altogether, these results and considerations show that the proposed process is worth further experimental validation and production of Mo medical isotopes by accelerators instead of nuclear reactor needs to be further studied experimentally, using low dimension monodisperse nanoparticles suspension.

Figure 1 -

 1 Figure 1-Stable Mo nanoparticle (15 nm diameter) defined with the WulffTheory and displayed in a vacuum cube of 50 nm length.

Figure 2 -

 2 Figure 2-Section of a 10 nm cube containing a 5 nm Mo NP and 28'758 water molecules randomly placed around the NP.

Figure 3 -

 3 Figure 3-Temperature evolution of the Mo NP. First, a Langevin thermostat allowed to rise the temperature to 400 K during 100 ps. Then, the thermostat was disabled and the Mo NP transfers his heat to the water.

Figure 4 -

 4 Figure 4-Comparison of different thermalization velocities of water.

Figure 5 -

 5 Figure 5-Number of defects (NPF) created for each PKA energy as a function of potential parameters used and electronic stopping power (Se).

Figure

  Figure Tree of possibilities after irradiation of the Mo NP according to the PKA initial position and velocity direction.

Figure 7

 7 Figure 7 illustrates the three stages occurring during the irradiation process of a 5 nm Mo NP in vacuum for a PKA energy of 100 keV.

Figure 7 -

 7 Figure 7-Three stages occurring for the irradiation of a 5 nm Mo NP in vacuum for a PKA of 100 keV.

Figure 8 -

 8 Figure 8-Number of sputtered atoms for a 5 nm NP (lower case) and for a 15 nm NP (upper case). Zero means that only the PKA escapes the NP (no sputtering). Total number of runs = 100.

Figure 9 -

 9 Figure 9-Fourth stage of a 5 nm NP irradiation with a 100 keV PKA. Small fragments are ejected (463 atoms on the initial 4,461 atoms) before the NP annealed itself.

Figure 10 -

 10 Figure 10-Number of sputtered atoms for the 19 keV PKA energy in water.

Figure 11 -

 11 Figure 11-Number of sputtered atoms for the 100 keV PKA energy in water.

Figure 12 -

 12 Figure 12-Section of different NPs close to each other and impacted by a 100 keV PKA

Table 1 -

 1 Parameters of the Mo-O potential calculated with the Lorentz

	Berthelot (L-B) or Kong mixing rules
		 (eV)	 (Å)
	L-B rules	0.002168	3.67945
	Kong rules	0.001562	2.8470

Table 2 -

 2 Surface energy of Mo bcc crystal determined with A&T potential

	Surface plane	Surface energy
		(eV/ Ų)
	(100)	0.132
	(110)	0.112
	(111)	0.148
	(210)	0.133
	(221)	0.140
	(320)	0.127
	(331)	0.133

Table 3 -

 3 Summary of the sputerring process occuring for a PKA of 19, 100 or 150 keV in a NP with a size of 5 or 15 nm (Mo NP in vacuum).

	NP size	EPKA (keV)	PKA escaping the NP (over all runs)	PKA escaping the NP solely (over all runs)
		19	92 %	17 %
	5 nm	100	100 %	33 %
		150	100 %	31 %
		19	68 %	5 %
	15 nm	100	99 %	11 %
		150	99 %	13 %

Table 4 -

 4 Summary of the sputerring process occuring for a PKA of 19 keV in a NP with a size of 5 or 15 nm for different water thermalization velocities

	NP size	Thermalization velocity	PKA escaping the NP (over all runs)	PKA escaping the NP solely (over all runs)
		100 dt	92 %	28 %
	5 nm	10 000 dt	93 %	22 %
		100 000 dt	93 %	21 %
	15 nm	100 dt	63 %	23 %

Table 5 -

 5 Summary of the sputerring process occuring for a PKA of 100 keV in a NP with a size of 5 or 15 nm for different water thermalization velocities

	NP size	Thermalization velocity	PKA escaping the NP (over all runs)	PKA escaping the NP solely (over all runs)
		100 dt	100 %	42 %
	5 nm	10 000 dt	100 %	23 %
		100 000 dt	100 %	23 %
	15 nm	100 dt	96 %	40 %

Table 6

 6

gathers all the results, highlighting the average value on all runs of the 100 Mo/ 99 Mo ratios i.e. the total number of Journal XX (XXXX) XXXXXX Panetier et al 9

Table 6 -

 6 Mean number of 100 Mo/ 99 Mo resulting of the sputerring process occuring for different parameters

	NP size	Thermalization velocity	19 keV PKA	100 keV PKA
		100 dt	3.62	1.76
	5 nm	10 000 dt 100 000 dt	10.22 12.65	5.62 5.73
		vacuum	43.21	18.87
	15 nm	100 dt vaccum	4.28 50.81	3.85 29.84

Table 7 -

 7 Summary of the values corresponding to the 95 th percentile of the sputtered atoms energy distribution and corresponding projected range (Rp) in water

	NP size	Thermalization velocity	Energy	Rp in water
		100 dt	16.9 keV	31.5 nm
	5 nm	10 000 dt	4.2 keV	14.2 nm
		100 000 dt	3.5 keV	13 nm
	15 nm	100 dt	15.8 keV	30.3 nm

Acknowledgements

The authors would like to acknowledge the Proof of Concept Program of the Joint Research Centre for financial support to the project M-Iso-Sep (Modelling of Isotope Separation by irradiation). The authors also thank M. Mommey and Y. Giraud (IP2I Lyon, France) for their help with the HPC, also J. Amodeo for fruitful discussions about NP behaviour and P. Serra Crespo for insights about Mo-99 production.