Stochastic Cluster Dynamics Simulations of Irradiated Materials

The study of irradiation damage and its accumulation and effects on material properties is an exceedingly complex problem that has attracted considerable atention over the last several decades. Its very nature involves a number of features that contribute to this complexity, such as:

Irradiation damage involves nonlinear interactions between defects. These require two- and many-body terms when formulating the kinetic evolution of the system.
Defect species move on radically different time scales depending on whether they are of vacancy or self-interstitial type. This makes it very difficult to study the long term evolution of the irradiated system if one has to resolve the fast dynamics constantly.
Irradiation brings about several convolutions in space and time. What happens in one region of the microstructure at a given time is related to what happens in a different region at a different time. The irradation history of the material matters considerably in a spatio-temporal sense.
Irradiation damage is a numbers game. It generally originates in the form of a displacement cascade, which produces just a handful of defects or defect clusters. Yet, the entire material's performance evolution may depend on the fate of these few defects accumulated over long times and large volumes. Miscounting or failing to account for small numbers of defects like these may lead to completely incorrect results.
Irradiation damage is the epitomy of multidisciplinarity in physics. Besides the general use of thermodynamic and kinetic principles, it combines elements of electronic structure, lattice dynamics, ballistics, alloy chemistry, mass transport, plasticity, dimensional changes, and mechanical behavior, among others, and it is inherently multiscale, as it spans many orders of magnitude in length and time scales.

Reaching timescales (or, equivalently, irradiation/implantation doses) of relevance for microstructural evolution requires approximations to reduce the complexity of the problem posed above. As well, often times, defect clusters can grow to very large sizes, contaning >10⁶ defects or so (particularly for high irradiation doses or long aging times). The most popular simplification used is the so-called mean field approximation, in which the interactions of a many-body system are replaced with an average or effective interaction. This reduces a multi-body problem into an effective one-body problem, resulting in significant computational gains at the expense of neglecting spatial correlations and stochastic fluctuations. In the context of defect population evolution, mean field rate theory (MFRT) provides solutions to a Master Equation (ME) that describes both growth and dissolution of the clusters due to reactions with mobile defects (or solutes), thermal emission of point defects from the clusters, and cluster coalescence if the clusters are mobile. Discretization of the ME in terms of the cluster species leads to a system of coupled partial differential equations (PDE) in which time is an explicit variable. The number of equations is the same as the number of cluster species (and/or solutes) in the largest possible cluster. Numerical integration of such a system is feasible on modern computers, although practicable solutions to the ME still make use of other numerical tricks to speed up the calculations. All of the above points, however, imply that the coeffcient matrices characterizing the PDE system will be very large, sparse in some regions of the cluster dimensional space, while concentrated in others, with elements that have widely disparate values, and which can be time-evolving. This represents a mathematical challenge that even the most formidable deterministic solvers can struggle with. Another limitation of MFRT models is that for multi-species systems (e.g., as it happens under fusion neutron irradiation conditions), the number of equations increases exponentially as the number of different species grows (a phenomenon known as 'combinatorial explosion'). This makes the application of these methods to systems with complex chemistries (such as when transmutation or alloy chemistry are important) highly impractical.

When conducting 'material-point' (bulk) simulations, one can reduce the problem yet a bit more to solving a set of ordinary differential equations (ODE) by dropping the spatially-dependent diffusion term from the equations. This eliminates the need for spatial discretizations to solve the spatial derivatives in the PDE system. In this case, the simulations truly become dimension-less and no notion of volume or finite length scale exists. This is generally the preferred calculation mode for irradiated structural materials.

Capturing spatial and time correlations and fluctuations requires solving every defect event explicitly. In reality, an 'event' has concrete time and space characteristics, and thus occurs in a spcific instant at a specific location. Other than atomistic methods, which are typically not capable of covering the range of spatio-temporal scales needed for damage accumulation, the most popular class of methods are 'event' and 'object' kinetic Monte carlo (kMC) methods, whether in their lattice or continuum variants. KMC simulations are inherently stochastic and hence introduce 'numerically' what otherwise occurs 'naturally', i.e., variations of defect processes in space and time. For this, one has to typically map the defect cluster kinetics to a Markov process, which itself makes use of several assumptions. Because kMC simulations evaluate events sequentially, one at a time, they are not subjected to exponential dimensional growth in cases of multiple species, which is a major advantage over standard deterministic methods. In the limit, kMC and MFRT simulations converge to the same solutions (in reality what converges is the average solution over a number of independent kMC simulations), but in most practical cases kMC simulations are limited in terms of the doses they can achieve. So the extra information that kMC methods can furnish is of course gained at the expense of computational efficiency.

There is, however, a way to take advantage of the beneficial aspects of both classes of techniques. One can solve the PDE/ODE system within the mean-field approximation but treating each term in the coefficient matrix as if it was an event rate. In this fashion, one naturally introduces the stochastic variability of explicit kMC methods but the mathematical form of MFRT methods is preserved. Event rates must defined for a specific material volume, much like for kMC simulations, but matrix coefficients are obtained as for MFRT simulations. Likewise, spatial derivatives (diffusion terms) can be converted to rates in the simulation volume using the divergence theorem and using standard discretizations as in the deterministic case. The result is a 'hybrid' method of sorts, which supresses combinatorial explosion, introduces spatio-temporal correlations naturally, and is capable of reaching irradiation doses much longer than pure kMC simulations. We call this method Stochastic Cluster Dynamics (SCD).

We are applying SCD to problems involving species such as He and H and how they interact with point defect clusters in tungsten-based plasma facing materials. Our research involves spatially-resolved simulations to mimic plasma exposure (He and H atoms) to tungsten surfaces. We have also been studying the effect of migration dimensionality of defect clusters on the kinetic coefficients of defect interactions that we use in the equations. We are also merging SCD with crystal plasticity models of deformation to simulate concurrent irradiation/deformation processes, such as hardening, swelling, and creep. We do this in a variety of materials such as iron, tungsten, or ceramics.