Revisiting the effect of shear stress on the γ฀α phase transition of Cerium under shock loading

The dynamic response of Cerium under high pressure and temperature is not clear due to the complexity of its phase diagram. Using large-scale atomistic simulations, we report the γ→α phase transition (PT) of Ce under shock loading. The PT behaviors predicted by present simulations are in good agreement with previous experiments, thus confirm the validity of the interatomic potential. While the simulated wave velocities of elastic precursor are larger than experiments, the PT wave velocities agree well with the LASL data. In the T - P phase diagram, the Hugoniot states generally agree with the lower-pressure extrapolation of the higher-pressure experimental data, but differ significantly from each other orientations. By examining the ratio between the shear stress and hydrostatic stress, it is found that the anisotropic response of PT kinetics might be dependent on the variance of shear stress with lattice orientations. The present study suggests that the effect of shear stress usually ignored in shock wave experiments has to be treated seriously when studying the low-pressure γ→α PT of Ce.


Introduction
Since discovered in 1803 simultaneously in Sweden and Germany, Cerium is one of the most amazing elements in the periodic table due to its complex phase diagram, which contains up to seven allotropic phases (γ, β, α, αʹ, αʺ, η, δ) [1][2][3][4][5][6][7][8].Meanwhile, due to the similarities such as crystal structure, melting temperature and phase transition (PT) behaviors, Ce is considered as a non-radioactive surrogate of δ-Pu alloys [9,10].Among the PTs of elemental Ce, the most intriguing is the isostructural γ→α PT, which is commonly believed to be associated with the variation of 4f electrons' states [1,5,[11][12][13][14].Despite numerous models [11][12][13][15][16][17], e.g. the Mott transition model [11] and Kondo volume collapse (KVC) model [12,13], have been proposed to describe the electronic dynamics during the γ→α PT, none of which can explain the published experiments [18][19][20][21][22] in one unified framework.The Mott transition model [11,23] hypothesizes the PT is dependent on the competition between the intersite f-f hybridization and the on-site f-f Coulomb interaction [24], i.e. f electrons are localized in γ phase but delocalized in α phase.The KVC model assumes f electrons hybridize with the spd electrons, and are localized in both γ and α phases [12,13].The triggering of PT in the KVC model [13] is associated with the change in the screening of the localized moments caused by the conduction spd electrons induced by the volume change.Detailed reviews of theoretical models, previous experiments and simulations can be found in Refs.[1,24,25].
During the low-pressure PT of Ce, the stress deviator is not negligible compared with the hydrostatic pressure, the variation of entropy and temperature induced by plastic deformation would contribute to the PT process [26].Thus, the evolution of shear stress has to be considered when studying the γ→α PT of Ce.Besides, Jeong et al. [27] estimated the lattice entropy change is about half of the total entropy change during the γ→α PT.Despite there have been a number of computational studies covering the melting curve [28] and γ→α PT [5,24,[29][30][31], these density function theory based simulations main focused on the evolution of electronic states at zero temperature and did not offer a glimpse of the lattice contributions and shear-induced deformation.Although classical molecular dynamics (MD) simulation could capture the lattice vibration at finite temperatures, few studies have been reported due to the difficulty of generating potentials that can describe the electronic picture precisely.
While static experiments [32][33][34][35][36][37][38] have identified 4 phases at zero pressure, several more solid phases at higher pressures, dynamic data are limited to the high-pressure liquid Hugoniot [39] and low-pressure melting [38] due to the difficulty of obtaining off-Hugoniot data over the wide range of pressures and temperatures.Jensen and Cherne [3] examined the dynamic response of Ce along the loading paths that span the γ→α PT, and determined the volume compression along the phase boundary pointing to a critical point at 1.648 ± 0.075 GPa.Recently, Hixson et al. [40] measured the shock temperature, stress and dynamic emissivity of Ce shocked from 8.4 to 23.5 GPa, and found that the samples were shocked from the γ into α phase when the shock stress is below 10.24 GPa.Zhernokletov et al. [41] determined the γ→α PT pressure as 0.83 ± 0.06 GPa among a number of shock wave tests ranging up to 3 GPa.The above review of reported experiments thus indicate that the more shock wave experiments are needed to provide enough data for constructing a multiphase EOS of Ce.
Considering the efficiency to control loading conditions, the present study uses the embedded atom method (EAM) potential formulated by Cadien et al. [28], to study the shock response of Ce via large-scale atomistic simulations.The results show that the isostructural PT observed in experiments can be recovered with the f electron state not changed, thus support the validation of the KVC mechanism.The anisotropy of shock response is attributed to the dependence of shear stress on the lattice orientations.

Methods
We performed MD simulations of the shock wave profiles in both single-crystal and nanocrystalline Ce undergoing piston-driven shock loading.All simulations were conducted by the open source code LAMMPS [42].As shown in Fig. 1 and Table 1, 5 types of single crystal samples with different lattice orientations were constructed.Based on extensive ab-initio calculations, the atomic interactions of single element are described by the EAM potential [28], which was obtained by fitting the potential energy surfaces of trivalent Ce (5s 2 5p 6 6s 2 5d 1 4f 0 ) and tetravalent Ce (5s 2 5p 6 6s 2 5d 1 4f 1 ) as well as their pseudo-alloys (Ce 1−  0 Ce   1 ).Previous studies [28] have shown that the melting behavior and lattice parameters yielded by tetravalent valence state EAM potential (Ce  1 ) are in good agreement with experimental values under high pressures.As inspired by the KVC model, in which the 4f electrons remain localized in both the γ and α phases, we use the f-electron localized (trivalent) potential (5s 2 5p 6 6s 2 5d 1 4f 0 ) in all present simulations.
The time step for the velocity-Verlet integration was set as 1 fs.Initially, the energy was minimized using the conjugate gradient algorithm.Then, the system was equilibrated within the microcanonical (NVE) ensemble, simultaneously with rescaling its global temperature to 300 K (if the initial temperature was not specified).Sequentially, the boundary condition along the x-axis direction was changed into shrink-wrapping, while keeping the y-and z-axis directions periodic.The shock loading was applied by assigning the "piston" an impact velocity u p (i.e. the particle velocity defined in experiments) along the x-axis direction (see Fig. 1).To output the profile of stress, temperature, density and particle velocity, the specimen was sliced into a series of bins along the x-axis direction.The dislocation extraction algorithm [43] was applied to identify dislocations generated during the loading process.The adaptive common neighbor analysis (a-CNA) pattern [44] and centro-symmetry parameters [45] were calculated for post-processing.The defective structures were visualized in OVITO code [46].

density profile
Initially, we calibrated the EAM potential [28] by comparing the simulation results to previous experiments [3] of elemental Ce and our recent shock wave experiments of Ce-5wt.%Laalloy [47].Fig. 2a shows the density profiles of Sam-SC-#1 sample simulated with particle velocities and initial temperatures corresponding to the experiments of elemental Ce [3].Consistent with experimental results, a sharp transition of density can be observed after the shock front in all six tests.The corresponding snapshots shown in Fig. S3 indicate that the lattice structure after the shock front is still of fcc type, thus the transition is the isostructural γ→α PT.Fig. 2b further indicates that the occurrence of the isostructural PT is accompanied with the emerging of the shoulder feature at the second peak of the radial distribution function (RDF).More comparisons with shock wave experiments can be found in Fig. S1 & S2 (see Supplementary Information).
The agreement between MD simulations and shock wave experiments [3] evidences that the electronic picture is well described by the f-electron localized (trivalent) potential [28] during dynamic compression.Fig. 3 shows the density profile of Ce loaded with various particle velocities corresponding to recent experiments of Ce-5wt.%Laalloy [47].As argued by Wang et al. [48], the La alloying should not change the electronic structure and lattice dynamics of Ce matrix significantly, due to the reasons that: 1) without f electron occupancy, La acts only as non-magnetic ion, which could tune the concentration of localized f electron of Ce; 2) La alloying could still maintain the fcc structure of Ce matrix.Considering the lack of Ce-La many-body potential, it is acceptable to model the Ce-5wt.%Laalloy with the Ce EAM potential due to the very few additions of La.For tests with particle velocity equal to 0.263 and 0.319 km/s (see Fig. 3b, c), sharp transition of density can be observed for all curves drawn at different moments, while the transition is not significant for the test with particle (piston) velocity of 0.087 km/s (see Fig. 3a).Fig. 3d shows that the wave profile is more complex than two-wave structure since the dislocation plasticity was activated under 1.045 km/s particle velocity loading.As shown in Fig. 4, we further choose two distinct regions with different density at t = 70 ps for corresponding cases shown in Fig. 3 to calculate the radial distribution function (RDF).The results indicate that for tests with the particle velocity equal to 0.263 and 0.319 km/s, a shoulder feature emerges on the second peak (at about 4.90 Å, indicated by the horizontal blue arrow) of the RDF of 700 < d < 800 Å region, where the isostructural transition occurs.In the study of shear-induced plasticity of bcc Ta under high pressure, Wu et al. [49] attributed this shoulder feature to the 1D plastic flow induced by shear and heating.However, as shown in Fig. 4c, both the common neighbor analysis (CNA) and dislocation extraction algorithm measurement show that no dislocations are generated during shock loading when the particle velocity is lower than 0.319 km/s.The MD results are thus consistent with the experiments, and clearly show that the γ→α isostructural PT occurs in the region after the shock front.As mentioned above, although the slight signal of density change was also captured for 0.086 km/s loading in Fig. 3a, the RDF calculation shown in Figure 8a does not support the existence of the isostructural PT.While significant variation of density can be observed in Fig. 3d, the RDF of 1.045 km/s shot (see Fig. 4d) shows that no PT was activated.Fig. 5 exhibits density profiles of all 7 shots reported in our recent paper [47], and corresponding atomic configurations at t = 70 ps.The simulation results agree with experiments, and show that the isostructural PT was only activated for tests with piston velocity equal to 0.263 and 0.319 km/s (i.e. the shot.5 and 6 in experiments [47]).The agreement between simulations of elemental Ce and experiments of Ce-5wt.%Laalloy might indicate that the La alloying does not change the electronic structure of Ce significantly.Besides, it is noted that the γ→α PT can be realized by the same f-electron localized potential, which is consistent with the KVC model and recent APRES measurement [8].Although the PT behavior can be qualitatively recovered by atomistic simulations, the difference of shock wave velocity (see Fig. 6 in following section) between present MD results and our recent experiments of Ce-5wt.%Laalloy [47] still indicates that the La alloying would affect the γ→α PT in thermodynamic space.γ→α PT was recovered again for the tests of Sam-SC-#2 sample with piston velocity equal to 0.263 and 0.319 km/s (see Supplementary information).The results thus indicate that the activation of isostructural PT is strongly dependent on the lattice orientation.Specifically, to relieve the locally high stress concentration induced by flyer impact, the specimen deforms via either the PT or dislocation multiplication.Since polycrystalline specimens were used in experiments [3,40,50], different lattice orientations are theoretically of equal probabilities to be loaded.

Shock Hugoniot
Fig. 6 shows the MD results of the u s v.s.u p relation, together with the experiments by Jensen and Cherene [3] and the LASL Shock Hugoniot data [51].The results show that the wave velocity of elastic precursor in single crystal specimens is dependent on the lattice orientation, and is generally larger than experimental values.In the atomistic study of Ti under shock loading, Zong et al. [52] also observed the similar results, i.e. the MD data of shock wave velocity is larger than experimental data and dependent on the loading direction.The simulation results of a nanocrystalline specimen with ~500 nm along the loading direction (~25.8 nm along the transverse direction, 9,643,806 atoms in total) are also larger than experimental values, might be attributed to the fact that the grain size of MD model is much smaller than that of polycrystalline experimental specimens.By fitting the LASL data [51] to the linear relation, we can obtain the fitting parameters as, C 0 = 0.87 km/s and λ = 1.9, which will be used in the following evaluation of Hugoniot temperature T H .It is found that the simulation results of wave velocity are in good agreement with the linear fitting of LASL data.The comparison between MD and experimental data shown in Fig. 6 thus indicates that the evolution of PT dynamics of Ce along the shock Hugoniot can be reasonably reproduced by the EAM potential [28] employed here.Fig. 7 shows the T-P phase diagram of Ce together with its Hugoniot curve and melting curve.
Combined with previous studies [25,28,36,40,50,53,54], the peak temperature (or pressure) data points obtained by present MD simulations mostly fall into the α phase regime.As marked by the dash blue curve in Fig. 7, the Hugoniot temperature can be evaluated via thermodynamic relations as, where the compression coefficient η can be related with the Hugoniot pressure P H , The heat capacity at constant volume can be written as [55], where the lattice vibration contribution is, and the electronic part is, More details about the theoretical formulas can be found in our previous study [55].As shown in the inserted panel of Fig. 7, the theoretical calculation of γ-phase Hugoniot curve is found to underestimate MD results.While most of the simulation data points agree with experimental data [40,50], partial results of Sam-SC-#3, #4 and #5 samples seem to deviate from them significantly.Specifically, at particle velocity u p = 1.045 km/s, the Sam-SC-#3, #4 samples are shocked into melting state, while the Sam-SC-#5 sample enters the mixed αʹ/αʺ regime.The measurements of T-P Hugoniot data further demonstrate that the anisotropic feature of the dynamic response cannot be neglected when studying the low-pressure phase diagram of Ce.

Discussions
As argued by Zong et al. [52], the role of stress deviators or shear stresses is often ignored in studies of shock-induced PT since in most cases the yield strength is much lower than PT stress.However, Pan et al. [26] have pointed out that the stress deviator is not negligible compared with the hydrostatic pressure during the low-pressure PT of Ce.The maximum shear stress during shock loading can be determined as following [56]: where   is the normal stress in the shock propagation direction, and   and   are transverse normal stresses.The effect of shear stress can be characterized by the so-called shear level parameter ϕ [52]: where  ℎ (= (  +   +   ) 3 ⁄ ) is the hydrostatic stress.For example, during the staticpressure-dominated PT in Bi and Fe, ϕ < 5% [52], while ϕ ≈ 30% for InSb and CdS [57][58][59].As shown in Fig. 8-a, b, c, it is found that the ϕ value of the tests loaded along the [251] and [120] lattice orientations (almost reaches ~45%) shows sharp transition if the PT occurs.However, in Fig. 8-d, e, f where no PT occurs, the ϕ value keeps almost constant after the shock front if dislocations were not generated.The results also show that the development of plastic regime leads to drastic decrease of the ϕ value.Revisiting the D-u p relation shown in Fig. 6, it is found that the lattice orientation dependence of shock wave velocity D generally shares the similar trend as that of the shear level parameter ϕ, i.e. the larger the ϕ value, the larger the D value.Specifically, the Sam-SC-#3 samples have the smallest ϕ and D values, while the Sam-SC-#4 and Sam-SC-#5 samples have the largest ϕ and D values.Besides, Fig. 8-c shows that the ϕ value of tests with initial conditions same as the experiments of Jensen and Cherne [3], might decrease with increasing temperature.We thus further perform tests with various initial Cadien et al. [28] Hixson et al. [40] Jensen et al. [50] Jayaraman et al. [54] Jayaraman et al. [54] melting curve [28] temperatures  0 (but all the same particle velocity 0.319 km/s) to evaluate the variation of ϕ as a function of  0 .Despite Koskimaki et al. [60] has pointed out that the single phase hP4 Ce can be made by thermally cycling the Ce sample between the room temperature and 4 K followed by annealing for long time at 348 K [25], the martensitic transformation from γ to β Ce actually depends on the purity and grain size of samples, and proceeds very slowly, almost never fully finishes.Part of the remained γ phase begins to transform into the α phase at ~ 100 K, while the β phase also starts to transform into the α phase with a temperature between 77 and 43 K [61].However, even after a large number of thermal cycles, neither the γ nor β phase will transform into the α phase upon further cooling.It thus could hypothesize that the γ phase can exist stably under lower temperatures.As shown in Fig. 9-a, a rectangular bin region was selected to detect the evolution of ϕ as a function of time.It indicates that the ϕ value increases drastically at t 1 due to the arrival of shock front, and deceases sharply at t 2 (marked by the red arrow) due to the occurrence of PT.The ϕ value at which the isostructural γ→α PT occurs is thus defined as ϕ c .Fig. 9-b shows that exponential type formulas can be fitted to describe the ϕ c v.s. 0 relation.Specifically, two different formulas, i.e. the three-parameter formula   =  1 +  2  − 0 and two-parameter formula   =  0 + (1 −  0 ) − 0 , where  0 ,  1 ,  2 ,  and  are fitting parameters (shown in Fig. 9-b), were used to fit the (T 0 , ϕ c ) data points.We further studied the ϕ c v.s. 0 relation of the Sam-SC-#1 and Sam-SC-#2 samples under different initial particle velocities, and found that the fitting results of three-parameter model show more dispersity under different loading velocities than that of two-parameter model (see Table 2).Specifically, the slight difference among the fitting values of α parameter indicates that there might exist a constant to describe the temperature effect in the two-parameter model.The fitting values of ϕ c0 parameter in the two-parameter model change significantly with sample type, but only varies slightly with particle velocity.It thus suggests that the high temperature limit of ϕ c required for the γ→α PT (i.e. the parameter  0 ) under shock compression is dependent on the lattice orientation.These results demonstrate that the anisotropy shown in the PT kinetics might be due to the variance of shear stress with lattice orientations.Indeed, the shear-induced anisotropy should not only be considered during the solid-solid PT but also essential to our understanding of the melting of metals under high pressure.Wu et al. [49] found that the bcc Ta undergoes a plasticity-mediated shear-induced structural transition and transforms into an anisotropic viscous plastic flow before becoming a liquid.Given the presence of shear in shock compression, the authors expected that the shock experiments go through a similar transition, e.g. the reported transition of Mo at 210 GPa [62].This argument might be able to explain why only the Sam-SC-#3 and Sam-SC-#4 samples can be shocked into melting state under 1.045 km/s (see Fig. 7).

Conclusion
The dynamic response of elemental Ce under shock compression was investigated by using MD simulations.Three main findings can be concluded as, with f electrons state unchanged, the isostructural γ→α PT found in shock loading experiments can be recovered in simulations of samples with specific lattice orientations; the shock Hugoniot behavior was found to be dependent on the shear level parameter ϕ, thus the lattice orientation of samples; under the same particle velocity, the critical shear level parameter ϕ c required for the γ→α PT decreases exponentially with increasing initial temperature T 0 .It is found that a two-parameter model can be used to fit the ϕ c v.s.T 0 data satisfactorily.The fitting results show that the high temperature limit of ϕ c changes significantly with lattice orientation, but almost not sensitive to the loading particle velocity.
In a word, the dynamic response of single crystal Ce samples is correlated with the shear stress, which is strongly dependent on the loading directions.However, it is not clear how the shear stress affects the evolution of electronic states during PT, which should be paid more attentions in future studies.

Fig. 6 .
Fig. 6.Hugoniot u s v.s.u p relation of Ce.The solid brown line represents the linear fitting of

Fig. 7 .
Fig. 7. T v.s.P phase diagram of Ce.The LDL and HDL represent the low-density liquid and high-density liquid, respectively[28].The inserted panel enlarges the region of [0.0, 7.5] GPa and [0, 750] K to show the present MD results clearly.

2 Fig. 9 .
Fig. 9. (a)A typical showing of the shear level parameter ϕ as a function of the time in a test of Sam-SC-#2 sample loaded under 0.263 km/s at the initial temperature of 500 K.A rectangular region containing several bins was selected to measure the evolution of ϕ.The elastic front arrives this rectangular region at t 1 , while the PT wave arrives at t 2 .(b) The critical shear level parameter ϕ c at which the PT occurs, as a function of the initial temperature  0 for Sam-SC-#1 sample loaded under 0.319 km/s.All tests used to extract the critical parameter ϕ c are performed under various initial temperature, but with other conditions, such as the specimen geometry and particle velocity set as the same.

Table 1 .
Lattice orientations of the Cartesian coordinate in various samples.
LASL Shock Hugoniot data, i.e. u s = 0.87 + 1.9u p in the range of 0.38 < u p < 2.2 km/s.The dash blue curve represents the second order fitting of the LASL data, i.e. u s = 0.768 + 2.133u p -0.113  2 . the