Molecular Dynamics Simulation of the Clustering of Minor Lead Additives in Liquid Sodium

A strong influence of minor lead additives on the liquid sodium microstructure is revealed in the molecular dynamics (MD) simulation of the Na0.98Pb0.02 alloy. The obtained results can be explained by the existence of lead-sodium clusters in liquid sodium built up by ionic bonds, Na+–Pb−, due to essential distinction of the alloy components in the electronegativity. On this reason, MD simulation of the Na0.98Pb0.02 alloy is carried out within the framework of a three-component bipolar model, Na + Na+ + Pb −, with Na↔Na+ recharging the nearest-neighbor particles of solvent in every 3 ps (an optimal period) during the numerical run.


Introduction
A nonideal solution of lead in liquid sodium, as an intermediate element in the alkaline-metal group, takes place in the whole composition range of Na-Pb alloy, but the strongest changes of thermodynamic properties are observed in the neighborhood of 20 and 50% at. of lead concentration [1][2][3].It is connected with the presence of two kinds of clusters in the alloy: (Na 4 Pb) n and (NaPb) m , having an ionic bond, Na + -Pb − , in diluted solutions and the covalent one, Na÷Pb, for the other compositions.They have explained this by the existence of Zintl's cluster, (Na + ) 4 (Pb 4 ) 4− [3][4][5][6][7].The anionic tetrahedron, (Pb 4 ) 4− , is considered as a united particle with sodium cations, Na + , on the four faces of this tetrahedron.
However such anion is not in a steady state in liquid sodium [6,7].Therefore the application of Zintl's cluster model [5] for understanding properties of the Na-Pb alloy is not enough.A new MD study as well as the experimental one of its microstructure and atomic dynamics by the method of neutron scattering is necessary.The knowledge about the Na-Pb alloy properties can be useful to develop a concept for improving the composition of the sodium coolant in particular by reducing its chemical activity in the environment and an automatic shut-down of sodium fires when sodium is used as a coolant of a fast nuclear reactor [8].

The Cluster Model of Na-Pb Alloy
It is known [9,10] that the mixture of sodium and lead relates to strongly interacting systems with covalent bonds that are illustrated by a long homologous series of compounds, Na n Pb m , where n and m are integers.They influence the structural and thermodynamic properties of the binary alloy [1-7, 11, 12].This effect can be displayed in analyzing the structural factors, S NN (q), S NC (q), and S CC (q), measured and calculated on various models [2,5,7].It is also characterized by a sharp change of the alloy entropy [1] or an increase of the alloy electric resistance [3,12] at x Pb ∼ 0.25 as well as a deviation of alloy composition fluctuations from the ideal mixture shown in Figure 1 as a function of S CC (q) q → 0 depending on the composition of the Na-Pb alloy [2].
One can see that the structural properties of Na-Pb alloy have distinctive features in the neighborhood of x Pb , which is equal to 0.2 and 0.5.This confirms the transitive character of sodium from easy lithium to heavy alkaline metals, Cs and K, and the existence of different clusters in the Na-Pb alloy such as (Na 4 Pb) k , (NaPb) l , and (NaPb 4 ) m .In that case, the cluster model of solid spheres [2,11] is the most realistic for considering the Na-Pb alloy in the whole range of composition.Moreover, the structural factor S(q) of Na 0.8 Pb 0.2 alloy calculated on the solid sphere model of Takeda et al. [11] or Hafner et al. [9] agrees with the The zero limit of structural factor, S CC (q → 0), as a function of the Na-Pb composition (solid line) at 803 K (1) and 698 K (2), the same limit for the ideal mixture of sodium and lead (dashed line).The structural factor, S(q), of the Na 0.8 Pb 0.2 alloy as an experimental result (line in black) [11] and the calculated ones on Hafner's model (line in red) at 698 K [9] and Takeda's model (line in blue) at 703 K [11].
experimental data [11] obtained by scattering cold neutrons in samples of this alloy at 703 K (see Figure 2).Thus, the structural factor of Na 0.8 Pb 0.2 alloy has a subpeak in the neighborhood of q equal to 1.2 Å which is caused by spatial correlation of clusters in the alloy.The application of indirect Monte-Carlo method [14] confirms the presence of spatial microheterogeneity of Na 0.8 Pb 0.2 alloy.
MD investigations of minor additives of lead in liquid sodium (x Pb < 0.2) are a few.One can mention the calculated structural factor, S(q), of Na 0.9 Pb 0.1 alloy obtained by Hafner et al. [9] and the experimental one obtained at neutron scattering [15] in the sodium alloy with 1.5 and 7.9% at. of lead.
At the same time, the character of heterogeneous structure of Na-Pb alloy remains unclear.Therefore we have applied a "partly ionized alloy" model [13,16] for MD simulation of the Na-Pb alloy as lead anions and sodium cations in the liquid matrix of "neutral" sodium because of the approbation of this model proves to be successful for describing such systems.

MD Model of Na-Pb Alloy
There is a technique used for MD simulation of the Na 1−x Pb x alloy of essentially different electronegative components which is based on the assumption that the additive one (Pb) is completely ionized at the electric neutrality of the system as a whole.This model describes a three-component system consisting of lead anions (Pb − ), equal number of sodium cations (Na + ), and neutral sodium atoms (Na) so that the total number N of particles in a MD cell is equal to Six pair potentials are necessary for describing the interaction of three types of particles: U NaNa , U NaNa(+) , U NaPb(−) , U Na(+)Na(+) , U Na(+)Pb(−) , and U Pb(−)Pb(−) .We accept that U NaNa(+) is equal to U NaNa which is the pair potential for the interaction of sodium particles in one-component system.For ionic pair potentials as U Na(+)Na(+) , U Na(+)Pb(−) , and U Pb(−)Pb(−) , the Born-Mayer-Huggins potential is used.It is applied earlier for simulating Na 0.5 Pb 0.5 [5] and used here for describing the interaction of neutral sodium particle and lead anion; U NaPb(−) for a charge of sodium is equal to zero.
In the offered model, the charge of sodium cation can be "bound" rigidly to the particle (Na + ) or handed to neighbor sodium atom by using the procedure of charge exchange among the nearest-neighbor particles (Na + ↔Na).This procedure can be a periodic or adaptive time function [17] for relaxing the charge perturbation of the system.This option allows approaching the three-component model of the Na-Pb alloy to a real electron exchange in the considered system.

Techniques for Numerical Runs
The MD simulation of Na 1−x Pb x alloy (where x = 0.20, 0.09, 0.02) is carried out by MDMC code [17] in the frameworks of NVT ensemble of particles when all the lead atoms and the equal number of sodium atoms are oppositely ionized (Pb − , Na + ) so that the alloy remains neutral as a whole.The periodic boundary conditions are used for a MD cube, and the interaction of particles is cut off in half of its edge.Ewald's summation is used for calculating the ionic pair potentials in the approximation [5].The equations of motion are solved by using Verlet's algorithm in the velocity form [18] and a time step of 2•10 −15 s.A numerical run is begun from random packing of particles in the MD cell, stopped at a given point of time for calculating required characteristics, and continued from the last atomic configuration.The temperature, pressure, kinetic and potential energy of the Figure 4: The structural factor, S(q), of the Na 0.8 Pb 0.2 alloy in the experiment (line in black) and numerical runs on models of Takeda (line in blue) [11] and Hafner (line in red) [9] as well as in numerical runs on the partly-ionized-alloy model [13] at the different values of Ω = 247, 260, 266 a.u.. system, and their numerical error bars have been defined every 50 time steps.The root mean square deviations of energy and pressure did not exceed 1%.The results are obtained by using Hasegawa's potential [19] for U NaNa = U NaNa(+) with a dielectric function [20], where R C is equal to 1.85 a.u.For the Born-Mayer-Huggins potentials: U Na(+)Na(+) , U Na(+)Pb(−) , U Pb(−)Pb(−) , and U NaPb(−) , the next ion diameters are used [5]: D Na(+)Na(+) = 3.77 Å, The Na 0.8 Pb 0.2 alloy is simulated in MD cells with the molar volume of Ω ∼ 247 a.u.[15] or Ω ∼ 260 a.u., and the Na 0.98 Pb 0.02 alloy is simulated in the MD cube with Ω ∼ 297 a.u.(the edge length is equal to 76.13 Å).
The numerical experiment included three stages: (1) generating the random particle configuration at the beginning of the run, (2) putting the MD model of system in thermodynamic equilibrium, (3) calculating some characteristics of particle dynamics and fixing snapshots of atomic configurations of the system in points of time from 2 up to 100 ps in every 2 ps during the run.
Each run taken separately is carried out with or without the option of charge exchange (see above) every 3 or 4 ps.In the first case (with the option of charge exchange), the system is balanced every time as it is shown in Figure 3 which demonstrates a total energy change of the Na-Pb alloy after which the charge exchange is carried out.It is visible that the MD simulation with the option of charge exchange of sodium particles in Na 0.98 Pb 0.02 alloy deduces the system in the thermodynamic balance earlier and at higher level of its energy.It specifies the greater configuration disorder of system than it is without the option.

Results of MD Simulation of the Na-Pb Alloy
As marked above, there is the prepeak in the vicinity of q ∼ 1.2 Å in the structural factor, S(q), of the Na 0.8 Pb 0.2 alloy (see Figure 2).This prepeak is caused by the spatial correlation of additive particles in liquid sodium.It proves also by molecular dynamics calculations in different models of the binary system (see Figure 4).It is visible that the results of molecular dynamics simulation of the Na-Pb alloy at 698 K in the model of partly ionized alloy [13] are closer to data of Hafner et al. [9].
On the contrary, such prepeak is absent in the S(q) curve of liquid sodium with the small concentration of lead, x Pb = 0.02, but at x Pb = 0.09, there is only a small bulge (see Figure 5).Its position in the model of partly ionized alloy [13] is in agreement with the same one of Hafner's model [9], but a shift of the main S(q) peak to the left is more correspondent to the data on the neutron scattering [11] in sodium alloy with 1.5 and 7.9% at. of lead.The 2% of lead atoms in sodium reduce diffusion mobility of sodium by 10% without changing the maximum position of vibration spectra.It specifies molecular bonds between ions in liquid sodium, the correlation of their Brownian motion, and a strong dissipation of oscillatory energy, for example, due to diffusive exchange of sodium particles in clusters (Na 2 Pb) n .
It is visible in Figure 3 that the thermodynamic relaxation kinetics of the particle charge exchange in every 3 and 4 ps is the same.Therefore, charge exchange periods less than 3 ps do not change the thermodynamic state of the microheterogeneous Na-Pb alloy.
As mentioned before, the MD simulation of the Na 0.98 Pb 0.02 alloy with the option of recharging sodium particles is established earlier and at higher level of total energy (see Figure 3) and characterized by the greater configuration disorder of the system that is illustrated in Figures 6-10.These runs start with exactly the same initial state.

Discussion of the Obtained Results
One can see in Figures 7, 8, 9, and 10 that recharging sodium particles on simulating the Na 0.98 Pb 0.02 alloy is the important option for adequate description of the clustering of the alloy, and the recharging period should not exceed 3 ps.At the same time, the clustering of the ionic components (Na + , Pb − ) of the alloy illustrated in Figures 6-10 remains as sodium cations and lead anions quickly form chain structure of a percolation cluster which then gradually is transformed in separate cluster groups of different particle density.Here, it is important to note that recharging sodium particles (Na↔Na + ) in every 3 ps and more often only reduces the cluster density but does not change the character of their reconstruction directed to the formation of steady-state colloidal particles (Na 2 Pb) n in a dilute solution of lead in liquid sodium.
These colloidal particles, as marked above, are characterized by molecular bonds between ions in liquid sodium and strong dissipation of oscillatory energy that is caused by an exchange of mobile sodium particles between clusters and the liquid matrix.

Conclusions
Using the three-component bipolar model, Na + Na + + Pb − , with periodic recharging of the nearest-neighbor particles of solvent (Na↔Na + ) has allowed us to carry out moleculardynamic simulation of the Na 0.98 Pb 0.02 alloy and to find (in the numerical runs) the strong influence of the minor additives of lead on the microstructure of liquid sodium.
The mechanism for clustering impurity particles in the dilute solution of lead in sodium is studied.It is shown that the essential distinction of electronegativity of the Na-Pb alloy components is the main reason for forming clusters on the base of ionic bonds, Na + -Pb − .Typical times of the colloid formation in any dilute solution of lead in sodium do not exceed 50 ps.

Figure 1 :
Figure1: The zero limit of structural factor, S CC (q → 0), as a function of the Na-Pb composition (solid line) at 803 K (1) and 698 K (2), the same limit for the ideal mixture of sodium and lead (dashed line).

Figure 2 :
Figure2: The structural factor, S(q), of the Na 0.8 Pb 0.2 alloy as an experimental result (line in black)[11] and the calculated ones on Hafner's model (line in red) at 698 K[9] and Takeda's model (line in blue) at 703 K[11].

Figure 3 :
Figure 3: Change of total energy of Na 0.98 Pb 0.02 alloy without recharging sodium atoms (line in black) and with this option every 3 ps (line in blue) and 4 ps (line in red) during the numerical experiment.

Figure 5 :
Figure5: The structural factor, S(q), of Na 0.9 Pb 0.1 alloy in Hafner's model (line in brown)[9] and in present work for pure sodium (line in blue) and Na 0.91 Pb 0.09 alloy (line in purple) as well as Na 0.98 Pb 0.02 alloy with 201 (line in fuchsia) and 378 (line in green) lead particles in the MD cell.

Figure 6 :Figure 7 :
Figure 6: Ionic configurations in the MD cell of 201 Na + (red balls) and 201 Pb − (blue balls) in 2 ps (a) and 4 ps (b) after the beginning of MD simulation without charge exchange of sodium particles.

Figure 8 :
Figure 8: Ionic configurations in the MD cell of 201 Na + (red balls) and 201 Pb − (blue balls) in 20 ps after the beginning of MD simulation without charge exchange of sodium particles (a) and with it (b).

Figure 9 :
Figure 9: Ionic configurations in the MD cell of 201 Na + (red balls) and 201 Pb − (blue balls) in 50 ps after the beginning of MD simulation without charge exchange of sodium particles (a) and with it (b).

Figure 10 :
Figure 10: Ionic configurations in the MD cell of 201 Na + (red balls) and 201 Pb − (blue balls) in 100 ps after the beginning of MD simulation without charge exchange of sodium particles (a) and with it (b).