Effects of Re on Vacancy Mobility in a Ni-Re System: An Atomistic Study

+e performance of modern Ni-based superalloys depends critically on the kinetic transport of point defects around solutes such as rhenium. Here, we use atomistic calculations to study the diffusion of vacancy in the low-concentration limit, using the crystalline fcc-framework nickel as a model. On-the-fly kinetic Monte Carlo is combined with an efficient energy-valley search to find energies of saddle points, based on energetics from the embedded atom method. With this technique, we compute the local energy barriers to vacancy hopping, tracer diffusivities, and migration energies of the low-concentration limit of Ni-Re alloys. It was estimated that the computed diffusion rates are comparable to the reported rates. +e presence of Re atoms affects the difference between the energy of the saddle point and the initial energy of point defect hopping. In pure Ni, this difference is about 1 eV, while at 9.66mol% Re, the value is raised to about 1.5 eV. +e vacancy migration energy of vacancy in the 9.66mol % Re sample is raised above that of pure Ni. Our findings demonstrate that even in the low-concentration limit, Re solute atoms continue to play a crucial role in the mobility of the vacancies.


Introduction
Ni-based superalloys are known to exhibit remarkable mechanical properties influenced by the interplay between dislocations and alloying elements [1,2]. Crystalline Nibased superalloys constitute mainly of two phases, c and c ′ , where the c phase is the random substitutional solid solution and the c ′ phase is the Ni 3 Al ordered phase with the L1 2 structure. e c ′ phase is the form of cuboidal precipitates distributed preferentially along the 100 〈 〉 crystallographic direction of the matrix c phase. Recently, a small amount of rhenium (Re) has been added to the superalloys to enhance the strength even further. Particularly, the important role of Re in the c matrix is to drastically increase the creep resistance of the superalloys (the "Re" effect) [3]. Typically and across a broad range of stress and temperature levels, the deformation in the Ni-based superalloys is limited to the plasticity within the thin c channels [1]. e c ′ phase acts as an obstacle to dislocation gliding; hence, dislocation mobility in the superalloys mainly relies on a climbing mechanism known to be associated by the diffusion of vacancies. As well accepted, the simultaneous flux of vacancies toward the dislocation core corresponds to the climb-up motion of the dislocation. In contrast, the flux away from the dislocation core leads to the climbing down of the dislocation. It is known that Re tends to segregate the c matrix and that it does not preferentially form a cluster in the Ni matrix [4][5][6]. e small addition of Re has been anticipated to slow down the vacancy diffusion and ultimately lead to the impediment of the dislocation climbing mechanism [5,7,8]. Recently, it has been found that an isolated substitution of Re for Ni atoms can enhance resistance to the [112](111) slip direction in the c phase [9]. Re is among the slowest-diffusing or immobile species within the c matrix [7,10,11].
is implies that the presence of Re may locally hinder the diffusion of vacancies where the spatial exchanges between chemical species and vacancy are essential.
Based on first-principle calculations and kinetic Monte Carlo (KMC) simulations, Schuwalow et al. [12] recently calculated the vacancy diffusion coefficients with the presence of the Re solutes in Ni. ey have analyzed that, within the dilute limit of Re solutes, the solutes only minimally influence the slowing down of vacancy mobility. However, they suggest that the solute concentration can still be locally increased, especially within the partitioning of the c/c ′ structure. It has been noted that highly concentrated and localized Re solutes are necessary for Re to have a sizable effect on impeding dislocation motion. More recently, Goswami and Mottura [13] considered the nondilute range of Ni-Re systems and studied the vacancy diffusion in Ni using a KMC algorithm and first-principle calculations. ey also employed a cluster expansion method designed for obtaining the total energy of the immense number of atomic configurations at the dilute limit. It was reported that, in the nondilute alloys of Re in Ni, solute-solute interactions have significant effects on the vacancy diffusivity [13].
It should be noted that the previous studies have decisively utilized the KMC methods, rather than conventional molecular dynamics (MD), to study this Re effect. Such a computational approach often ensues in order to overcome the issue of timescales where physical properties of materials owing to phenomena, such as microstructural evolution and creeps, evolve over the experimental timescale (microseconds to seconds). It is well accepted that, though valuable insights can be garnered, the conventional MD methodology lacks the atomistic fidelity capable of reaching laboratorytimescale [14][15][16]. In the past decades, the kinetic Monte Carlo (KMC) or rate theory methods were often employed in the scientific community to overcome this issue [17][18][19]. However, such studies require inputs of atomistic details on the defect interaction mechanism [20], which in turn limits their predictive capability and prevents a faithful understanding of the microstructural evolution and microchemistry, including those of the Re effects.
To overcome this timescale issue, many modern dynamic calculations rely on the concepts of all-atom activation/relaxation techniques based on the concepts of escaping from free-energy minima or metadynamics introduced by Laio and Parrinello [21].
is calculation approach assumes no particular atomistic trajectory in the system evolution, but a few activation parameters. A few examples showing the success of the activation/relaxation techniques include the modeling of long-timescale dynamics in the unfaulting mechanism of trapped atom clusters in bcc Fe [22], void nucleation and growth in bcc Fe [20,23], dislocation-defect interactions and their influence on the strain rate [24], anisotropic diffusion of point defects in hcp Zr [25], and thermally activated deformation in metallic glass [26]. Similar studies have been carried out in other material systems and proved to give new insights into the material kinetics [27][28][29][30][31][32][33].
In this paper, we aim at employing a modified metadynamics algorithm that assesses the vacancy diffusion in Ni with some presence of Re. Unlike the previous works in [12,13], no assumptions about any atomic trajectory are made a priori. e system configurations are autonomously generated based on the concepts of escaping the energy minima during the evolution of the atomic configuration due to diffusion. e employed energy-landscape algorithm for the evolution of vacancies in the Ni-Re systems was the autonomous basin climbing (ABC) method [34]. e ABC method is an algorithm based on the all-atom activation/ relaxation approach that explores and reconstructs the system's potential energy surface. e reaction barriers of the trajectories are then computed via a nudged elastic band (NEB) calculation which connects the transition trajectories. e governing kinetic law at the experimental timescale is then identified using the time interval derived from the kinetic Monte Carlo calculation. Our primary goal is to seek better understanding of the vacancy diffusion mechanism in the Ni-Re system at the rates comparable to those found in the laboratories and real-life applications. e recent works by Goswami and Mottura [13] and Schuwalow et al. [12] precisely focused on this particular alloying system from perspectives that are quite different from ours. us, our findings should present an opportunity for constructive comparisons.

Computational Details
e employed ABC method, developed by Kushima et al. [34] in computing the viscosity of supercooled liquids, is based on Laio and Parrinello's concept of escaping from free-energy minima [21]. e original work by Laio and Parrinello added penalty functions to the free energy of a set of reaction coordinates in an N-particle system in order to force the system to evolve toward new rare-event configurations. e main improvement of the ABC method over that of Laio and Parrinello was that it adds a penalty function to the potential energy of the entire 3N-dimensional space with no knowledge of the reaction coordinates and with no need to acquire a set of forces from local MD simulations. erefore, the ABC method is considered extremely useful in identifying the reaction paths and low-probability configurations without full knowledge of the reaction paths or final configurations a priori.
To study the vacancy motion in the Ni-Re system, a single vacancy site has been placed within the periodic boundary conditions cell. is vacancy migration process was considered an escape from potential basins that push the entire atomic configuration toward new rare-event configurations. A synopsis of our approach and key ingredients to explore the energy landscape for the evolution of the vacancy diffusion is briefly discussed below.
(i) Neighboring-searching with the ABC method: our simulation has employed the ABC method with some extensive steps as in [25]. Neighboring states were autonomously generated from any current atomic configuration by the activation of the configuration on the 3N-dimension potential energy surface. e activation was implemented by using a set of activation functions followed by relaxations in order to sample multiple transition pathways. e devised calculation procedure carries out the following cycle of operations [1]. From a local configuration at an initial minimum r min , a Gaussian activation function ψ 1 (r) � W exp(− (r − r min ) 2 /2σ 2 ) is added to the total energy such that a system energy was Φ 1 p � Φ + ψ 1 . e parameters W and σ are constants that determined the strength and spatial extent of the activation, respectively [2]. From this excited state, the structure is relaxed to minimize Φ 1 p using the BFGS algorithm in order to obtain a new atomic configuration [3]. e cycle is repeated over sufficiently a large number of times until a new local minimum is identified. e total energy of the system is, therefore, progressing as ., depending on the number of cycles undertaken in this step [4]. After a local minimum has been identified, the nudged elastic band (NEB) method as implemented in [35] is introduced in order to collect transition-state configurations (also discussed in the next paragraph). e key in employing the NEB method in this step is for a further application of more penalty Gaussian functions at the determined transition state. e overall atomic system is, therefore, prevented to undergo the same reaction pathway as already determined [25], rendering minimal wasteful computational effort [5]. All the abovementioned processes are cycled until sufficient amounts of minima surrounding the initial minimum have been sampled for transitionstate calculations.
(ii) Determination of transition states: with a suitable ABC algorithm, the nudged elastic band (NEB) method [35] has been used to determine the activated configuration associated with the transition and corresponding activation barriers (E a i ). e subsequential KMC simulation uses this energy barrier as an input. is energy data are converted into transition rates and organized in the form of a look-up transition rates catalog. (iii) Kinetic Monte Carlo: after the ABC sampling and NEB calculations, a set of neighboring minima and the associated energy barriers are known. e kinetic Monte Carlo method was used to push the reaction coordinate forward along one of the available reaction paths to the nearby rare-event states. e probability that such an event would occur within the time t is associated with exponential distribution [19,36], where k tot is the total rate of migration events that can occur. It is given by the summation k tot � events i k i of all the rates k i associated with the neighboring minima. e evolution time t is calculated by inverting the exponential probability distribution, p(t) � k tot exp(− k tot t). e synopsis of the devised KMC method can be summarized as follows [1]; from an initial time t � 0, an initial o-state is chosen [2]; the N o possible transition rates r oi are listed as the reaction may progress from o into a generic state i [3]; the cu- [4]; a random number u ∈ (0, 1] is then generated in order to represent the statistical fluctuations of the transitions. e event that sat- All steps are then cycled until sufficient timescale samples have been achieved. (iv) Transition-state theory (TST): with calculated energy barriers from (ii), each transition rate, described in (iii), can be calculated with the harmonic approximation of the TST as where k * is an attempt frequency. As an approximation, we used the value of k * � 10 13 s − 1 . e time update (t + Δt) in step (iii), in accordance with principles in probability theory, is expressed as where r ∈ (0, 1) is another random number drawn from a uniform distribution. (v) Diffusion coefficient determination: based on the random-walk theory, the vacancy diffusivity coefficient can be calculated using the mean square displacement of a vacancy position and time information in (iii) and (iv). Exploiting the randomwalk theory, the diffusivity coefficient can be expressed as By tracking a displacement of vacancy during KMC simulation in every KMC step, the vacancy diffusion coefficient can be direct-forwardly extracted. To make a direct comparison of the resultant D * with experiments, the KMCbased D * is scaled according to the following expression: where c * vac is the vacancy of the simulation box and c vac, real is the vacancy concentration of the the corresponding real system at the same temperature. is correction is introduced since we only allow for one vacancy within the supercell, making the vacancy content in the simulated system much higher than in real alloys. e reflection back to the physical vacancy concentration and diffusion coefficients is achieved by the abovementioned rescaling factor [18,37]. Since Re is present in low amount, c vac, real can be set to the experimental vacancy concentration in pure Ni. e experiment data determined using differential dilatometry measurements [38] were used to fit the parameter c vac, real as a function of temperature. Based on this, c vac, real � 8.4 × 10 − 6 , 2.67 × 10 − 5 , 7.17 × 10 − 5 , and 1.69 × 10 − 4 are found for T �1200 K, 1300 K, 1400 K, and 1500 K, respectively [18].
In this work, the ABC algorithms were implemented inhouse and on top of the available Atomistic Simulation Environment (ASE) framework [39]. e calculation steps that require the total energy of the systems used ternary (Ni-Al-Re) potential based on the embedded atom potential developed by Du et al. [40,41]. e supercell for the atomistic simulation was cubic with 2a o × 2a o × 2a o dimensions and containing 32 lattice sites in the fcc structure. Here, the equilibrium lattice parameter was a o � 3.542Å.
e Ni-Re alloying contents of 3.22 mol% Re, 6.44 mol% Re, and 9.66 mol% Re were represented by putting one, two, and three Re atoms (and one vacancy) in the supercell of 32 lattice sites, respectively. In NEB steps, relaxation was performed until forces from all images were below 0.1 eV/Å. e parallelized version of this algorithm in the LAMMPS software [42] has been employed to facilitate this NEB calculation.

Transition State of the Vacancy Hopping Process.
e fcc Ni is an anisotropic material in which there are multiple adjacent locations on which a vacancy can migrate. It is believed that the presence of a Re solute atom can influence dislocation climbing as it resists the migration of vacancy. e underlying vacancy evolution is the positional exchange between a vacancy position and its neighbor atoms, as shown in Figure 1.
e close-packed stacking of the fcc structure with a vacancy at the center is shown in the figure. Note that, in Figure 1, a basal plane (111) is shown along with top-stacking and bottom-stacking planes.
ere are two migration paths for the vacancy diffusion, one on the basal (111) plane and one out of the basal plane. For example, the exchange of the atom/vacancy pairs vacancy ⟶ 8 and vacancy ⟶ 7 constitute a basal plane diffusion, while the exchange vacancy ⟶ 12 constitutes a diffusion out of a basal plane.
To investigate the effect of rhenium solute on the vacancy diffusion, a Re atom has been placed at the 9th location in Figure 1, while the centered vacancy is allowed to migrate freely. Diffusion energy barriers for the first neighbor shells in the proximity of the centered vacancy are summarized in Table 1. In the table, the migration mechanisms and corresponding energy barriers of point defects were identified by our ABC method. For example, in the case of vacancy migration in pure Ni, i.e., self-diffusion, our energy barrier (0.98 eV) is consistent with other reported values obtained using other computational approaches (1.01 eV in [12], 0.997 eV in [43], 1.22 eV in [18], 1.2 eV in [7], and 1.04 eV in [44]). e overall results revealed direct quantitative comparisons of the activation energies of point defects migration in the pure Ni and Ni-Re systems.
In the pure Ni system, the energy barriers acquire the same values across all in-plane and out-of-plane transitional paths. However, with a Re atom added to the 9th location in Figure 1, the ABC algorithm reveals that the vacancy finds itself harder to migrate to the 9th location where the Re atom is located. However, there are lower energy barriers for outof-plane migration of vacancy toward the 1st and 12th locations. Note that the higher-energy transition (vacancy ⟶ 9) has been found within 4 atomic activation/ relaxation steps within the ABC calculation, as compared to 25-35 ABC steps for the out-of-plane vacancy migration.
is, in part, indicates that there are some intermediate small transition steps that the resolution of the ABC steps may not be able to catch. is issue is often found during the simulations at low temperatures and identified as a "sandtrap" limitation [45], where chemical species spend a large number of steps moving locally without accomplishing significant displacements. Hopping back and forth, the vacancy in the "sandtrap" is confined between multiple valleys. e saddle point of the energy surface between these small valleys is much smaller than that necessary for the escape with longer displacements. Small activation energy in the ABC method is needed to capture this behavior, although the cost of computational time renders the simulation inefficient. We believe that the origin of the "sandtrap" is physical and due to interactions among the point defects, Re, and Ni atoms. Hence, no attempt to forcefully taking the vacancy out of traps has been made in our ABC calculations of the effective migration barriers.
Due to this "sandtrap" issue, other nonequivalence mechanisms (i.e., vacancy ⟶ 4, vacancy ⟶ 5, vacancy ⟶ 6, vacancy ⟶ 3, and vacancy ⟶ 2) are not reachable by the present ABC approach. However, as comparisons, we have shown energy barriers of these nonequivalence mechanisms calculated from classical NEB methods in Table 1 (values in parentheses).
ese NEB calculations were made assuming the final atomic structure of all exchanges is that of the vacancy ⟶ 1 mechanism. Since the final structures of these calculations were not fully relaxed, the NEB values were reported only as complementary results. It can be seen that the motion of vacancy, both in-plane and out-of-plane, is made harder with the presence of a nearby Re solute. e exchange between the Re atom and the point defected is further highlighted in Figure 2. e detailed pathway was obtained for an NEB calculation. e blue sphere represents the Re atom, while the vacancy is represented as the smaller pinkish sphere. Figure 2 shows the energies of all NEB images, including the first replica as shown in Figure 2 e transition-state configuration corresponds to the state with the higher energy where the Re atom is in between its original site and vacancy site. As alluded in a previous paragraph, the Re center enhances the energy barrier of the Re-vacancy exchange, leading to the overall reduction of point defect diffusivity. e influence of rhenium atoms to the diffusion coefficients will be discussed in the following sections.
By the ABC method, energy differences in each step during the atomic hopping can be ordered, as shown in Figure 3. It can be seen that the general trend of increasing difference between the initial energies before the atomic jump and the energies of saddle points can be seen. In pure Ni, this difference is about 1 eV, while at 9.66 mol% Re, the value is raised to about 1.5 eV. is finding is consistent with the data in Table 1, where the presence of Re is associated with the higher energy barrier. e overall energy barrier is raised by the presence of Re. An example of an atomic trajectory of a vacancy and Re atoms is shown in Figure 4. It can be seen that, in the course of the simulation, Re atoms rarely moved while the vacancy was quite mobile. It should also be noted that the location of the vacancy was found throughout the lattice, hinting that no net attraction/repulsion of vacancies by individual solute atoms was anticipated.

Effect of Re Concentration on Vacancy Diffusion
. By using the technique described in the methodology section, we calculated the bulk tracer diffusivity in a Ni-Re alloy in lowconcentration limits. Four difference values of Re concentrations were picked in this investigation, namely, the pure Ni or 0 mol% Re, 3.22 mol% Re, 6.44 mol% Re, and 9.66 mol % Re. We deliberately excluded the extremely high Re composition in this investigation to avoid the discussion on the Re solute-solute interactions, which is not the primary objective of this work. e diffusion prefactor (D o ) and the migration energy (Q vac ) have been evaluated by fitting the numerical results to the following Arrhenius relation: As a comparison and validation of our approach, Figure 5 displays the computed Ni self-diffusion coefficients in the temperature range of 1200-1500 K against available experimental [46,47] and simulation results [18,43]. At the temperature range considered, it can be seen that the predicted vacancy diffusion coefficients yields the effective energy barrier in the same order as those reported in [18,43,46,47]. On the other hand, our diffusion prefactor (D o ) is an order of magnitude higher than first-principle results.
is difference could stem from either (1) the electronic interactions pertained to the density functional theory used in the references [18,43], while in the present work, an embedded atom potential is utilized, or (2) the difference in state sampling since in our calculations, all atoms are fully relaxed during the course of the ABC algorithm. However, as seen from extrapolating curves labeled Table 1: Elementary vacancy diffusion pathways observed by the ABC method and their corresponding energy barriers. e position numbers are depicted in Figure 1. e corresponding energy barriers calculated using classical NEB methods are reported in parentheses in the unit of electron volt. e values in the last column are taken from a [12], b [43], and c [7]. Ni (Present) and Bakker (Exp.) in Figure 5, the predicted effective energy barrier Q vac and prefactor D o are in good agreement with the experimental data in [47]. Numerical values emphasizing this latter point are tabulated in Table 2.
In Figure 6, the calculated diffusivities in the investigated compositions of 3.22 mol% Re, 6.44 mol% Re, and 9.66 mol% Re are displayed in comparison to the diffusivities of vacancy in pure Ni. Generally, it can be observed that the diffusivity tends to decrease with increasing Re content. In contrast to a recent report [12], our findings show that, within the low-Re content limit, even a few Re solutes may affect the mobility of vacancies as previously theorized in [5,6]. ough small in value, the reduction in diffusion coefficient can be seen more clearly at lower temperature ranges (1300-1200 K). us, the influence of Re is expected to be more pronounced in the c channel within the c/c ′ structure where partitioning of Re atoms to form higher-concentration Re clusters has been observed. As alluded earlier, the slight changes in the vacancy diffusion coefficient are due to the intrinsic Ni mobility in the lattice. In the presence of Re solute that moves more slowly within the lattice than Ni atoms, the resultant vacancy mobility will be low. With complementing observations from Figure 4, the diffusion retardation mechanism agrees with the picture of simple-point obstacles, where the nearing-zero interaction cross section between Re solutes and vacancies is observed.
To see the spreading of the diffusion coefficients, we collected a set of diffusion coefficients, computed from   equation (3), in the last 50 out of the total 100 KMC steps. As a result, a total of 50 diffusion coefficients are shown in the form of histograms in Figure 7. As is evident from Figure 7, the estimated diffusion coefficient at low temperature exhibits a relatively narrow range of values: within 1.72 × 10 − 14 and 1.94 × 10 − 14 m 2 s − 1 at T � 1200 K. At the temperature T � 1300 K, the range slightly expands to 1.0 × 10 − 13 − 1.14 × 10 − 13 m 2 s − 1 . e distribution range of the case T � 1400 K is similar to that of T � 1300 K. e range is the widest at T � 1500K, where it spans from 1.78 × 10 − 12 to 2.06 × 10 − 12 m 2 s − 1 . At higher temperatures (T � 1300 − 1500 K), the preponderance of the distribution (about 70% of the number of diffusion coefficients) is in a narrow range around the modes. It should be noted that, in all distributions, the modes shift to higher values as the temperature increases, reflecting the higher vacancy mobility at higher temperatures. e overall distribution of diffusion coefficients, therefore, show no unusual behavior where irregularly fast or slow diffusivities are identified.
Effective activation energy barriers Q vac and diffusion prefactors D o are summarized in Table 2, where the increases in both Q vac and D o with higher Re content are revealed. Using the energy barrier for self-diffusion of Ni as a reference, adding 9.66 mol% Re reduces Q vac albeit by a small amount (about 5% when compared to Q vac of pure Ni).

Alfonso Zhang
Figure 5: Self-diffusion coefficients for pure Ni at T �1200 K, 1300 K, 1400 K, and 1500 K and its comparison with available experimental data [46,47]. Available results derived from first-principle calculations are also shown [18,43].   Journal of Chemistry However, it transpires that Re is effective at reducing vacancy mobility in the c phase of Ni-based superalloys. e change in the energy barrier is also shown in Figure 8. e continual rising of Q vac is anticipated at the higher Re content beyond the limiting value (9.66 mol% Re) in this study. Considering the change in the diffusion prefactors, the increasing diffusion prefactors, together with the increasing Q vac results in the reduction of overall diffusion coefficients in the temperature range 1200-1500 K considered. It should be noted that the calculations presented in this work do not consider the full effects of vibrational entropy, and hence, we do not expect the change in hightemperature intercepts of Figure 6, across all considered Re contents, to represent true results. Lastly, it should also be noted that the trend of increasing Q vac with Re content (shown in Figure 8) may only be a characteristic of a lowconcentration limit of the Ni-Re system. A fuller investigation at higher concentrations should require the consideration of solute-solute interactions within the KMC model.

Conclusions
e presence of Re is critical to the performance of Ni-based superalloys. In this work, atomistic calculations based on an activated saddle point search method, namely, the ABC method, was utilized to study the diffusion of vacancy in the low-concentration limit. e energetic terms are obtained through an embedded atom potential, while the material kinetic was represented by the kinetic Monte Carlo process. By examining the vacancy diffusion in pure Ni, it was argued that the computed diffusion coefficients are comparable to available experimental data and data obtained from firstprinciple calculations. e presence of Re atoms affects the difference between the energy of the saddle point and the initial energy of point defect hopping. In pure Ni, this difference is about 1 eV, while at 9.66 mol% Re, the value was raised to about 1.5 eV. e effect of Re atoms on vacancy mobility becomes more pronounced at lower temperatures. When considering the corresponding Arrhenius relation to the temperature dependence of the diffusion coefficients, the effective energy barrier of vacancy in the Ni with 9.66 mol % Re was raised by 5% percent over that of pure Ni. In contrast  to a recent report, we argued that the presence of Re solutes, even in small amounts, aids the reduction of vacancy diffusion. e trend of increasing effective energy barriers Q vac is anticipated at the higher Re contents beyond values (0 to 9.66 mol% Re) in this study.
Data Availability e datasets of this research are available and will be provided upon request.

Conflicts of Interest
e authors declare no conflicts of interest.