Surface and Size Effects on the Behaviors of Point Defects in Irradiated Crystalline Solids

We present an elaborate study of the surface and size effects on the transient and steady-state behaviors of point defects in irradiated solids. In this investigation, both pure Ni and binary Ni-Cr were utilized as model systems. We utilize the spatially-resolved rate-theory (SRRT) modeling approach, and directly account for the effects of dose rate, production bias, and defects recombination, reactions with volumetric sinks, and diffusion to surface sinks. Several simulations were conducted to investigate the effects of these parameters in both coupled and decoupled manners. In the presence of production bias, the effects of surface and size persist even as the surface to volume ratio decreases. This was associated with a surface-induced and size-regulated instability. This instability is only triggered above a critical size between 100 and 500 nm. The critical size decreases with increasing dose rate, increasing production bias, or lowering the temperature. Moreover, this instability results in a pattern that favors the separation of vacancies and interstitials. Once this pattern develops, anomalies in the dependence on size for the transient and steady-state concentrations of point defects and the surface/boundary sink strength are observed. These anomalies tend to render irradiation damage more severe. For pure Ni, it was shown that vacancy supersaturation increases with size, and the rate of increase also rises with size. For the binary Ni-Cr system, it was shown that the magnitude of enrichment/depletion of Ni/Cr at the boundary increases with size, and the width of the enrichment/depletion layer also increases with size. The results obtained here agree well with experimental observations in irradiated materials such as the formation of void denuded zones adjacent to grain boundaries and the size and temperature dependence of the radiation resistance of nanomaterials. The size-dependent behaviors reported here also shed new light on the radiation tolerance of nanomaterials, i.e., the irradiation-induced instabilities are suppressed in such materials. Lastly, the implications of the results obtained here on the development of efficient reduced order models or the utilization of ion irradiation as a surrogate to neutron irradiation are discussed.


INTRODUCTION
Nuclear energy is an important source for the production of emission-free electricity. It has the highest power density and reliability among all forms of energy. Several new reactor designs are being developed to render nuclear energy more efficient, safe, and economic. The main limitation in deploying these novel designs is the qualification of novel materials capable of serving in such harsh environments.
Reactor materials operate under extreme conditions of irradiation, high temperature, high stress, and corrosive media. In such environments, these materials experience microstructural changes that result in the degradation of their physical, mechanical, and thermal properties. Some of the main consequences of microstructural changes are radiation embrittlement, irradiation creep, irradiation growth, void swelling, and high temperature helium embrittlement. Irradiation effects on materials need to be understood to predict their performance and hence facilitate the advancement of nuclear energy applications (English and Hyde, 2003;van der Laan et al., 2012;Ryabikovskaya et al., 2021;Yingling et al., 2021).
Understanding irradiation effects on in-reactor materials is very challenging. First, it is very expensive to run in-reactor experiments. Second, the handling of the resultant radioactive material demands special treatment and equipment. Third, the relatively low dose rate experienced in normal operation requires performing long term experiments that might need to run for years to collect critical data. Fourth, it is extremely difficult to set up reactor experiments to simulate the materials response to transient and abnormal conditions such as in accidents. Fifth, it is also a cumbersome task to separate the combined effects of irradiation, temperature, and stress in such experiments. Lastly, the data obtained from the experiments carried out in a specific reactor might be irrelevant to new reactor designs with different operating conditions. Two main approaches are currently in use to overcome these difficulties. In the first approach, ion irradiation is utilized to mimic neutron irradiation. Ion irradiation experiments can be conducted in ion-beam laboratories and using accelerators. Moreover, because of their high dose rate, accelerated irradiation experiments can be performed in a matter of hours or days (Krasheninnikov and Nordlund, 2010;Was et al., 2014;Shao et al., 2017;Sun et al., 2017;Zhang et al., 2018). Nonetheless, because of this large difference in dose rate, the extrapolation of the ion irradiation data to assess the behavior of materials under neutron irradiation is questionable. The second approach leverages the modern computational tools to implement physics-based and data-driven models to simulate the response of materials to neutron irradiation. Nonetheless, because of the inherently multiscale nature of the problem, conducting high fidelity simulations is still unfeasible.
While the two approaches currently have their limitations, they have demonstrated several successes in elucidating irradiation effects on materials (Was et al., 2002) examined 304SS and 316SS cladding steels samples to evaluate whether proton irradiation can be utilized to simulate neutron irradiation.
They irradiated the samples with neutrons at 548 K and protons at 633 K to offset the difference in dose rate. As a result of their investigations, it was shown that radiationinduced segregation (RIS) behavior at the grain boundary and the defect kinetics under neutron irradiation can be replicated using proton irradiation (Was et al., 2002). Nonetheless, large discrepancies between ion and neutron irradiation are also reported in other situations (Was, 2017). For the second approach, a divide and conquer strategy is often employed to decouple the distinct irradiation-induced phenomena that usually span multiple length and time scales. Moreover, different atomistic, mesoscale, and macroscale models are often utilized to relate the evolution of the atomistic-and micro-structure of irradiated materials to the resultant change of their properties (Was, 2017). Nevertheless, the decoupled models proved to be still difficult to solve because of their multiphysics nature.
The computational approach is, however, very promising given the ever-increasing computational power. Moreover, it can be used to either guide the ion irradiation experiments to better simulate neutron irradiation or extrapolate their results to neutron irradiation (Xu et al., 2012). Classically, several approximations were utilized to lower the computational cost associated with numerical simulations. For instance, in ratetheory models that simulate long term irradiation effects, a homogenization scheme is often utilized, which ignores the underlying microstructure and spatial resolution of irradiation defects (Was, 2017). Clearly, this limits the model to obtain predictions only about the averages of quantities of interest. However, other models that relax these assumptions were introduced based on the so-called spatially resolved rate-theory (SRRT), which requires solving coupled partial differential equations representing different types of defects (Thompson, 1974;Bennett, 1975;Franklin, 1975;Warburton and Turnbull, 1975;Ullmaier and Schilling, 1980). Several recent investigations have implemented these models to gain more insight on irradiation effects on materials (Short et al., 2016) solved rate theory equation by accounting for injected interstitials and investigated the changes in the spatial void swelling profiles . Studied proton irradiation induced dislocation loops in modified 310S steels. They used a series of rate equations for different defect types including self-interstitial atoms (SIAs), vacancies, interstitial clusters, vacancy clusters, dislocation lines, and grain boundaries (Choi et al., 2020) used rate theory to understand irradiation creep behavior of ironbased alloys exposed to stress-applied environment. They coupled rate theory equations with a model for creep and analyzed the behavior at different dose rates (Saidi et al., 2021) presented an approach to study the effect of dose rate on radiation damage accumulation in polycrystalline systems in the presence of internal and grain boundary sink strengths. They used both proton and electron irradiation in their study and calculated the relation between internal sink strength to the grain boundary sink strength. All these studies demonstrated the importance of resolving the spatial dependences of both the underlying microstructure and irradiation defects for developing better understanding of the behavior of materials under irradiation.
Here, we follow the SRRT approach to investigate the surface and size effects on the transient and steady-state kinetics of point defects in irradiated solids. Moreover, we take into consideration the effects of dose rate, production bias, and defects recombination, reactions with distributed/volumetric sinks, and diffusion to localized/surface sinks. While there are a few attempts in literature to study the surface and size effects (Yang et al., 2010;Demkowicz et al., 2011;Xu et al., 2012), they are limited to small sizes and ignored one or more of the main irradiation conditions or defect mechanisms accounted for here. Indeed, this work demonstrates the importance of considering the combined effects of surface, size, dose rate, production bias, and temperature on the response of materials to irradiation. Specifically, we show that a few transient and steady-state irradiation-induced instabilities could be observed only if all these factors were accounted for simultaneously in the simulations. Furthermore, these instabilities lead to the development of different patterns in the steady-state concentrations of point defects that can either increase or decrease their ability to cluster into extended defects, and hence determine the extent of irradiation damage in materials.
This paper is organized as follows. In section 2, the SRRT approach is outlined and its numerical implementation is detailed. In section 3, the results are presented and discussed. Lastly, in section 4, concluding remarks are drawn and recommendations for future studies are given.

Balance Equations for Point Defects
The concentration of single vacancies and interstitials (e.g., point defects), C i and C v , are described mathematically by the chemical rate balance equations given in following sets of equations, There are different terms with different rate coefficients as seen in Eq. 1. The first term on the RHS is a source term, representing defect generation with a rate constant of K 0 . Defect production can be either uniform or non-uniform depending on irradiation type. The second term is a reaction term, referring to the recombination process between interstitials and vacancies with a rate constant of K iv . The third term is also a reaction term, representing the reaction between the defect x (interstitial or vacancy) and the sink s with a rate constant of K xs .
where r iv and r xs are interaction radii for the regarding reaction.
Interaction radii are usually taken as a multiple of the lattice constant, a, e.g., 10a. In this study, the sink concentration, C s , is assumed constant and uniformly distributed over domain. The last term represents the diffusion of defects (Motta et al., 2017;Was, 2017). Lastly, b is a production bias factor to be discussed next.
In principle, atomic collisions create vacancy-interstitial pairs (Frenkel Pairs), and they might recombine to restore the equilibrium state. Nonetheless, irradiation cascades usually have an underlying structure with vacancy rich core and interstitial rich periphery. This leads to the formation of small vacancy and interstitial clusters. Interstitial clusters are usually more thermally stable than their vacancy counterparts, which usually dissolve and inject vacancies back to the matrix. Eventually, this results in the presence of more single vacancies surviving the cascade stage than single interstitials. Moreover, the presence of dislocations in materials before irradiation (because of processing) may also increase this imbalance in point defects because of the preferential absorption of dislocations to interstitials. For models that are concerned with the long term irradiation effects, as in the case here, this defect imbalance is usually captured using a bias factor in the production rate (Motta et al., 2017;Was, 2017). We follow the same procedure here by introducing the bias factor, b, in the production rate of vacancies.
Another common form for the point defect equations is usually written by using sink strength, k 2 , definition for the defect-sink reaction term as follows, Subsituting Eq. 3 into Eq. 1 gives This definition of sink strength is beneficial in understanding the point defect-sink reactions. It is also important in evaluating the differences between the homogenized rate-theory and the spatially-resolved rate theory. We will utilize it here also to assess the difference in the predictions between the two theories for different microstructures and irradiation conditions.

Non-Dimensionalization of Balance Equations for Point Defects
To simplify the implementation of the equations and improve the convergence of solution, non-dimensionalization of the above equations was performed. Specifically, the following nondimensionalization scheme and scaling factors are introduced.
where Ω is the atomic volume in m 3 , l is a characteristic length scale in m, and ω is a characteristic time scale in s, ω l 2 /D i . After substituting Eq. 5 into 1 and reorganizing, the final form of the non-dimensionalized dynamical system becomes, defining new non-dimensionalized parameters, the final form is, The non-dimensionalized form of Eq. 4 can be derived easily from Eq. 8 by using Eq. 3.
where k 2 x is non-dimensionalized sink strength for defect x, k 2 x l 2 k 2 x .

Model Parameters and Numerical Implementation
The Multiphysics Object-Oriented Simulation Environment (MOOSE) is used as a tool to solve the balance equations for point defects numerically. MOOSE is an open-source framework that provides users with several capabilities to solve partial differential equations (PDEs). It contains a multitude of built-in features such as different linear and non-linear solvers, different integration schemes, mesh and time adaptivity, mesh generation and post-processing options, and many more (Permann et al., 2020) MOOSE has been utilized before to investigate irradiation effects and coevolution of microstructure and properties in materials (Ahmed and El-Azab, 2018;Ozturk et al., 2018;Tonks et al., 2018;Abdoelatef et al., 2019;Badry et al., 2019;Ozturk, 2019;Ahmed and El-Azab, 2020;Badry, 2021). We conducted 1D simulations in this study. The kinetic equations were solved in both Cartesian and spherical coordinates to represent planar and spherical grains. The trends were similar, and the differences were minor. Hence, for most of the results presented here, Cartesian coordinates were used unless stated otherwise. The domain boundaries are assumed to represent free surfaces or grain boundaries. Hence, they are modeled as perfect and neutral, i.e., point defects take on their equilibrium values at the boundaries. Sinks such as voids and dislocations, which are present in the bulk, are assumed uniformly distributed with a concentration of 1e18 m− 3 . 1D first-order Lagrange elements were used for spatial discretization and a second-order implicit scheme was utilized for time integration. A grid size in the range of 0.1-1 nm was used in the simulations. Mesh sensitivity studies were conducted, and the results reported here are confirmed to be independent of grid size within the range mentioned. A length scale, l, that is much smaller than the lowest grain size was used in our simulations. Since the smallest grain size was 5 nm, a length scale of 0.1 nm was used in all simulations.
Ni-based alloys are used in many nuclear reactors because of their excellent corrosion resistance and high mechanical properties. For example, Ni-based alloys containing 15-30% Cr are generally used in the cooling systems of reactors due to their corrosion resistance. Ni-Cr based alloys are also candidates for structural materials in molten-salt reactors. Therefore, here we selected Nickel as a base material for this study. Ni properties used in the simulations: lattice parameter, a, is 0.352 nm, atomic volume, Ω, is 0.01206 nm 3 , migration energy for interstitial, E m,i , is 0.30 eV, migration energy for vacancy, E m,v , is 1.30 eV, formation energy for interstitial, E f,i , is 4.27 eV, formation energy for vacancy, E f,v , is 1.60 eV, pre-exponential diffusion coefficient for interstitial, D 0,i , is 1e-7 m 2 /s, and preexponential diffusion coefficient for vacancy, D 0,v , is 6e-5 m 2 /s (Walgraef et al., 1996;Zhao et al., 2016). In most of the simulations presented here, a temperature of 773 K was assumed unless otherwise explicitly stated. The temperature dependence of diffusion coefficients follows the well-known Arrhenius relationship. Reaction rate constants, K iv , K is and K vs are calculated using Eq. 2. The interaction radii, r iv , r is and r vs are 3.52 nm. The calculated reaction rate constants, K iv , K is and K vs are 4.8958e-17 m 3 /s, 4.8498e-17 m 3 /s, and 3.7795e-21 m 3 /s.

RESULTS
The defect production term, K 0 , depends on the type of the incident radiation. Its magnitude and spatial dependence are determined by the interactions of incoming particles with the target material. Neutrons have an energy spectrum that varies by a couple of orders of magnitude. Defect generation is possible at different energies for neutrons. Also, because of their charge neutrality, they are able to travel through materials to farther distances corresponding with their energies. Thus, they cause a uniform or flat-like defect production profile in target materials. Electrons and protons are not able to penetrate as much as neutrons do, but for small domain sizes (<5 μm), their defect production can be assumed uniform as well. On the other hand, ions have a relatively narrow energy channel width, which is almost monoenergetic. They also cannot penetrate through material excessively because of their high electronic energy loss in collisions. This characteristic behavior of ions causes a non-uniform defect production profile. We have conducted simulations for both uniform and non-uniform irradiation. While there are several quantitative differences, the main qualitative trends of interest are similar. Therefore, here, for Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 684862 the sake of brevity and clarity, we will discuss only the results of uniform irradiation. Distinct irradiation types lead to different production rates of point defects. Nonetheless, for investigations concerned with long term irradiation effects, a range of dose rates representing different irradiation types is usually employed. For instance, a range of 1e-4 to 1e-2 dpa/s for energetic electrons/protons/ions and a range of 1e-8 to 1e-6 dpa/s for neutrons are usually quoted. The defect production bias, b, also strongly depends on the type of irradiation. This is because of the different cascade structures and densities produced by different types. Neutrons and heavy ions result in highly dense cascades, leading to high bias (because of the higher stability of small interstitial clusters relative to their vacancy counterparts). Electrons, on the other hand, result in isolated point defects, with very low bias. Protons and light ions lie between those two limits (Was, 2017).
In most rate-theory based studies, most of the focus is usually on the effects of temperature, dose rate, and sink densities. The effects of size and isolated/discrete surfaces are often ignored. Only a few studies (Yang et al., 2010;Demkowicz et al., 2011) investigated to some extent these effects. However, those studies were limited to small size grains (e.g., nanograins) and ignored the effect of production bias. As will be evident from our study here, surface and size effects play a major role in determining both the transient and steady state profiles of point defects in irradiated solids. Moreover, it will be demonstrated that there is a critical size at which a nonlinear instability is triggered, which in turn results in a specific pattern that favors the separation of vacancies and interstitials and hence promotes irradiation damage. That critical size is influenced by does rate, production bias, and temperature. Therefore, considering the influences of temperature, size, dose rate, production bias, and their combined effect is required in any future studies of irradiation effects and ignoring any of these principal parameters is unjustifiable.

Surface and Size Effects on the Transient and Steady-State Behaviors of Point Defects in Irradiated Ni
To validate the numerical implementation of the balance equations of point defects in MOOSE, the solution to a simple well-known benchmark problem is reproduced first. In this problem, analytical solutions to the point defect concentrations in a spherical grain with a perfect and neutral boundary/surface are derived. To simplify the derivation, the production bias and mutual recombination are ignored. For such conditions, the steady state analytical solutions for a spherical grain of radius a are found as below (Heald and Harbottle, 1977).
As evident from Figure 1, for a spherical grain with a radius of 500 nm, MOOSE simulation results (dots) agree very well with the analytical solutions (lines).

Concentrations
The defect concentrations can follow different patterns depending on their properties. Interstitials diffuse to sinks faster than vacancies, which lowers their steady-state concentration in the bulk. Figure 2 shows the evolution of the average fractions of point defects until steady state is reached. In these simulations, a production bias of 5% is assumed and the effects of grain size and dose rate are explored. As expected, longer times are required to reach steady state for larger grains because of the greater diffusion distances. In most cases, also as expected, higher average fractions of defects are associated with larger grains. This is due to the higher surface to volume ratio for smaller grains, which is usually quoted as the origin of radiation tolerance in nanocrystalline materials (Was, 2017;Yang et al., 2010;Demkowicz et al., 2011). However, there are a few exceptions to this latter rule. The transient vacancy concentrations in some of the smaller grains are actually higher than their large counterparts. For instance, at the low dose rate of 1e-6 dpa/s (Figure 2A), the transient vacancy concentration at 500 nm is higher than at 5,000 nm. Moreover, for the high dose rate of 1e-3 dpa/s ( Figure 2B), the transient vacancy concentration at 50 nm is higher than at both 500 and 5,000 nm. This anomaly in the dependence of vacancy accumulation on grain size was observed before in the study by (Yang et al., 2010), where lattice kinetic Monte Carlo simulations were employed. This was attributed to the competition between grain boundary absorption of defects and bulk recombination. Both processes have different length and time scales. Nonetheless, the steady state concentrations satisfy the general rule that smaller grains tend to accumulate fewer defects, which also agrees with (Yang et al., 2010). The other exception that we report here is a different anomaly related to the accumulation of interstitials. Above a critical grain size, both the transient and steady-state concentrations of interstitials reduce with increasing grain size. For a dose rate of 1e-6 dpa/s, the critical size is about 1,000 nm ( Figure 2A). The critical size is about 250 nm for a dose rate of 1e-3 dpa/s. It is worth noting that for the high dose rate the average interstitial fraction at 5,000 nm is almost the same as its counterpart at 5 nm, while the corresponding vacancy concentration is about four orders of magnitude higher than its counterpart at 5 nm. This anomaly in the interstitial behavior was not detected in (Yang et al., 2010;Demkowicz et al., 2011) because of the limitation to small sizes (<50 nm) in those studies.
The effect of grain size on the point defect kinetics resembles to some extent the effect of temperature. If compared to the classical results from homogenized rate-theory (Was, 2017), which ignores surface, diffusion and size effects, as the grain size changes from small to large, the results show similar trends as if the temperature changes from high to low. This is due to the fact that both factors affect the diffusion kinetics of defects.
Nonetheless, the effect of size is even more complicated because of the anomalies discussed above.
The anomaly in the transient behavior of vacancies with grain size is similar to the cases reported in (Yang et al., 2010). However, the nature and origin of the anomaly in the behavior of interstitials are different. This is evident from the fact that this new anomaly persists even under steady-state conditions. Hence, a close examination of this anomaly is in order.
The anomaly in the dependence of interstitial accumulation on grain size is associated with an instability that once triggered, different spatial profiles of vacancies and interstitials are established. This is captured in Figure 3 that depicts the evolution of the interstitial and vacancy profiles with time for a 500 nm planar grain with 1e-3 dpa/s dose rate and 5% production bias. Initially, the concentrations of both defects increase in the center and reduce at the surfaces, which act as sinks. However, as time progresses, vacancy and interstitial concentrations start to show opposite trends, where vacancies tend to accumulate in the center and interstitials tend to accumulate by the surfaces. This instability appears effectively in relatively large grains (>100 nm) with any nonzero value of production bias (even as low as 0.001%). The initiation of this instability is attributed to the fact that production bias breaks the symmetry of the balance of point defects leading to non-uniform bulk recombination and non-uniform losses to surfaces (effectively transforming neutral sinks to biased), which eventually results in the development of distinct steady-state profiles for the different point defects.  (Yang et al., 2010). Second, transient and steady-state interstitial concentrations decrease with grain size above a critical size (e.g., 1,000 nm for the low dose rate and 250 nm for the high dose rate). For the high dose rate, the average steady-state interstitial concentration at 5,000 nm is almost equal to its counterpart at 5 nm, while its corresponding vacancy concentration is, however, four orders of magnitude higher than its counterpart at 5 nm.
Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 684862 The opposite trends of vacancies and interstitials arise because of their distinct dynamics. Interstitials are produced at a lower rate than vacancies. Also, interstitials are faster than vacancies. If one runs a hypothetical simulation with those conditions reversed, the steady-state spatial profiles of vacancies and interstitials will also be reversed.
Since the diffusion coefficients are temperature dependent, the D i /D v ratio changes with temperature. To investigate the effect of diffusivity ratio/temperature, several simulations were conducted for different temperatures with a production bias of 5% and a production rate of 1e-3 dpa/s in a 1,000 nm grain. The results are presented in Figure 4. As evident from the figure, there is a critical ratio (≈33 for Ni) below which the instability is suppressed. In other words, there is a critical temperature only below which the instability develops.
This instability is also sensitive to grain size, dose rate, and production bias. This is demonstrated in Figures 5, 6. In Figure 5, the effect of grain size is shown for two different dose rates with a production bias of 5%. In all cases, the steady-state vacancy concentration increases with grain size and is highest in the center. The vacancy concentration gradients in the vicinity of surfaces increases with dose rate. On the other hand, the steadystate interstitial concentration initially follows the same trend as the vacancy concentration up to a critical size, above which the trend is reversed, i.e., interstitial concentration decreases with grain size. This critical size is 250 nm for a dose rate of 1e-3 dpa/s and 1,000 nm for a dose rate of 1e-6 dpa/s. This is consistent with Figure 2, and hence this instability is indeed the origin of the size anomaly discussed above. The production bias also plays an important rule. Using 3D plots, Figure 6 demonstrates the combined effects of size, dose rate, and production bias.
Generally speaking, the critical size reduces with dose rate and production bias.
It is worth mentioning that such instability will render irradiation damage more pronounced. This is because the pattern that develops based on this instability favors the separation of vacancies and interstitials. Therefore, it will enhance the possibility of clustering of different point defects into distinct extended defects. Namely, vacancies at the center of the domain will cluster into voids or vacancy loops, while interstitials close to the boundaries will form dislocation loops. It is also interesting to note that this instability is suppressed in nano-grains (with sizes less than 100 nm). This provides a new insight into the stability of nano-materials, i.e., such materials support the development of point defect patterns that limit the formation of extended defects.
It is unfeasible to experimentally measure the densities of single point defects in a direct way. However, our results here can be indirectly validated by the well-known existence of void denuded zones close to grain boundaries. Since the instability in our results show accumulation of interstitials and depletion of vacancies adjacent to GBs, it is expected that void formation will be suppressed under such conditions.
Additionally, the dependence of the severity of radiation damage on grain size predicted here agrees with the recent experimental study by (Mao et al., 2019). In that study, the sink strength of nano-GB is quantitatively investigated in the nano-grained Cu and dilute Cu-W alloys with the average grain size ranging from 20 to 500 nm by 300 keV He ions irradiation at room temperature and 673 K. The irradiation induced void volume ratio, size and distribution are confirmed to strongly depend on the grain size and irradiation temperature. It was found that there is a critical size below which the samples are completely radiation resistant and almost no voids form. Moreover, such full radiation tolerance was only achieved at the high temperature but not the low temperature. Those findings are clearly in agreement with our prediction here. Lastly, it is noteworthy that this instability is different from the Turing instability usually invoked to explain void/bubble superlattices formation in irradiated solids (Krishnan, 1980;Ghoniem et al., 2001;Gao et al., 2018;Noble et al., 2020). Turing instability can be inferred from linear stability analysis and depends on the ratio and/or anisotropy of diffusion coefficients. The current instability, however, is nonlinear, surface induced, and size regulated.

Vacancy Supersaturation
Void formation is one of the most studied processes in irradiated solids (Krishnan, 1980;Ghoniem et al., 2001;Was, 2017;Gao et al., 2018;Noble et al., 2020). This process often takes place whenever a vacancy supersaturation develops in irradiated solids, i.e., when the vacancy concentration is much higher than the equilibrium concentration. Vacancies tend to cluster together into voids to lower the overall free energy of the system. Void formation usually leads to void swelling. The presence of voids and the dimensional change associated with void swelling drastically alter the thermo-mechanical properties of materials.
In this section, we investigate the effects of surface, size, dose rate, and production bias on the development of vacancy supersaturation in irradiated solids. An effective vacancy supersaturation, S v , can be estimated from the following equation (Was, 2017), where C e v is the equilibrium vacancy concentration. The effective vacancy supersaturation increases with vacancy concentration and reduces with interstitial concentration. Figure 7 presents the effects of dose rate, production bias, and size on vacancy supersaturation. First, the spatial dependence of vacancy supersaturation as function of size for two different dose rates is captured in the 2D plots. Following the vacancy concentration (recall Figure 3), the maximum vacancy supersaturation is largest in the center. In general, vacancy concentration increases with dose rate, production bias, and size. For the same size, temperature, and production bias, the dependence on dose rate is almost linear. However, the dependence on size or production bias is non-linear. This is clear from the 3D plots that shown the combined effects of size and bias. The effect of size is detailed in Figure 8 that depicts the change of the maximum vacancy supersaturation with size at a dose rate of 1e-6 dpa/s and 5% production bias. As obvious from the figure, the development of vacancy supersaturation is suppressed in submicron grains, while facilitated in large grains. Not only vacancy concentration increases with grain size, but its rate of increase also climbs with size. It is noteworthy that the increase of vacancy concentration becomes more apparent above the critical size. For instance, the ratio between the maximum vacancy concentration at 500 nm compared to 50 nm is 67.6, while the corresponding ratio at 5,000 nm compared to 500 nm is 89.97. This supports our earlier prediction that the bias-enabled, surface-induced, and size regulated instability and developed pattern exaggerate irradiation damage.

Grain Boundary Sink Strength
The sink reaction term in Eq. 4 represents the total loss of point defects to sinks in the bulk such as voids and dislocations. It can be written in expanded form as below, FIGURE 5 | Effect of grain size on the steady-state concentration profiles of point defects at 5% production bias. Upper row: interstitials at (A) low dose rate (1e-6 dpa/s), (B) high dose rate (1e-3 dpa/s). Lower row: vacancies at (C) low dose rate (1e-6 dpa/s), (D) high dose rate (1e-3 dpa/s).
Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 684862 9 Similarly, sink strength of a grain boundary can be calculated by considering total flow of point defects to that grain boundary.
where F is total flow at grain boundary, Z gb is sink strength of grain boundary, D is diffusion coefficient of the point defect, and C 0 is concentration of point defect at the domain center.
The total flow at boundary is equal to multiplication of surface area and current at boundary as below, The total grain boundary sink strength, k 2 gb is calculated by multiplying the sink strength of a grain boundary, Z gb , by the grain density, ρ gb , (Heald and Harbottle, 1977).
where d is diameter of the grain, d 2a.
For the benchmark problem of a spherical grain described earlier, grain boundary sink strengths, Z gb , and k 2 gb , are calculated analytically by substituting Eq. 11 into the expressions above. First, the concentration at the grain center and total flow at the FIGURE 6 | An illustration of the combined effects of size, dose rate, and production bias on the steady-state profiles of interstitials. Change of interstitial concentration profiles with production bias for (A) low dose rate (1e-6 dpa/s) at 50 nm grain, (B) high dose rate (1e-3 dpa/s) at 50 nm grain, (C) low dose rate (1e-6 dpa/s) at the critical grain size (≈1,000 nm), (D) high dose rate (1e-3 dpa/s) at the critical grain size (≈250 nm).
Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 684862 grain boundary are expressed in terms of the production rate and the bulk sink strength, k 2 , as follows Heald and Harbottle (1977).
Then, grain boundary sink strengths are given by Eq. 18 predicts the same individual (for a grain boundary) and total sink strengths for both point defects in the presence of neutral bulk sinks. However, dislocations are known to be biased towards the absorption of more interstitials than vacancy because FIGURE 7 | The effects of dose rate, production bias, and size on vacancy supersaturation. Upper row: Vacancy supersaturation profiles for 5% production bias for (A) low dose rate (1e-6 dpa/s), and (B) high dose rate (1e-3 dpa/s). Lower row: Maximum vacancy supersaturation for different values of production bias and grain size for (C) low dose rate (1e-6 dpa/s), and (D) high dose rate (1e-3 dpa/s).
Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 684862 11 of their stress field. In that case, grain boundary sink strengths for different defect types can be calculated by using appropriate diffusion coefficient, D, bulk sink strength, k 2 , and concentration at center, C 0 . The analytical expressions are then given by the following equations, The numerical solutions from MOOSE for the benchmark problem were utilized to calculate the sink strengths of point defects for a spherical grain with a 500 nm radius. For simplicity, only neutral bulk sinks were considered. The differences between MOOSE results and the analytical expressions given in Eq. 18 were less than 0.02%. This is to be expected given the excellent agreement reported before for the concentration profiles (recall Figure 1).
We then investigate the effects of dose rate, production bias, and size on the boundary sink strength. Again, for the sake of simplicity, distributed sinks in the bulk are assumed neutral. Nevertheless, the individual and total sink strengths for vacancies and interstitials are different for non-zero values of production bias. This is demonstrated in Figures  9-12. Figure 9 shows the effect of grain size on both the individual and total boundary sink strength. The individual FIGURE 8 | The effect of size on maximum vacancy supersaturation at a dose rate of 1e-6 dpa/s and 5% production bias. (A) whole size range, (B) submicron sizes (C) micro sizes. The maximum vacancy supersaturation increases with grain size and the rate of increase also rises with grain size. The effect of grain size is more pronounced above the critical size (1,000 nm). The fitting equations are valid for sizes between 5 and 5,000 nm.
Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 684862 sink strengths for both vacancies and interstitials increase with size in agreement with the predictions of Eq. 18. However, their values start to deviate from each other above the critical size because of the instability and the resultant patterns mentioned before. The accumulation of interstitials close to the boundaries results in larger flow of interstitials to the surfaces, effectively transforming them into biased sinks. For the total boundary sink strengths for defects, they initially decrease with size as predicted form Eq. 18. Nonetheless, above the critical size, it increases with size for interstitials, while it continues to decrease with size for vacancies. This is consistent with the anomaly of the dependence of the steady-state concentration of interstitials on size discussed earlier. The dose rate influences the critical size, and hence the critical point where the boundary sink strengths for interstitials and vacancies become different. The dependence of boundary sink strengths on production bias are presented in Figures 10, 11. Figure 10 depicts the effect of production bias on the individual boundary sink FIGURE 9 | Change of grain boundary sink strengths with grain size for 5% production bias. Upper row: Individual grain boundary sink strengths, Z gb,i and Z gb,v , for (A) low dose rate (1e-6 dpa/s), and (B) high dose rate (1e-3 dpa/s). Lower row: Total grain boundary sink strengths, k 2 gb,i and k 2 gb,v , for (C) low dose rate (1e-6 dpa/s), and (D) high dose rate (1e-3 dpa/s). Because of the size-and dose rate-dependent instability, the individual and total boundary sink strengths for interstitials are higher than vacancies above the critical size. This indicates that the boundary changes from neutral to biased at the critical size.
Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 684862 strength. Below the critical size, the boundary sink strengths for both types of defects are the same and show no dependence on production bias. However, once the critical size is reached, the boundary sink strength for interstitials increases with bias, while it decreases with bias for vacancies. The same trends are manifested in the behavior of the total boundary sink strengths shown in Figure 11. The combined effects of size and production bias on boundary sink strengths for both defects are captured in Figure 12 using 3D plots. Again, the differences in the behaviors of the defects appear only above the critical size, which is dependent on dose rate and production bias (e.g., irradiation type). Therefore, in contrast to the predictions of Eq. 18, the boundary sink strengths not only depends on the material and its underlying FIGURE 10 | Change of individual grain boundary sink strengths with production bias and grain size. Upper row: Interstitials, Z gb,i , for (A) low dose rate (1e-6 dpa/s), and (B) high dose rate (1e-3 dpa/s). Lower row: Vacancies, Z gb,v , for (C) low dose rate (1e-6 dpa/s), and (D) high dose rate (1e-3 dpa/s). Below the critical size, the boundary sink strengths for both defects are equal and constant (independent of the value of production bias). However, at and above the critical size, the boundary sink strength for interstitials increases with bias, while it decreases with bias for vacancies.
microstructure, but it also depends on the irradiation type and temperature.

Size Effect on Radiation-Induced Segregation in Ni-Cr
The size-dependent instability reported here is expected to render irradiation damage more severe as discussed in the last section. It was demonstrated that the vacancy supersaturation rises sharply above the critical size (recall Figure 8). This will in turn facilitate the clustering of vacancies into voids and/or vacancy loops. The clustering process, though, is highly influenced by material properties such as binding energies of point defects, surface energy, and stacking fault energies. However, a more direct effect for the increase of point defect concentrations and gradients is the radiation-induced segregation (RIS) of solute atoms to surfaces/boundaries (Was, 2017). Depending FIGURE 11 | Change of total grain boundary sink strengths with production bias and grain size. Upper row: Interstitials, k 2 gb,i , for (A) low dose rate (1e-6 dpa/s), and (B) high dose rate (1e-3 dpa/s). Lower row: Vacancies, k 2 gb,i , for (C) low dose rate (1e-6 dpa/s),and (D) high dose rate (1e-3 dpa/s). The trends for the dependence of the total boundary sinks on production bias are similar to the trends for the individual boundary sink strengths discussed in Figure 10.
Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 684862 15 on whether solute atoms exchange preferably with vacancies or interstitials, they will be depleted or enriched at the surface/ boundary, respectively.
Here, we investigate the size effect on the RIS in a Ni-Cr binary alloy. Ni-Cr based alloys are one of the most important materials for current and future advanced nuclear reactors. These alloys show high strength, high resistance to swelling at high temperatures, and superior corrosion performance. Nonetheless, despite their many advantages, radiation induced segregation in these alloys is known to significantly alter their properties (Allen et al., 1997;Was et al., 2002;Allen et al., 2007;Wharry and Was, 2013). Therefore, studying the process of RIS in Ni-Cr based alloys is crucial for advancing the utilization of these alloys and improving the safety and performance of reactors. To investigate radiation-induced segregation in alloys, the inverse Kirkendall model is usually employed (Allen and Was, 1998;Wharry and Was, 2014;Was, 2017). In that model, in addition to the balance equations of point defects (Eq. 1) solved above, the continuity equations for the atomic species are coupled and solved concurrently. Specifically, for a binary Ni-Cr system, the balance equations are (Barr et al., 2015), The time derivatives in Eq. 20 represent the evolution of the concentration of interstitials, C i , the concentration of vacancies, C v , and the concentration of Cr atoms, C Cr . On right hand side, J i , J v , and J Cr are fluxes of interstitials, vacancies and Cr atoms. The second term in Eq. 20 is a source term representing defect generation with rate a constant of K 0 . Third term, K iv C i C v , is a reaction term, referring to the recombination between interstitials and vacancies. The last term is also a reaction term, representing the reaction between the defect x (interstitial or vacancy) and the sink s with a rate constant of K xs . The diffusivity of element k (Cr or Ni) for defect x, d x,k , is defined by the following equation: where d 0 x,k is pre-exponential factor of Ni or Cr for vacancy or interstitial, E m x,k is the migration energy of Ni or Cr for the vacancy or interstitial, T is temperature and k is the Boltzmann constant.
We expressed the partial diffusion coefficient by D dN, where N is atomic fraction and d is diffusivity (Eq. 21). The spatial dependence resides in the N. In contrast, the d is composition independent and contains the kinetic and diffusion information (Wiedersich et al., 1979).
where d v,Cr and d i,Cr are vacancy and interstitial diffusivities for Cr, d v,Ni and d i,Ni are vacancy and interstitial diffusivities for Ni. Therefore, fluxes take on the following expressions, where J v and J i are the fluxes of vacancies and interstitials, J Cr is flux of the Cr atoms. Ω is atomic volume.
In this study, a Ni-5Cr binary alloy system was used as model alloy. Ni and Cr properties used in the radiation induced segregation model: lattice parameter is 0.352 nm, pre-exponential factor for Ni interstitial diffusivity 5.04e-8 m 2 /s, pre-exponential factor for Cr interstitial diffusivity 3.20e-7 m 2 /s, pre-exponential factor for Ni vacancy diffusivity 1.85e-4 m 2 /s, pre-exponential factor for Cr vacancy diffusivity 2.26e-4 m 2 /s, activation energy for Ni interstitial diffusivity 0.30 eV, activation energy for Cr interstitial diffusivity 0.37 eV, activation energy for Ni vacancy diffusivity 1.16 eV, activation energy for Cr vacancy diffusivity 1.10 eV, vacancy formation energy 1.79 eV, interstitial formation energy 4.0 eV (Barr et al., 2015). For all simulations in this section, the Cartesian coordinate system was used to represent planar geometry.
We simulated radiation induced segregation model at an irradiation temperature of 773 K, the most common temperature used for RIS studies (Allen et al., 2007;Wharry and Was, 2013;Barnard and Morgan, 2014;Barr et al., 2015), and dose rates of 7e-6 and 1e-3 dpa/s. In one of the simulations, temperature was varied between 473 and 1073 K to study its effect. Vacancy and interstitial concentrations were determined as the concentrations at thermal equilibrium for the initial conditions of system. Ni and Cr concentrations are taken as the nominal value in the alloy composition. We assumed that the alloy atoms are randomly distributed throughout the grain at the initial time and bulk sinks are uniformly distributed with a concentration of 1e18 m− 3 . We conducted several simulations to investigate the effects of size, production bias, dose rate, and temperature on RIS in Ni-5Cr. In all of these simulations, the boundary is assumed neutral and perfect, and zero flux boundary condition was used for Cr. Due to the symmetry of the problem, only one boundary (at x 0 in the figures) is considered while at the other far end (representing the grain center), zero flux boundary conditions are considered for point defects and Cr. The equations were solved numerically with the initial and boundary conditions, and 1D simulations were performed.
Our model predictions for the cases with large grain sizes and without production bias agree well with the numerical results and experimental data presented in (Barr et al., 2015). However, we demonstrate here a strong dependence of RIS on size and production bias, which was not explored in (Barr et al., 2015). The quantitative effects of size and production bias on RIS in Ni-5Cr are presented in Figures 13, 14. In Figure 13, the effect of production bias is shown for two different grain sizes at a temperature of 773 K and a dose rate of 1e-3 dpa/s. As evident from the figure, the depletion of Cr at the boundary becomes more pronounced as the production bias increases. Both the boundary Cr concentration is reduced and the depletion layer is widened with production bias. This is consistent with the fact that Cr exchanges preferably with vacancies, since the concentration of vacancies rises with production bias. However, it is noteworthy that the effect of production bias is more apparent in larger grains. For instance, in the case of 500 nm size with 10% production bias, the steady-state vacancy concentration is one order of magnitude higher than its corresponding value in the absence of production bias. Conversely, the steady-state interstitial concentration is one order of magnitude lower than its corresponding value in the absence of bias. This is because of the size-dependent instability discussed before, i.e., the instability is only observed in large grains. In Figure 14, the effect of size on RIS in Ni-5Cr is explicitly investigated for two different dose rates in the presence and absence of production bias. The size effect saturates above 250 nm in the absence of production bias. Nonetheless, in the presence of production bias, the depletion of Cr increases with size. Moreover, the rate of increase rises with size. It is also clear that the magnitude of depletion and its size dependence are amplified at higher dose rates. At the low dose rate of 7e-6 dpa/s, the model predictions lie within the results obtained before (Barr et al., 2015), where Cr concentration at the boundary reported to be between 2 and 4%. Nonetheless, the predictions here for 20% bias show inverse dependence on grain size, i.e., Cr concentration at the boundary becomes smaller with increasing size. Moreover, at the high dose rate of 1e-3 dpa/s, the Cr boundary concentration becomes negligibly small at sizes above 250 nm. However, in this case, the treatment of the boundary as a perfect sink is probably no longer valid.
Lastly, the effect of temperature on RIS in Ni-5Cr was investigated. Figure 15 summarizes the results of simulations conducted for a 500 nm size irradiated at 1e-3 dpa/s with two different production biases. First, the dependence of point defects on temperature is obvious, i.e., as temperature increases, the steady-state concentration of vacancies decreases, while the interstitial counterpart rises. Therefore, the difference between the point defect concentrations is reduced at high temperature. Moreover, at the highest temperature of 1073 K, the instability is absent ( Figure 15A). However, it could still happen at this temperature for higher dose rate, higher production bias, or larger size. Hence, RIS is the least observed at the highest temperature. As the temperature decreases, RIS becomes more apparent with higher segregation of Ni and depletion of Cr at the boundary as the production bias increases. RIS is predicted here to be maximum at the lowest temperature of 473 K, while it is often reported to be maximum at intermediate temperature (Barnard and Morgan, 2014;Barr et al., 2015). This could be attributed to the uncertainties associated with interstitial diffusivities at lower temperature and the validity of the perfect sink assumption, as was discussed in (Barnard and Morgan, 2014). Nevertheless, the combined effect of size and production bias on RIS is clear.

CONCLUSION
In this work, we conducted a detailed study of the effects of surface and size on the transient and steady-state behavior of point defects in irradiated solids. Pure Ni and Ni-Cr were utilized as model systems. The balance equations for point defects were solved simultaneously using a fullycoupled and fully-implicit scheme implemented in the MOOSE framework.
It was found that in the presence of production bias, the effects of surface and size are pronounced even at low surface to volume ratio. This was attributed to a surface-induced and size-regulated instability. This instability appears above a critical size between 100 and 500 nm. The critical size decreases with increasing dose rate and/or production bias or lowering the temperature. Our detailed analysis revealed that this instability is triggered due to the fact that production bias breaks the symmetry of point defects flow to the surface/boundary, effectively transforming it into a biased sink. Below the critical temperature and at the threshold size, the imbalance in the flow of point defects is enough to create a pattern that favors the separation of vacancies and interstitials. Once this pattern develops, anomalies in the dependence on size for the steady-state concentrations of point defects and the surface/boundary sink strength are observed. Specifically, it is shown that once this instability is initiated, while the steady-state vacancy concentration continues to increase with size, its interstitial counterpart starts to decrease with size. On the other hand, while the surface/boundary sink strength for vacancies continue to decrease with size, its interstitials counterpart starts to increase with size. Hence, the developed instability and the resultant patterns will promote irradiation effects in materials. Indeed, this was demonstrated here by investigating vacancy supersaturation in Pure Ni and radiation-induced segregation in Ni-Cr. In the former, it was shown that vacancy supersaturation in Ni increases with size and the rate of increase also rises with size. In the latter, it was shown that the magnitude of enrichment/depletion of Ni/Cr at the boundary increases with size and the width of the enrichment/ depletion layer also increases with size.
It is worth noting that the above-mentioned instability and the resultant pattern and anomalies are absent in nanomaterials (with sizes <100 nm). This gives a new perspective on the resistance of these materials to irradiation as opposed to the susceptibility of conventional materials (with sizes >100 nm) to irradiation. In conventional materials, the developed pattern of point defects facilitates the clustering and formation of extended defects. Conversely, the exact opposite is observed in nanomaterials.
The results obtained here agree well with experimental observations in irradiated materials. For instance, the patterns appearing in the steady-state concentration profiles of point defects and the resultant pattern in the vacancy supersaturation agree well with the well-known observation of the development of void denuded zones in the vicinity of grain boundaries in irradiated polycrystalline solids (Was, 2017). Notably, our model can capture this behavior without the added assumption of 1D migration of interstitials, which was previously thought as a necessary condition to explain void swelling and its heterogeneity near the interfaces in materials. Moreover, the dependence of irradiation tolerance in materials on size and temperature demonstrated here is in agreement with the comprehensive experimental study conducted by (Mao et al., 2019) on irradiated Cu and Cu-W alloys.
It is also noteworthy to point that models/simulations that ignore defects diffusion to localized sinks (e.g., surfaces/ boundaries), bulk recombination of point defects, or production bias are incapable of capturing this size-dependent effects. Such reduced-ordered models or lower-fidelity simulations will not only have inaccurate quantitative predictions but may also have completely wrong qualitative predictions. Therefore, precautions must be taken in developing/conducting such models/simulations and interpreting their results.
It is also important to discuss the limitations of the current study. First, the treatment of the production bias here is simplified, i.e., its value was considered as unknown, and a parametric study of its effect was conducted. However, exact values for different irradiation types and conditions can be obtained through a multiscale modeling approach that simulate the development of irradiation cascades and the clustering of point defects (Was, 2017). Second, the evolution of extended defects and their heterogeneity were ignored. The coevolution of point and extended defects must be accounted for to develop a comprehensive understanding of irradiation effects in materials. Third, further microstructure descriptors in addition to size should be incorporated into the model since grain morphology, interfacial energy, and misorientation are expected to play important roles in the response of grain boundaries to irradiation (Was, 2017). Alleviating these limitations will be the focus of upcoming studies.
Lastly, we want to reflect on the implications of this study on the on-going debate on the utilization of ion irradiation to mimic neutron irradiation to accelerate the qualification of novel materials for advanced reactors. It is clear from all the results presented here that there is a strong and non-linear quantitative dependence, regardless of size, on production bias, dose rate, and temperature, i.e., irradiation type and conditions. Moreover, when surface and size effects are considered, qualitative differences in the irradiation response of materials to different irradiation types are expected. In fact, we proved that, in addition to its known dependence on material and size/microstructure, the surface/boundary sink strength is dependent on irradiation type. Hence, the generalization of predictions based on ion irradiation data to neutron irradiation is generally inaccurate. However, the utilization of physics-based models and high-fidelity simulations, as presented here, can be used to bridge the gap between the different irradiation types and improves the generalization. This can be accomplished by leveraging ion irradiation data to validate the models and then using the models to simulate the response under neutron irradiation.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: MOOSE is an open source code that can be freely downloaded from https:// mooseframework.inl.gov/. The input files and plots for this study are available for download from a public repository located at https://github.com/abdurrahmanozturk/SRRT. All data generated during this study can be reconstructed by running MOOSE input files.

AUTHOR CONTRIBUTIONS
AO, MG, and KA contributed to conception and design of the study. AO, MG performed the numerical simulations. AO, MG, and KA wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.