A simple model for the stochastic evolution of defects in a material under irradiation is presented. Using the master-equation formalism, we derive an expression for the average number of defects in terms of the power flux and the exposure time. The model reproduces the qualitative behavior of self-healing due to defect recombination, reaching a steady-state concentration of defects that depends on the power flux of the incident radiation and the material temperature, while also suggesting a particular time scale on which the incident energy is most efficient for producing defects, in good agreement with experimental results. Given this model, we discuss the integral damage factor, a descriptor that combines the power flux and the square of the irradiation time. In recent years, the scientific community involved in plasma-facing materials for nuclear fusion reactors has used this parameter to measure the equivalent material damage produced in experiments of various types with different types of radiation and wide ranges of power flux and irradiation time. The integral damage factor is useful in practice but lacks formal theoretical justification. In this simple model, we find that it is directly proportional to the maximum concentration of defects.

## I. INTRODUCTION

The study of materials exposed to irradiation is a broad topic whose applications are fundamental for medicine,^{1} astrophysics,^{2} and the design of nuclear reactors.^{3,4} The effects of radiation on materials are varied, at both the atomic and the microstructural scale, and depend mainly on the energy of the incident particles (neutrals, ions, electrons, and neutrons) impacting the surface. These effects include heating, embrittlement due to ion implantation,^{5} cracking and formation of voids,^{6} local melting and recrystallization or amorphization,^{7} evaporation and ablation of various types, and even different degrees of ionization. At the atomic level, radiation damage includes the formation of point, line, and volumetric defects: clusters of vacancies or interstitials, dislocation loops, stacking faults, twins, and the formation of stacking fault tetrahedra.^{8} Interestingly, it is well established that several of the line and volumetric defects can be traced at the atomic level to the accumulation of vacancies and interstitial defects due to radiation damage events.^{9,10}

A problem when designing and building a nuclear fusion reactor is selecting candidate materials for its plasma-facing components. The essential constraint is that such materials must resist extreme heat fluxes, together with high fluxes of neutrons, ion beams, and He and H isotopes such as deuterium. To test these materials, one must have access to experimental facilities that reproduce similar conditions to those expected on the materials in magnetic confinement (MC) fusion devices such as the projected ITER tokamak^{11} or inertial confinement (IC) fusion experiments such as those carried out at the National Ignition Facility.^{12,13}

The degree of damage suffered by a material under irradiation depends on several variables, including the interaction time between the plasma and the material, the peak power, and the energies deposited into the material. In practice, these quantities can vary significantly in different environments from MC to IC fusion experiments. In MC, most notably in ITER, energy loads in the divertor associated with *edge localized modes*^{14} of 100 J/cm^{2}–300 J/cm^{2} are expected to have durations of 0.1 ms–0.5 ms, with 10^{3} pulses per shot and frequencies of 0.5 Hz–2 Hz.^{15} In the case of IC experiments, energy loads on the chamber wall of 3 J/cm^{2}–6 J/cm^{2} are expected with durations of 0.2 *µ*s–1 *µ*s and a frequency of around 5 Hz.^{16}

Therefore, it is highly desirable to have a universal measure of radiation damage that characterizes the surface erosion in all these different environments. Different descriptors have been proposed, one being the *integral damage factor* (IDF)^{3,4,17–23} defined as

where *Q* is the power flux, *E* is the energy deposited on an area *A*, and *t* is the time of interaction with the material. There is experimental evidence that a similar degree of deterioration in the material is generated in different environments when the IDF is of the same order of magnitude.^{18} For instance, melting was observed when tungsten was (i) irradiated in the RHEPP ion accelerator with 4.5 J/cm^{2} for 200 ns, (ii) irradiated with the QSPA Kh-50 plasma gun^{24} with 150 J/cm^{2} for 0.5 ms, and (iii) exposed to an energy flux of 550 J/cm^{2} for 3 ms in the JUDITH electron accelerator; in all three cases, the IDF was 7 kW s^{1/2}/cm^{2}–10 kW s^{1/2}/cm^{2}.^{18} Considering the aforementioned conditions expected in ITER and IC experiments, in both cases the expected IDF in the first wall of the reactor is around 10 kW s^{1/2}/cm^{2}.

Among other irradiation facilities and devices (e.g., plasma guns, electron and ion accelerators, pulsed lasers), plasma focus (PF) devices have been used in recent years to probe these conditions because they reach similar IDF values through the high-power flux densities that they generate, thereby enabling material damage to be studied. After pinch compression, a PF device creates a plasma shock that is ejected axially and, when a material is exposed to it, can concentrate energies of 0.01 J/cm^{2}–100 J/cm^{2} for interaction times of 10 ns–500 ns, depending on the distance from the anode to the target.^{21–23,25–27} These conditions correspond to IDFs of 1 W s^{1/2}/cm^{2}–10^{4} W s^{1/2}/cm^{2}. Thus PF discharges can generate an equivalent material IDF to that expected in the first wall of a nuclear fusion reactor (∼10 kW s^{1/2}/cm^{2}).^{18} Although PF devices with storage energies of megajoules were originally used to produce these IDF values,^{28} in recent years PF scaling laws^{29–31} have allowed sub-kilojoule^{32–37} and sub-hundred-joule PF devices^{29,38–45} to be built, and the scaling laws have been extended down to less than 1 J.^{46–49} In particular, work done in the Plasma Physics and Nuclear Fusion Laboratory at the Chilean Nuclear Energy Commission^{30,31} has helped to show that it is possible to scale the PF in wide ranges of energy and size for given ion density, magnetic field, plasma sheath velocity, Alfvén speed, and temperature. Encouraged by these results, table-top PF devices were used to study how the plasma shock ejected axially after the pinch compression affects materials; it was found experimentally that a table-top PF device of hundreds or even only a few joules can produce the same energy densities and IDF values^{23,27,50} as those produced by megajoule PF devices, simply by adjusting the distance between the source and the sample. In this way, it is possible to test candidate plasma-facing materials in a table-top plasma device under the conditions expected in nuclear fusion reactors.

Thus, regardless of the wide varieties of effects and experimental setups, it seems possible to have a global descriptor with which to estimate the degree of damage. In this sense, the IDF *F* appears to be a useful indicator for comparing results from experiments of different types, but unfortunately a proper theoretical explanation of this empirical fact is lacking. Although some arguments have been put forward, there is clearly a need for a sound theoretical foundation for *F* from the atomic phenomena of the creation and recombination of defects.

In the present work, based on the continuous-time master equation, we propose a simple kinetic model that estimates the fraction of defects generated for a given input power flux *Q* and exposure time *t*. This model represents a starting point for understanding how *F* depends on the fraction of defects, giving a physical explanation for this empirical descriptor.

## II. DERIVING MODEL FOR DEFECT FORMATION FROM MASTER EQUATION

We consider a material with *N* atomic sites in area *A*, which defines a surface atomic density of *σ* = *N*/*A*. The state of the sample at time *t* is then characterized by the number of defects *n* out of the total number of sites *N*, which defines a defect fraction *x* = *n*/*N*. However, because the creation and annihilation of defects is a stochastic process, we describe the state of the sample by the probability *P*(*n*|*t*) of having *n* defect sites at time *t*. Alternatively, we can use the fraction of defects *x*.

Our starting point for building a model is the continuous-time master equation^{51,52} for stochastic processes, namely,

which represents local conservation of probability between states. Here, *W*(*n*, *m*) is the transition rate from the state with *n* defects to the state with *m* defects. We assume that defects are created or destroyed one at a time, which implies that transitions occur only between neighboring states, that is,

where *ω*_{+} and *ω*_{−} represent the transition rate amplitudes. Under this assumption, the master equation simplifies to

Now we introduce some assumptions about the *ω*_{+} and *ω*_{−} transition rates. Defects can be created by the incident radiation, but they can also be generated thermally when the temperature *k*_{B}*T* = 1/*β* is sufficiently high, at a rate that is proportional to the fraction of non-defective or ideal sites remaining in the sample. Therefore, the total rate of defect creation must have the form

where *E*_{v} is the energy of defect formation. Furthermore, we know that the energy available to generate defects in time interval Δ*t* and area *A* is *QA*Δ*t*, so the potential number of defects in that area is

In other words, *ω*_{Q} must be proportional to the ratio *Q*/(*σE*_{v}), and we can write

where *η* is a dimensionless factor that accounts for the incident energy that does not create defects, for instance, due to channeling mechanisms.^{53}

Defects are annihilated via recombination, which is also a thermal process and proportional to the current number of defects. This is similar to what is known as a “death process” in the theory of stochastic processes.^{52} We postulate that its rate has the form

Taking the continuous limit by assuming *N* ≫ 1 and replacing *n*/*N* by the continuous variable *x*, we see that the probability density *P*(*x*|*t*) of having a fraction of defects *x* at time *t* follows the continuity equation

with the velocity field *μ* defined as

Taking the first moment of the continuity equation (7) gives (see Appendix A for details) a dynamical equation for the evolution of the average number of defects $xt$, namely,

the solution to which is

where we have defined the constants

In this solution, we now assume that at *t* = 0 the average fraction of defects agrees with the thermal value,^{54,55} given by *x*_{eq}(*β*) = exp(−*βE*_{v}). We can also see clearly that in the limit of no radiation (*Q* = 0), the average fraction of defects must be time-independent, which therefore fixes

in Eq. (10). This implies

Using these conditions and the definitions given in Eq. (11), we obtain *ω*_{T} and *ω*_{R} as

where *τ*_{0} is a characteristic relaxation time of the system. Finally, we obtain the average fraction of defects as

with

where *Q*_{0} = *σE*_{v}/(*ητ*_{0}) is a natural unit of power flux for the system. In natural units, Eq. (13) can be written as

where the exposure time *t* is measured in units of *τ*_{0} and *Q* is measured in units of *Q*_{0}.

Equation (15) is our main result. It describes the evolution of the average fraction of defects as a function of irradiation time *t* starting from a thermal equilibrium state at temperature *T* = 1/(*k*_{B}*β*). It has two free parameters: (i) the temperature-dependent relaxation time *τ*_{0}, and (ii) the fraction *η* of radiation energy used to generate defects. For this expression to preserve the constraint $xt\u22641$ at all times, the inequality

must be met, which naturally gives a range of applicability of the model in terms of the maximum power flux that it can describe.

## III. DISCUSSION

In the following, we explore our model by varying its free parameters to understand how they affect the predicted number of defects in the material and how they are related to the IDF *F*. Varying the power flux *Q* and the material temperature *β*, we estimate the fraction of defects that can be formed for different exposure times *t*. We explore the limitations and range of applicability of our model.

We introduce numerical values for the parameters defined in Sec. II, taking the vacancy-formation energy *E*_{v} and the surface density *σ* as those of bcc tungsten,^{56} namely,

Given these material parameters, we choose reasonable values for the adjustable parameters *η* and *τ*_{0}, which define the value of *Q*_{0}:

Note that the value of *τ*_{0} is on the timescale required for the recombination of mobile defects, which is more than 100 ps.^{57}

### A. Numerical exploration of model

In Fig. 1, we show the evolution of the fraction of defects, given by Eq. (15), as a function of exposure time for different values of material temperature *β* = 1/*k*_{B}*T* and power flux *Q*. At low temperatures (*β* = 10 in units of 1/*E*_{v}), condition (16) allows a wide range of *Q*, indicated by the different dashed blue lines. The fraction of defects approaches zero at *t* = 0 for any value of *Q* because $xt=0=xeq(\beta )\u22480$ for low temperatures. The number of defects increases with the exposure time, reaching a steady state with a maximum number of defects of

for long exposure times. When the power flux is increased, the steady state is reached more quickly because there is more energy available to create defects efficiently, and therefore shorter exposure times are required. The amount of defects reached in equilibrium also increases. However, when the power flux is negligible (*Q* = 0), we obtain the time-independent solution $xt,Q,\beta =xeq(\beta )$, which means that the material is already in the steady state of equilibrium and has a constant number of defects equal to that created by thermal excitation only.

When the temperature is increased for a given power flux (e.g., *Q* = 0.5 *Q*_{0} in Fig. 1), the initial fraction of defects is already greater than zero at *t* = 0 because of thermal excitations. Of course, the initial temperature of the material cannot be arbitrarily large for a given power flux because Eq. (16) places a limit on our model. Given this condition, the maximum temperature that can be explored for a given value of *Q* is *β*_{min} = 1/*k*_{B}*T*_{max} = ln(1 + *Q*)/*E*_{v}, or equivalently, we can only explore power fluxes below *Q*_{max} = 1/*x*_{eq}(*β*) − 1 for a given temperature *β*. This behavior is independent of the time and power-flux scales because these quantities have been normalized to the natural units of *τ*_{0} and *Q*_{0}, respectively.

In previous work,^{27} we showed experimentally that a tungsten sample located 15 mm from the anode top of our PF received a power flux of *Q*^{*} = 9.2 × 10^{6} W cm^{−2} = 4.76 × 10^{−4} *Q*_{0}, which interacted with the sample for *t*^{*} = 75 ns = 250 *τ*_{0}. In Fig. 2, we show the fraction of defects that we predict for this experiment with our model, using a temperature *β* = 150/*E*_{v} (*T* = 300 K), which leads to $x(t*)=4.74\xd710\u22124$. This value is consistent with melting, as we observed in our experiment, because a fraction of defects of between 10^{−4} and 10^{−3} is required to melt a material.^{54,58–60}

In Fig. 3, we show the asymptotic values to which the fraction of defects converges in the limit of long exposure time, $limt\u2192\u221ext,Q,\beta =xmax$, given by Eq. (19). The forbidden region is determined by Eq. (16), and the combinations of *Q* and *β* that lead to the same asymptotic value $xmax$ are indicated by the contour lines. When the power flux *Q* and temperature *β* are such that they lie on the critical (black) line that separates the forbidden region from the region of allowed values, defined by $Qmax=e\beta Ev\u22121$, the fraction of defects will always converge to $xmax=1$. This is, of course, an unphysical situation because having a defect in every single lattice site would violate material stability criteria. In practice, this will not happen because *βE*_{v} < 10 corresponds to temperatures *T* > 4000 K which is well above the melting temperature of most materials. Therefore, reasonable values of *β* are expected to lie far from the forbidden region, which implies that the values of *Q* can span a wide range.

In a recently published article,^{61} the kinematics of lattice defects during irradiation were modeled through the evolution of the distribution functions for different defect types. Our model captures the phenomenology of the concentration of vacancies and interstitials, which converges asymptotically to a constant value with time. According to this study,^{61} aggregates of vacancies are formed more efficiently under irradiation, which is consistent with Monte Carlo simulations and experimental data.^{62} These populations of defects can be measured indirectly through resistivity measurements^{63} and provide valuable information about the nature of their interaction and evolution.

### B. Damage factor *F*

In the quest for a single descriptor for estimating the degree of material damage in different irradiation experiments, the IDF (sometimes known as the *heat flux factor*) defined by

has been proposed.^{3,4,17–23} Linke *et al.*^{18} argued that the same macroscopic damage events (e.g., roughening, cracking, melting) are observed at similar values of *F*, while Fujitsuka *et al.*^{17} suggested that the material surface temperature is proportional to $Qt$, which gives a measure of surface erosion resistance. This hypothesis was later shown to be valid only for constant power density (square pulses) and up to 2500 W s^{1/2}/cm^{2}, so *F* should not be used to describe surface temperature in non-square heat pulses.^{64} In addition to the observations in experiments using plasma guns,^{24} electron and ion accelerators,^{16} and PFs,^{20–23,25–27} melting and roughening were also observed in laser irradiation^{65} for *F* > 3600 W s^{1/2}/cm^{2}. Therefore, over a wide range, *F* seems to work as a simple method for normalizing the thermomechanical effects of irradiation with pulses of different durations.^{66}

From the observation that the same microstructural damage is observed for a given *F* in different experiments, we begin by testing the hypothesis that the average fraction of defects $x$ is a function of *Q* and *t* only through $Qt$; that is,

where *f* is a function to be determined. We can also construct a natural unit for the IDF, namely,

The first point from simple inspection of Eq. (15) is that the IDF hypothesis [i.e., the assumption that $xt=f(Qt)$] does not hold for arbitrary *Q* and *t*. However, as will be seen shortly, interpreting the IDF slightly differently makes it useful for describing material damage by radiation.

To explore the model in terms of the IDF *F*, we consider isolines of fixed $F=Qt$, that is, where the flux *Q* can be expressed as a function of the exposure time *t* by

Replacing *Q*(*t*, *F*) in Eq. (15), we obtain $xt$ as a function of *t* along these contour lines, which are shown in Fig. 4. Here we see that with increasing exposure time *t* for a given *F* (therefore decreasing *Q*), there is an optimum exposure time *t*^{*} that maximizes the average number of defects produced. This is expected because recombination mechanisms dominate for longer exposure times. Interestingly, this maximum value of $xt$, denoted in the following by *x*_{max} and corresponding to an exposure time *t* = *t*^{*}, is in fact determined solely by the value of *F*, i.e., *x*_{max} = *x*_{max}(*F*), as shown in Fig. 5.

Moreover, for *F*/*F*_{0} ≪ 1, this maximum fraction of defects is linear with *F*, as shown in Fig. 6. We can understand this by maximizing $x$ with respect to *t* while replacing $Q=Q(t)=F/t$ with fixed *F*. The resulting maximum average fraction of defects $x*$ is implicitly only a function of *F*. After some algebra, the first derivative of $x*$ with respect to *F* at *F* = 0 is given by

where *u*_{0} is the non-zero solution of the transcendental equation $exp(u02)=1+2u02$ (i.e., *u*_{0} ≈ 1.120 906 423), giving *α* = 0.638 172 686. In this way, for a fraction of defects below 2%, the linear approximation

holds, where *α* ≈ 0.64. For melting to occur, it is well known that a fraction of defects between 10^{−4} and 10^{−3} is needed,^{54,59} so the linear approximation should be good enough in the case of radiation-induced melting. Beyond this linear approximation, the exact equations are given in Appendix B. A more familiar scale for the damage in a material can be produced by defining an “effective temperature” as a function of the maximum fraction of defects with given *F*. For this, we identify

leading to

In other words, $Teffmax(F)$ is the temperature needed to achieve the same fraction of thermal defects in the material as that produced by incident radiation at an IDF *F*. Of course, this scale depends on the material through *E*_{v}, the energy cost of producing a defect. Figure 7 shows this effective temperature computed using Eq. (27) as a function of *F*.

We can also associate an arbitrary fraction of defects with an effective temperature through $xt=exp(\u2212Ev/kBTeff)$ for a given *F*. We show this effective temperature *T*_{eff} as a function of exposure time *t* in Fig. 8 along the lines of constant *F*, where we observe how the effective temperature attains the maximum $Teffmax$ that we discussed in Fig. 7. Because $xt(F)<xmax(F)$ for any exposure time *t* at a given *F*, we have $Teff(F)<Teffmax(F)$. Therefore, for a fixed value of the IDF in Fig. 7, we have a range of smaller possible values of $xt$ that are associated with smaller effective temperatures *T*_{eff}, with $Teffmax$ being the maximum temperature (number of defects) that can be produced for a given *F*. Note that these values are degenerate because, as can be seen in Fig. 4, two different values of *t* are compatible with the same value of $xt$; hence, for a given *F*, two values of *t* are possible for the same effective temperature *T*_{eff}.

## IV. CONCLUSIONS

We have developed a simple model for the kinetics of the fraction of defects as a function of the input power flux *Q* and the exposure time *t*. The model considers the creation of defects due to the incident radiation and also their recombination due to thermal mobility. Interestingly, the model is irrespective of irradiation type (e.g., electrons, ions, plasma shocks, heat flux, electromagnetic radiation, laser light), depending on only the energy of defect formation and not on the specific type of defect. In this way, we obtained a formula [Eq. (15)] for the average number of defects in terms of the incident power flux, the material temperature, and the irradiation time.

We predict that for any given power flux, the material reaches a steady state after a long exposure time, where the fraction of defects stays constant, because of the equilibrium between the creation and recombination of defects. As the material temperature increases, the material reaches a new steady state with more defects. The number of defects that our model predicts is consistent with observations of melting in PF experiments for the same power flux and exposure time, in terms of the expected fraction of defects for melting to occur.

According to the assumptions of our model, the maximum fraction of defects increases as the IDF increases, and the dependence is roughly linear below *F*/*F*_{0} ∼ 0.1, which means that the creation of defects is not as efficient for high values of IDF. Of course, *x*_{max} cannot increase indefinitely, so we may consider that our model is no longer valid above some threshold of defect concentration. Above this threshold, our model breaks down and a different phenomenology, other than defect formation, must be considered.

The description of the IDF from this model supports the observations from experiments of various types with different types of radiation and wide ranges of values for the power flux and irradiation time. More interestingly, it also supports the use of table-top PF devices to test candidate plasma-facing materials under the conditions expected in nuclear fusion reactors. In this way, it is possible to perform these tests in small-scale laboratories, including table-top experiments, which was previously possible only in large experimental facilities.

We also acknowledge that our model does not include other type of defects, such as divacancies or clusters of vacancies, dislocations, or the effect of grain boundaries. However, it constitutes a starting point that allows us to understand the connection between atomic processes and macroscopic measurements of damage induced by radiation in the context of plasma-facing materials, as well as the physical meaning of the IDF descriptor.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the support received from ANID PIA ACT172101, ANID FONDECYT 1171127, and (iii) the bilateral Project No. CONICYT ACE-01 Chile, ANPCyT-PICT2697 Argentina, and CRP IAEA Contract No. 20370. F.G.-C. acknowledges the support received from CONICYT postdoctoral fellowship Grant No. 74160058.

### APPENDIX A: FIRST MOMENT OF ADVECTION EQUATION

We begin with the one-dimensional continuity equation for a random variable *x* ∈ [0, 1] with drift coefficient *μ*(*x*):

Multiplying both sides by *x* and integrating from *x* = 0 to 1, we have

In the first term, we exchange the partial derivative in *t* and the integration, giving

The second term can be solved using integration by parts, namely,

because of the boundary conditions *P*(*x* = 0|*t*) = 0 and *P*(*x* = 1|*t*) = 0, finally giving

### APPENDIX B: NONLINEAR MODEL

We take the average fraction of defects in Eq. (15) and, after defining $u(t)=t$ and eliminating *Q* by using *F* = *u* · *Q*, write it as

with *u* given also implicitly as a function of *F* by the transcendental equation