# Levy flight of atoms in collision cascades

David Simeone, Laurence Luneville, Jean-Pierre Both

► To cite this version:

David Simeone, Laurence Luneville, Jean-Pierre Both. Levy flight of atoms in collision cascades. EPL - Europhysics Letters, 2008, 83, pp.56002. 10.1209/0295-5075/83/56002. cea-05467177v2

HAL Id: cea-05467177

https://cea.hal.science/cea-05467177v2

Submitted on 11 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Levy flight of atoms in collision cascades

D. Simeone¹(a), L. Luneville² and J. P. Both³

¹CEA, DEN, SRMA, LAZM, Equipe Mixte MFE - F-91191 Gif-sur-Yvette, France, EU

²CEA, DEN, SERMA, LLPR, Equipe Mixte MFE - F-91191 Gif-sur-Yvette, France, EU

³CEA, LIST, SSTM, LPSS - F-91191 Gif-sur-Yvette, France, EU

**Abstract** – In this work we analyze the ballistic phase of displacement cascades. Monte Carlo simulations performed with realistic cross-sections were used to calculate the probability density function, $pL(l)$, describing the length traveled by a particle between two successive collisions within the binary collision approximation framework. From the analysis of $pL(l)$, this work shows that the trajectories set displays a bifractal behavior. The comparison of different MC simulations clearly shows that this bifractality is only due to the long tail of $pL(l)$. The unusual form of $pL(l)$ implies that atoms displaced in the collision cascade during irradiation follow a Levy flight. Since the Levy law is scale invariant, an anomalous diffusion of atoms at large distances is then expected under irradiation and may explain experimental results.

**Introduction.** – The ion solid interaction is of significant interest to both academic and industrial researchers. Ion implantation revolutionized the microelectronic industry offering a control over the number and depth of dopant atoms in semiconductor materials. One of the main characteristic features of ion irradiation [1,2] is an athermal diffusion of atoms during the penetration of an incident ion in a solid. During the first step of a collision cascade, incident particles with a kinetic energy below few MeV transfer mainly their energy to the target atoms through elastic collisions. As the kinetic energy of knocked-on particles remains large, most of the collisions occurring in the target can then be treated as separate binary collisions and the binary collision approximation (BCA) can be applied [3,4]. At the end of the collision cascade, the kinetic energy of recoils is always small. The collisions get a more collective behavior with many atoms suffering small random displacements in a limited volume called the displacement spike [5,6]. The molecular dynamics (MD) method is nowadays extensively used to describe the slowing-down of particles in the displacement spike [7,8]. MD simulations allow to calculate the spontaneous recombination volume and then the amount of stable defects. However, MD simulations of a displacement cascade initiated by high-energy particles (above few keV) is very time consuming [8]. The contribution of high-energy events in a cascade is therefore out of reach in MD. This analysis points out that the BCA remains then an useful theoretical framework to study the migration of atoms at large distance generated during a displacement cascade. Many models first developed in the 1960s by Lindhard, Scharff and Schiott (LSS theory) [3] and in the 1970s by Winterbon, Sigmund and Sanders (WSS theory) [9] were then built up to describe the displacement of atoms within the BCA using simple interatomic potentials. According to these models [10–12], atoms involved in a cascade display a Gaussian displacement distribution at low temperature, where defect-mediated diffusion is negligible. A ballistic diffusion coefficient can then be derived to describe the effect of irradiation on the diffusion of atoms in solids. However, experimental profiles induced by irradiation of light targets by high-energy particle, i.e. few thousands of keV, exhibit large tails [13–15] which cannot be understood within this diffusion approximation framework.

The aim of the present work is to estimate the impact of the displacement of atoms at large distance on the set of trajectories in the displacement cascade and on the diffusion of atoms at large distances. To reach such a goal, the universal Biersack-Ziegler-Littmark (BZL) potential, extensively used to describe the interatomic pair potentials at short distance in MD, allowed to build trajectories of atoms in motion in a displacement cascade in a realistic way within the BCA.

In the first part of this work, we describe in detail the elementary probability density functions used in our MC simulations. In the second part, the fractal property of the trajectories set extracted from these MC simulations is then analyzed within the multifractal formalism framework. In the last part of this letter, we then study the diffusion of atoms at large distances induced by collision cascades, using the fractal properties of collision cascades.

**Monte Carlo simulation of a displacement cascade within the BCA.** – To generate trajectories of atoms in motion in a displacement cascade, the energy transfered between atoms during a collision as well as the lengths traveled by atoms after a collision need to be known. Within the BCA, the probability density function for an incident particle with energy $x$ to loose an amount of energy ($x - y$) during a collision is equal to $k(x, y) = d\sigma(y,x)/dy * 1/\sigma(x)$, where $d\sigma(y, x)$ and $\sigma(x)$ are the differential and total cross-sections, respectively. The function $k(x, y)$ is no more than the conditional probability of the distribution of $y$ provided $x$. In this work, MC simulations are performed using a realistic BZL differential cross-section derived from the universal Biersack-Ziegler-Littmark interatomic potential. This approximation permits to define accurately the energy transferred to recoils over a large energy range [16,17]. The strength of this inter atomic potential depends only on the average ionic charge, $Z$, between incident particles and atoms at rest in the solid. Moreover, simpler interatomic potentials were also used to compare our results with previous works [6,10,18,19]. Even if these approximations describe crudely the energy transferred between atoms during collisions, they can be considered as guidelines to analyze MC simulations performed using the BZL approximation. The hard-sphere (HS) differential cross-section has often been used as a convenient starting point to study the slowing-down of atoms in motion in a displacement cascade [6]. Using such a toy model, the energy lost by atoms during a collision is analytically calculated. However, $k(x, y)$ does not evolve with $y$ for this approximation. The power law (PL) differential crosssection in conjunction with the inverse power potential gives a more realistic description of the energy transferred by an incident particle to an atom at rest in the solid because $k(x, y)$ depends on $y$ as expected. However, the PL differential cross-section describes the transfer of energy between atoms only in a limited range of energy since the strength of the interatomic pair potential $m$ which describes the binary collision is a constant [5].

To estimate particle trajectories in a displacement cascade within the BCA, the length traveled by a particle after a collision, $l$, needs to be calculated. Neglecting memory effects from a collision to the next due to spatial correlation between scattering atoms, the conditional probability density of $l$ for a given $y$ value is an exponential law [20]:
$$p(l|y) = N \sigma(y) e^{-N \sigma(y) l}, \quad (1)$$

where N is the number of scattering atoms per unit volume. In all our MC simulations, the length l traveled by a particle of energy y after a collision is sampled using the inversion method from eq. (1). Moreover, the function k(x, y) allows to sample the energy y lost by a particle during a collision. Since no analytical formulation for the BZL differential cross-section is available, the rejection method [21] is used to sample y. For the HS and PL approximations, analytical expressions for differential cross-sections exist:

$$d\sigma(y, x) = C_m x^{-2m-1} f((x - y)/x) d(x - y), \quad (2)$$

where Cm and m are constants independent of the energy. For the HS approximation, $f((x - y)/x)$ is equal to 1 and y varies from 0 to x. For the PL approximation, $f((x - y)/x)$ is equal to $((x - y)/x)^{-m-1}$ and y varies from 0 to $x(1 - \epsilon)$ where $\epsilon$ is a cut-off factor insuring a finite value for the total crosssection. As $k(x, y)$ can be analytically calculated from these differential cross-sections, y can be sampled using the inversion method for these two approximations [21].

This double sampling in energy and in length allows then to determine accurately the trajectory of each displaced atom in the displacement cascade. The simulation stops when the energy of all atoms in motion decreases below the displacement energy threshold, xd. In all our MC simulations, the numerical values Cm, N, and $x_d$ were equal to $1 \cdot A^2 \text{eV}^{2m}$, $1 \cdot A^{-3}$ and 10 eV, respectively. For the PL approximation, MC simulations were performed using different values of $\epsilon$ ranging from $10^{-5}$ to $10^{-3}$. Only small l values evolve slightly as a function of $\epsilon$. However, the nature of the trajectories set does not depend on $\epsilon$. For this reason, $\epsilon$ was set to $10^{-4}$ in this work.

To investigate the impact of the strength of the interatomic pair potentials on the geometric properties of the trajectories set, MC simulations were performed. In these simulations, Z then varies from 1 to 99. In the same way, MC simulations were performed with different m values, corresponding to different energy ranges [7] for the HS and PL approximations. For that, m varies from 0 to 1 to scan both low- (below 10 keV) and high-energy (above few hundreds of keV) domains.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Comparison of $p_L(l)$ extracted from MC simulations performed with $x_0 = 0.5 \text{ MeV}$ and $m = 1$ for HS (black triangles) and PL (open dots) approximations. MC simulations were performed with $Z = 2$ for the BZL approximation (black stars). The functions $p_L(l)$ calculated for different cross-sections follow the same power law over many decades (dashed line).</figcaption></figure>

**Multifractal analysis of collision cascades.** – The pattern of the trajectories set extracted from MC simulations is composed of separated domains in which the main part of collisions occurs and exhibits a fractal property. Previous authors have already observed the fractal behavior of trajectories in a displacement cascade [6,10,16–19]. Cheng was the first one to point out the fractality of a deterministic idealized collision cascade using the HS approximation [6]. However, the random nature of both y and l in its idealized displacement cascades was neglected. To describe more accurately the random nature of y, different authors analyzed the geometric properties of the displacement cascades using the PL and the BZL approximations [16–19]. However, they did not take into account the random nature of l. Some authors take into account the random nature of l. In their works, the probability density function describing the evolution of l is not linked to the physics of binary collisions [16,17]. In our work, eq. (1) describes the density probability function of l as a function of the energy y of the particle after a collision without introducing any cut-off for l. Since the length distribution is the physical parameter of the problem, the total cross-

section $\sigma(x_0)$ of incident particles of energy $x_0$ and the number of atoms per unit volume in the solid, N, are used to define the characteristic length of the problem $L = 1/(N\sigma(x_0))$. Summing over all energies transferred to atoms of the medium during the collision cascade, the probability density function associated with the length traveled by an atom between two successive collisions, $p_L(l)$, can be derived from eq. (1):

$$p_L(l) = \int_{x_d}^{x_0} p(l|y)g(x_0, y)d y,$$

where $g(x_0, y)$ is the probability density function associated with the number of atoms in motion with an energy y generated by an incident particle of energy $x_0$ in a displacement cascade [22].

Figure 1 shows the evolution of $p_L(l)$ as a function of l L derived from MC simulations. The functions $p_L(l)$ exhibit a large tail over many decades. Since the function $x_d/\gamma^2$ is equal to $g(x_0, y)$ for the HS approximation and fits $g(x_0, y)$ over a large energy range for the PL approximation, $p_L(l)$ was analytically calculated from eq. (3). The function $p_L(l)$ is then equal to $b/(l^{1+q_c})$ over a large length range, where $q_c$ is a numerical constant and $b = x_d/x_0 q_c r(1 + q_c)L^{q_c}$ for these two approximations. Since the PL cross-section describes on a first approximation the cross-section derived from the BZL interatomic potential, the function $b/(l^{1+q_c})$ is also used to approximate $p_L(l)$ extracted from MC simulations for the BZL approximation. Figure 1 clearly displays then that $p_L(l)$ is a power law of l over a large range for all the approximations (dashed line) since l for a given y value is sampled by eq. (1).

<figure><img src="image_3.png" /><figcaption>Fig. 2: Variation of $Z(q, L)$ extracted from MC simulations performed with different incident energies as a function of L using $m = \frac{1}{2}$ for HS (black triangles) and PL (open dots) approximations, and $Z = 15$ for the BZL (stars) approximation. The linearity of these curves insures the fractality of the trajectories set.</figcaption></figure>

To study in detail the impact of the large tail of $p_L(l)$ on the fractal property of the trajectories, we introduce the normalized moments $Z(q, L)$ [23] describing the different moments of $p_L(l)$:

$$Z(q, L) = \frac{\int_{0}^{+\infty} p_L(l) l^q dl}{N_c^{-1} \left( \int_{0}^{+\infty} p_L(l) l dl \right)^q}, \quad (4)$$
where $q$ is a real number varying from $-\infty$ to $+\infty$ and $N_c$ is the total number of collisions occurring in the cascade. $N_c$ is extracted from MC simulations.

Figure 2 shows the evolution of $Z(q,L)$ as a function of $L$ for different $q$ values. For large $L$ values, $\text{Log}(Z(q,L))$ is proportional to $\text{Log}(L)$. This point insures that the trajectories set clearly possesses a fractal behavior [24] for all the approximations. This fractality of the trajectories set is not linked to the peculiar form of the crossoctions. Moreover, fig. 2 clearly shows that the slopes of the different lines do not depend on the energy of the incident particle but vary only with $q$. This variation of the slopes with $q$ is the signature of the multifractality of the trajectories set [24]. From this analysis, we deduce that $Z(q,L) \sim L^{(q)}$ for large $L$ values, i.e. for irradiations by few hundreds of keV particles.

To analyze the multifractality of collision cascades, the evolution of the exponent $\tau(q)$, called the mass exponent [24], is plotted as a function of $q$ in fig. 3. The evolution of the mass exponent vs. $q$ is composed of two straight lines for all the approximations. Since $\text{pl}(l) \sim b/(l^{1/qc})$ over many decades, as shown in fig. 1, $\tau(q)$ exhibits two different slopes defined as the Holder exponents [24] for $q < q_c$ and $q > q_c$ insuring the bifractal

<figure><img src="image_4.png" /><figcaption>Fig. 3: Evolution of the mass exponent $\tau(q)$, extracted from MC simulations as a function of $q$ for the different approximations (black triangles for HS with $m = \frac{1}{2}$, open dots for PL with $m = \frac{1}{2}$, stars for BZL with $Z = 90$). The kink of the two slopes at $q_c$ is the proof of the bifractal property of the trajectories set.</figcaption></figure>

behavior of collision cascades. A kink of the slopes appears then clearly for $qc$.

<table><caption>Table 1: Calculation of different coefficients $\tau(q)$, $q_c$ and the two Holder exponents $\alpha_1$ and $\alpha_2$ defining the fractality of the trajectories set for different approximations. For the HS and PL approximations, these coefficients can be explicitly calculated and depend only on $m$. For the BZL approximation, these different coefficients extracted from MC simulations are fitted to functions of $\text{Log}(Z)$. For $q > q_c$, a constant $Z_c$ equal to $0.6 \pm 0.05$ is needed to fit $\tau(q)$. The correlation coefficients, $\rho$, assessing the quality of fits are equal to 0.99 for all fits.</caption><thead><tr><th></th><th>HS and PL</th><th>BZL</th></tr></thead><tbody><tr><td>$q_c$</td><td>$\frac{1}{2m}$</td><td>$\frac{1}{2}\text{Log}(\frac{Z}{Z_0})$</td></tr><tr><td>$\tau(q) (q < q_c)$</td><td>$\frac{1-q}{2m}$</td><td>$\frac{1-q}{2}\text{Log}(Z)$</td></tr><tr><td>$\tau(q) (q > q_c)$</td><td>$q - \frac{q}{2m}$</td><td>$q - \frac{q}{2}\text{Log}(Z) + \frac{\text{Log}(Z_0)}{2}$</td></tr><tr><td>$\alpha_1$</td><td>$\frac{1}{2m}$</td><td>$\frac{1}{2}\text{Log}(Z)$</td></tr><tr><td>$\alpha_2$</td><td>$\frac{1}{2m} - 1$</td><td>$\frac{1}{2}\text{Log}(Z) - 1$</td></tr></tbody></table>

This analysis shows that the origin of the bifractality is only due to the particular power law form of $\text{pl}(l)$ as in studies on aggregation processes [25]. For the HS and PL approximations, $\tau(q)$, $qc$ and the two Holder exponents can be analytically calculated. These functions do not depend on the energy and are only functions of $m$ which characterizes the strength of the interatomic pair potentials. For the BZL approximation, the same analysis was performed and the evolution of these functions is also independent of the energy in agreement with MC simulations performed in this work. As expected from the fractality, these functions evolve only with $Z$ which defines the strength of the BZL interatomic potential. Using empirical approximations for the BZL interatomic pair potential [26], the strength of the BZL differential crossoction is roughly a function of $\text{Log}(Z)$. Such a function has then been used to fit the evolution of $\tau(q)$, $qc$ and the two Holder exponents for this approximation. Table 1 summarizes the different variations of $\tau(q)$, $qc$ and the two Holder exponents for all the approximations.

The mass exponent $\tau(q)$ allows also to define the generalized dimensions $D(q) = \tau(q)/(1-q)$ of the trajectories set [24,27]. From $Z(0, L)$, the LSS path length [18,19,28], Ltot and the number of collisions $N_c$ previously used to point out the fractal behavior of the trajectories in the collision cascades [6,10] can be easily calculated. In our work, these two values are proportional to $L^{D(0)}$. The fractal dimension associated with Ltot and $N_c$ is then $D(0) = 1/(2m)$ in agreement with the fractal dimension of idealized cascade displacements already determined in previous works [6,10,18,19,28] for the HS and PL approximations. However, fig. 3 clearly shows that the fractal properties of collisions cascades can then not be defined by an unique fractal dimension [6,19,28]. Moreover, the fractal behavior of displacement cascades dictates that the fractal dimension $D(0)$ must not depend on the energy and cannot be used as a criterion to define the sub cascade formation [6,29].

**Diffusion of atoms at large distances.** – The probability density function, $\text{pn}(z|y_0)$, associated with the total length, $z$, traveled by a particle with an initial energy $y_0$, after $n$ collisions can be written as [30]

$$p_n(z|y_0) = \prod_{j=1}^{n} \int_{0}^{+\infty} \int_{0}^{y_{j-1}} k(y_{j-1}, y_j) p(l_j|y_j) d\mathcal{l}_j \delta\left(z - \sum_{j=1}^{n} l_j\right), \quad (5)$$

where $l_j$ and $y_j$ are the length traveled by the incident particle and its energy after the j-th collision, respectively. The function $k(x, y)$ insures that the lengths $l_j$ are not independent. This equation describes a non-Markovian random walk process with correlation between the lengths of the subsequent steps.

The function $\text{pn}(z|y_0)$ is similar to the relocation function in the WSS theory [31]. To take into account the displacements of all atoms in the displacement cascade, the function $\text{pn}(z|y_0)$ needs to be averaged by $g(x_0, y_0)$, the amount of atoms set in motion directly or indirectly by the incident particle. This procedure allows to introduce the probability density function, $\text{pn}(z)$, associated with the total length, $z$, traveled by a particle after $n$ collisions whatever its trajectory:

<figure><img src="image_5.png" /><figcaption>Fig. 4: Evolution of $p_{100}(z)$ as a function of $z$ extracted from MC simulations. Simulations were performed with incident particles of energy 0.5 MeV and $m=1$ for the HS approximation (black triangles) and for the PL approximation (open dots). For the BZL approximation (black stars), MC simulations were carried out with $Z=2$. For all the approximations, $p_{100}(z)$ tends to a Levy law and exhibit long tails proportional to $\frac{1}{z^{1/2}}$ (dashed line).</figcaption></figure>

This coarse-graining procedure insures that $\text{pn}(z)$ depends only on the energy of the incident particle. For $Z > 32$ or $m < 1/4$, i.e. $qc > 2$, $(l)$ and $\sigma(l)$, the mean value and the standard deviation of $\text{pl}(l)$ exist. From the central limit theorem [30], $u_n = (zn-n(l))/(\sqrt{n\sigma(l)})$ tends in law to a Gaussian function, $\text{LG}(0, 1)$. On the other hand, for $Z < 32$ or $m > 1/4$, i.e. $qc < 2$, these values are large and $u_n = zn/(n^{1/qc})$ or $u_n = (zn-n(l))/(n^{1/qc})$ tend in law to a Levy function $\text{Lb}(qc, 0(u))$ [32] for $qc < 1$ and $qc > 1$, respectively. Figure 4 clearly displays that $p_n(z_n)$ tends to a Levy function $p(z)$ and justifies this analysis. Assuming that all the $y_j$ are independent random variables, eq. (6) reduces to a convolution
product of n independent variables associated with the same probability density function pL(l). For $qc < 2$, the density probability function associated with zn tends then to a Levy law $p(z)$ [30]. Since $k(x, y)$ is a decreasing function of y, the random variables yj and $y_{j+r}$ can be considered as independent for large r values. This point insures that zn follows a Levy flight as shown in fig. 4.

From this analysis, it appears clearly that replacing eq. (1) by its mean value [6,10,28,29] introduces an artificial cut-off for pL(l). This cut-off forbids the appearance of the large tail of pL(l) and then the description of p(z) by a Levy law.

The function $p(z)$ ties the geometry of displacement cascades with the ability of a displacement cascade to shift atoms. Since $g(x_0, y)$ can be determined from MC simulations for all our approximations, the deviation angle between the trajectory of each atom in motion and the incident beam, $\theta$, can be evaluated within the BCA. The probability for $\theta$ to be above $1^\circ$ is then inferior to 10% in the ballistic phase of collision cascades induced by incident particles of energy $x_0$ above few keV. Since only highenergy particles with a small deviation angle contribute to significative z values, z is then roughly equal to the projected range. The function $p(z)$ describes the physics of collision cascades and determines the long time evolution of an atomic profile, $c(z, t)$, under irradiation. This equation is analogous to the Kolmogorov-Fokker-Planck equation used to describe the macroscopic kinetic behavior of macrovariables in chaotic processes [33]. For a given experiment, the nature and the energy of particles in motion define $qc$. From this analysis, $qc$ plays the role of a critical parameter [34] for this kinetic equation. Neglecting sputtering effects and overlapping of collisions cascades [35], the effect of an homogeneous irradiation at constant flux, $\phi$, on $c(z, t)$ is given by a one-dimensional balance equation for $t>tc \sim 10^{-9}$ s, the last of the collision cascade:

$$ \partial c(z, t) \partial t = \phi \sigma(x_0) + \infty \int_0^\infty p(\Delta)c(z - \Delta, t) - p(\Delta)c(z, t)d\Delta, c(z, 0) = \delta(z - 0). \quad (7) $$

Since $p(z)$ displays two distinct behaviors as a function of $qc$, the evolution of $c(z, t)$ as a function of $qc$ needs to be studied. For $qc > 2$, $p(z)$ is a Gaussian law and eq. (7) is reduced to a diffusion equation. The concentration profile $c(z, t)$ follows also a Gaussian law. The second moment of $c(z, t)$, defined by $\langle |z|^2; t \rangle$, exists and is proportional to the fluence $\Phi = \phi t$ in agreement with experimental data [10,12]. A diffusion coefficient proportional to the flux $\phi$ can then be extracted from this Gaussian profile. For $qc < 2$, the leading term of $p(z)$ is $b/z^{1+qc}$ for large z values. For $(z/L) > (x_d/x_0)^{1+qc} \sim 10^{-2}$, $c(z, t)$ is then approximated by a Levy law $L^{\alpha(x_0)\Phi b}_{qc,0} (z/t^{1/qc})$. The function $c(z, t)$ exhibits then a large tail which could explain large tails of impurities profiles measured at low temperature in many experiments [13–15]. From the scaling properties of the Levy law, the r-moments of $c(z, t)$, $\langle |z|^r; t \rangle$ can be calculated for $r<qc< 2$. These moments are proportional to $\Phi^{r/qc} \langle |z|^r; 1 \rangle$. It appears then that $\langle |z|^2; t \rangle$ is large and proportional to $\Phi^r$ with an Hurst exponent [24], $H = 2/qc$, always greater than 1 for long irradiation times. This point insures that $\langle |z|^2; t \rangle$ is no more a linear function of $\Phi$ for $Z < 32$. We expect then that an anomalous diffusion must then take place in mixing experiments when light targets, $Z < 32$, are irradiated by few hundreds of keV particles at low temperature. Since the diffusive mixing occurring in the displacement spikes that dominates at the end of atom trajectories takes place over small distances largely inferior to $L(x_d/x_0)^{1/qc}$ [8], this anomalous diffusion should then not be washed out during the displacement spike and could explain experimental results [13-15].

**Conclusion.** – In this work, we investigate the geometric properties of collision cascades within the binary collision approximation using a physically sounded function to handle the physics of binary collisions. The universal BZL interatomic pair potential permits then to perform accurate MC simulations of the trajectories set taking into account the displacement of atoms at large distances. Introducing the natural length for the problem, L, we clearly show that the probability density function associated with the length traveled by a particle between two successive collisions, pL(l), exhibits a large tail. This work clearly ties the long tail of pL(l) to the fractal nature of trajectories. Applying the multifractal analysis formalism, a study of the generalized moments insures that the trajectories set is bifractal. MC simulations performed with different cross-sections assess that this bifractality does not depend on a differential peculiar form of cross-sections. The long tail of pL(l) is only responsible for this bifractality. This bifractality of collision cascades permits to establish a macroscopic kinetic equation for the evolution of an atomic profile $c(z, t)$. For $Z < 32$, this profile follows a Levy law. This unusual form of $c(z, t)$ could explain the long tail of experimental profiles observed in many experiments [13–15]. From the self-similarity property of the Levy laws, we expect, in this case, an anomalous diffusion of atoms at large distances. Such an anomalous diffusion could be evoked to explain the unusual diffusion of impurities in nanometric materials under irradiation [36]. Moreover, its impact on the calculation of the ballistic jump frequency, extensively used to study the stability of solids under irradiation, has to be taken into account [37].

***

The authors thank Y. Limoge and E. Bouchaud for many valuable discussions and continuous encouragements.

## REFERENCES

[1] Martin G. and Bellon P., Solid State Phys., 53-54 (1997) 1.

[2] Krasnochtchekov P., Averback R. S. and Bellon P., Phys. Rev. B, 75 (2007) 144107.

[3] Lindhard J., Schraff M. and Schiott H. E., Mat.-Fys. Medd. K. Dan. Vidensk. Selsk., 33 (1963) 1.

[4] Lindhard J., Nielsen V. and Schraff M., Mat.-Fys. Medd. K. Dan. Vidensk. Selsk., 36 (1968) 1.

[5] Sigmund P., Sputtering by Particles Bombardment (Springer Verlag) 1981.

[6] Cheng Y. T., Nicolet M. A. and Johnson W. L., Phys. Rev. Lett., 58 (1987) 2083.

[7] Smith R., Jakas M., Ashworth D., Oven B. and Bowyer M., Atomic and Ion Collisions in Solids and at Surface (Cambridge University Press) 1997.
[8] Averback R. and Diaz de la Rubbia T., Solid State Phys., 51 (1997) 281.

[9] Winterbon K. B., Sigmund P. and Sanders J. B., Mat.-Fys. Medd. K. Dan. Vidensk. Selsk., 37 (1970) 1.

[10] Cheng Y., Mater Sci. Rep., 5 (1990) 45.

[11] Kacsich T., Weber T., Bolse W. and Lieb K., Appl. Phys. A, 57 (1993) 187.

[12] Bolse W., Mater. Sci. Eng., 12 (1994) 53.

[13] Besenbacher F., Bottiger J., Nielsen S. and Whitlow H., Appl. Phys. A, 29 (1982) 141.

[14] Banwell T., Liu B., Golecki I. and Nicolet M., Nucl. Instrum. Methods, 209-210 (1983) 125.

[15] Wang Z., Westendorp J. and Saris F., Nucl. Instrum. Methods, 209-210 (1983) 115.

[16] Kun F. and Bardos G., Phys. Rev. E, 55 (1997) 1508.

[17] Kun F. and Bardos G., Phys. Rev. E, 50 (1994) 2639.

[18] Winterbon K., Radiat. Eff., 60 (1982) 199.

[19] Winterbon K., Urbassek H., Sigmund P. and Gras-Marti A., Phys. Scr., 36 (1987) 689.

[20] Amsel G., Basttistig G. and L'Hoir A., Nucl. Instrum. Methods B, 201 (2003) 325.

[21] Rubinstein R., Simulation and the Monte Carlo Methods (John Wiley and Sons) 1981.

[22] Sigmund P., Rev. Roum. Phys., 176 (1972) 1079.

[23] Halsey T. C., Jensen M. H., Kadanoff L. P., Procaccia I. and Shraiman B. I., Phys. Rev. A, 33 (1986) 1141.

[24] Feder J., Fractals (Plenum) 1988.

[25] Amitrano C., Coniglio A. and di Liberto F., Phys. Rev. Lett., 57 (1986) 1016.

[26] Ziegler J. F., Biersack J. P. and Littmark U., The Stopping and Range of Ions in Solids (Pergamon Press) 1985. [27] Grassberger P., Phys. Lett. A, 97 (1983) 227.

[28] Moreno-Marin J., Conrad U., Urbasseck H. and Gras-Martín A., Nucl. Instrum. Methods B, 48 (1990) 404.

[29] Rossi F., Parkin D. and Nastasi M., J. Mater. Res., 4 (1989) 137.

[30] Feller W., An Introduction to Probability Theory and its Applications (John Wiley and Sons) 1966. [31] Sigmund P. and Gras-Marti A., Nucl. Instrum. Methods, 182-183 (1981) 25.

[32] Bouchaud J. P. and Georges A., Phys. Rep., 195 (1990) 127.

[33] Shlesinger M., Zaslavsky G. and Klafter J., Nature, 363 (1993) 31.

[34] Frisch U., Turbulence and Predictability in Geophysical Fluid Dynamics and Climat Dynamics, edited by Ghil M., Benzi R. and Parisi G. (North Holland) 1985.

[35] Collins R., Marsh T. and Jimenez-Rodriguez J. J., Nucl. Instrum. Methods, 209-210 (1983) 147.

[36] Meldrum A., Boatner L. A. and Ewing R. C., Phys. Rev. Lett., 88 (2001) 25503.

[37] Enrique R. A. and Bellon P., Phys. Rev. Lett., 84 (2000) 2885.