Stochastic simulation of nano-scale surface and near surface region modification: high-temperature blistering and thin films formation
Preprint, Inst. Appl. Math., the Russian Academy of Science

Стохастическое моделирование наноструктурных изменений поверхности и приповерхностного слоя: высокотемпературный блистеринг и формирование тонких пленок

A.L.Bondareva, G.I.Zmievskaya
(Бондарева А.Л., Змиевская Г.И.)

Russian Academy of Science, Keldysh Institute of Applied Mathematics
ИПМ им. М.В.Келдыша РАН

Moscow, 2005

This work is partially supported by Russian Science Support Foundation, RFBR grant № 02-01-01004, the program "Nanoparticles and nanotechnology" by Department of Mathematical Sciences RAS 3.5 and scientific school 1388.2003.2


This paper deals with examination such nanotechnology important processes as blistering and thin films formation. These processes change not only substrate properties, but also plasma behaviour. Phenomena are studied by numerical simulation. Kinetic and stochastic models of fluctuating stage of examined processes are created. The models base on kinetic theory, Brownian motion model and stochastic analog method. Fluctuating stage is very short, but it is very important for development blistering as well as thin film formations. The method of second order of accuracy for solution of stochastic differential equations is modified for examined problems. Received physical results coincide with experimental data and can be basis for laboratory experiments.


Данная работа посвящена рассмотрению таких важных в нанотехнологическом плане процессов как блистеринг и образование тонких пленок на поверхности материала. Эти явления изменяют не только свойства материала подложки, но и свойства пристеночной плазмы. Изучение данных процессов проводится с помощью математического моделирования. Созданы кинетические и стохастические модели флуктуационных стадий рассматриваемых процессов, которые основаны на кинетической теории, модели Броуновского движения и методе стохастического аналога. Флуктуационная стадия является чрезвычайно быстропротекающей во времени, но играет существенную роль для дальнейшего протекания как блистеринга, так и процесса формирования тонких пленок. Модифицирован метод второго порядка точности решения стохастических дифференциальных уравнений для применения к рассматриваемым задачам. Полученные физические результаты совпадают с экспериментальными данными и могут служить основой для проведения лабораторных экспериментов.

Development of high-temperature blistering under the influence of plasma or ion beams on solid surface.



Understanding of processes of interaction between solids surfaces and plasma has significance for development technologies which have concern with cosmophysics, controlled thermonuclear fusion, nanotechnologies. Such problems as blistering, flaking and physical processes, which take place under irradiation of solids surfaces by ions beams. All these processes result in solids properties change, plasma pollution and plasma properties change. Computer simulation of interaction between plasma and solids surfaces takes on special significance in connection with high price, labor-consuming and complexity of laboratory experiments with plasma.

Fluctuation stage of high-temperature blistering is examined in this work. Formation of gas-vacancy pores, which were named blisters, into Ni crystal lattice under the influence of He ions is discussed. He ions have energy from 10 keV to 100 MeV, radiation dose is from 1016 to 1019 ions/sm2, the temperature of sample material T is , Tmelt is melting temperature of Ni. Fluctuation stage is very short-range. Its duration is approximately 10-4 sec. But this stage is defined determinates all following peculiarities of blistering [1-8].


Model of blistering. Kinetic and stochastic equations of model.

          Stochastic model of fluctuation stage of high-temperature blistering has been suggested. The model under discussion is based on Brownian motion model. Blistering is considered as first-order phase transition on its fluctuation stage. Bubbles have size approximately of several angstorm and this kind of defects is considered Brownian particle with sphere form and variable mass. Blisters can interact with each other, with solid lattice and with solid surface. Let us use scheme of splitting on physical processes: bubbles formation and its stochastic motion. The evolution of bubbles presented as a superposition of the stochastic processes of size increase and bubbles stochastic motion in crystalline lattice. It is possible since processes of bubble size increase and bubble migration in lattice have appreciably different time scales. Characteristic time for blisters size increasing is 10-9 sec, typical time for its migration in lattice is 10-8 sec. Kinetic equations of Brownian motion model can be solved by the method of stochastic analog.

The main idea of the method is change of kinetic equations its stochastic analogs and solve stochastic differential equations [9]. Authors apply not only splitting on physical processes but splitting on coordinate too [1-8], because we have taken into account the interaction of all defects with surface. Received kinetic equation and its stochastic analogs display below. Let us present our model by the system of two kinetic partial differential equations of Fokker-Planck-Kolmogorov and Smoluchovskiy-Kramers kinds.

The first equation under consideration (Kinetic Fokker-Planck-Kolmogorov equation for evaluation of blister size) is follow:


g is number of He atoms in bubble,  is distribution function,  is diffusion coefficient in phase space of all possible bubble sizes,  is Gibbs potential of bubble (cluster) formation.

The equation Ito in Stratonovich form, which is equation of stochastic analog of Fokker-Planck-Kolmogorov equation, is presented below

Tk is full time of computation,  is stochastic function, . The same equation can be formulated for vacancies into lattice, but we are examined only He atoms for simplicity of presentation.




We have taken into account the difference between chemical potentials of two phase (), surface tension on bubble –metal surface (b), elastic force of lattice reaction (c) [3-8], inequality between locations in lattice points and internodes (), relcases in crystalline lattice () [3-8]. Nbr is number single relcases,  is energy of single relcase,  is binding energy in lattice. , gcr is critical size of bubble, , . The choosen of initial state of clusterization depends on fluctuation instability of system.


The Gibbs potential has form which is presented on figure below.



The dependence of Gibbs potential (Joule/atom) from the both blister size (g) and position in crystal lattice (r). Blister size is measured in number of helium atoms in bubble. r is measured in lattice parameter of Ni (lattice parameter of Ni a=3.5 Å). The break on g »80 corresponds the first break of due to increasing of blister size.


One-dimensional Gibbs energy versus cluster size. This figure illustrated the region of instability of DF for fluctuation-dependent stochastic process   {g(t), t ³ 0}.



Cluster migration into lattice can be modelized by 3-dimentional Brownian motion. The kinetic equation for motion of blisters is following


   ,   .

We consider that blister comes out surface and dies on it if z£ z2.

Fig. 3 presents examined region schematically.


x,y,z are measured in lattice parameter, zmin is surface, zmax = 2 Rp, Rp is middle depth of projection run, origin of coordinates is in point                          {x=0, y=0, z=0},  is kinetic distribution function,  is diffusion coefficient in lattice space, Mg is mass of blister,  is dissipative factor,   is potential of indirect interaction between bubbles by way of acoustic phonons and Friedel oscillation of electron density in the case of metal substrate [10].

where ar, br, cr are fitting coefficient of model.


The stochastic analog of the kinetic equation for one coordinate is following:

Em is energy of migration, dW is increment of Wiener stochastic process. For z coordinate U=Ubb+Ubs, Ubb is potential of interaction between bubbles and Ubz is interaction potential between bubble and surface.



The fusion of blister is considered in model under discussion. The fusion of two bubbles is made approximately. Two bubbles interflow if following condition is realized: distance between centers of mass bubbles is less than sum of blisters radiuses and some model parameter Df. 0 £ Df £a, a is lattice parameter [3-8].

Fig. 5

Used model allows to find the distribution functions of bubbles from size and pozition in lattice at different moment of time.

Such characteristics as middle size of bubbles, porosity of layers, tension in layers, number of blasted blisters and dependence of these values from time can be find as a result of processing of these distribution functions.
The modified Artem’ev method [11, 3-8] has been used for solution of stochastic differential equations.


The numerical scheme is presented below:


here ADg –operator of size change; ABr –operator of lattice broken; AUD- operator of diffusion, bubbles interaction and interaction between bubble and surface; AFus- operator of fusion; ASurf – operator of exit on surface, reflection from it or destruction on it.


We used 106 trajectories of stochastic process for receiving of physical results, which are presented below.


Two dimensional distribution function, initial time


Blister radius, Å


Distance from surface,

lattice parameter




Two dimensional distribution function,

moment of finish of calculation


Blister radius, Å



Distance from surface,

lattice parameter



Dispersion of blister size

N-number of trajectories




Time, step of algorithm



Fig. 10


Time, step of algorithm


Distance from surface to center of blister,

lattice parameter


Fig. 11


Time, step of algorithm




The porosity of layer with number j is calculated as            blister with number i belongs to j-layer

where V is volume of all sample, Vj   is volume of layer with number j, gi is size of blister at examined time moment, g0 is size of blister at initial time, f(g,z,t) is distribution function at examined time, f(g) is distribution function at initial time, z is distance between surface and centre of blister, z is measured in lattice parameter, z=0 is surface under irradiation.

initial time                                        moment of time finish of calculation

Fig. 12                                                 Fig. 13

Stress in point (x,y,z) is calculated as

z is distance between surface and centre of blister, z is measured in lattice parameter, z=0 is surface under irradiation.

                          The dependence of stress from layer

                   Fig. 14

Stress, s0


Stress, s0


initial time



         time of finish of calculation

                   Fig. 15

Stress, s0


zÎ[0;0.25 Rp]

t=0.25 Tfinish


coordinate y,

lattice parameter



coordinate х,

lattice parameter


                   Fig. 16

Stress, s0


zÎ[0;0.25 Rp]



coordinate х,

lattice parameter


 y coordinate,

lattice parameter


Fig. 17



The dependence of blister size at finish of calculation from different temperatures


Blister size,

Number of He atoms in blister


Temperature, К


                      Fig. 18

Blister size



Change of themperature during calculation


                      Fig. 19



The most interesting results of these experiments are:

1.                 the new kinetic and stochastic model of fluctuation stage of blistering are created;

2.                 the Artem’ev method is modified for solution of system of stochastic equation with functional-coefficients;

3.                 the complex of programs for simulation of fluctuation stage of blistering is elaborated;

4.                 distribution functions of bubbles from sizes and coordinates are nonequilibrium during fluctuation stage;

5.                 blisters chains are formed athwart to incident flux of ions;

6.                 quasi-lattice of bubbles is observed;

7.                 blister size can reach 12 Å during fluctuation stage, decrease of growth rate after bubble size 10 Å is connected with increase of crystal lattice damage;

8.                 bubbles size reaches maximum when relation of solid temperature to melting temperature of materials is 0.47;

9.                 the bubbles migration into the direction of surface under irradiation if bubble radius is less than 5 A, in the other case bubbles stop;

10.             the greatest porosity and greatest tensions are observed on depths of ~ 0.85 Rp and ~ 0.35Rp, Rp is middle depth of projection run;

11.             approximately 15% of blisters destroyed on solid surface during fluctuating stage.


Thin films formation under the influence of plasma or ion beams on solid surface.



          Numerical simulation of adatoms clusterization of solids surfaces under plasma influence is interesting for creation of thin films and covers with necessary behaviour. Release coatings, anticorrosion covers, nano-functional, resistant to pollution and ultrahydrophobic coverings are examples of coatings with defined properties. The study of nano-capsules and ions implantation into near surface layers, interstitial atoms and formed clusters migrations from solids to surface is very important for creation of self-repair materials and coves. The fluctuation stage of thin films formation is examined in this paper.


Model, kinetic and stochastic equations of model.

The modified by ions substrate is choosen metal W. Let us consider this problem using stochastic approach, similar has been used in blistering model, which had been presented in [1-8]. The surface of metal substrate has been contacted with vapour of Ni. For example the thin film material can be examined liquid Ni metal. Thin films behaviour differs from behaviour of solids consisting from same material. Cover formation includes adsorption, creation of new phase islands, increase/decrease of them sizes, motion of new phase islands on surface and others processes. The initial fluctuating stage of thin films formation is of great important, parameters of processes during this stage determine the behaviour of covers in many respects. The duration of this stage is approximately 10-4 sec. The island of new phase (cluster of adatoms) consists of deposited atoms of evaporated material predominantly and small number of implanted atoms which went on surface as result as diffusion on material lattice. Increasing or decreasing of island (i.e. cluster of adatoms) depends on fluctuations during its sizes stochastic changes and jump-like fusion of clusters. Initial stage of cover formation has been considered as the geterogeneous first-order phase transformation (from vapour to liquid on the substrate), here have not be examined chemical reactions during fluctuation stage of this phase transition. Let us model the stochastic diffusion of islands on the surface as a brownian motion adapted for flat coordinate configuration. Following previous experience of computer simulation [1-8] we can use stochastic analog approach as well physical processes splitting. The sizes of islands formation and stochastic migration of islands processes have the different time of development. So, typical time for change of island size 10-8 sec and for migration on substrate surface 10-7 sec. The moving of clusters on surface is realizes under exposure of surface potential and long-range indirect potentials of interaction of clusters each with other.

The phase transition on the surface can be formulated using fundamental Leontovich equation, which is presented by the system of two kinetic equations of Kolmogorov-Feller and Smolukhovskii-Kramers kinds (which have been received after splitting procedure of problem):



where S a is source of vapour which generats ion with fa - maxwell ion function, which is characterized by temperature 2500 K, g is the number of atoms which is consisted in island-clusters,   is the diffusion coefficient in the space of cluster sizes;  is the bubble size distribution function – the probability to find the cluster with size g in interval of values of g  [g,g+Dg],  is the Gibbs energy, Mg is the cluster mass, g is constant of friction, distribution function  is the islands space function  is the position of cluster mass centre in orthogonal coordinates system:  = -200,  = 200,  = -200,  = 200,  is the potential of long-range clusters interaction between them through phonons and oscillation of electron density. The form of potential is similar [3-8], firstly this interaction had been formulated by [10] in problem of light defect clusterization into lattice.


ar, br, cr is model parameters.

The non-linear functional coefficient of equation  is stochastic diffusion coefficient of islands on surface, D0 is coefficient diffusion of adatom on surface, Em is bonding energy of adatom with surface,  is parameter of model.

The Gibbs energy looks like following:

Fig. 20 The scheme of new phase island on solid surface.


where ,  is difference of chemical potential of phases (vapor and liquid), ,  are surface tension between Ni vapour-liquid Ni in island, liquid Ni island-W, Ni vapour-W, c is coefficient of elastic lattice reaction,  is the energy required for breaking of a single bond with lattice, in our case it is value from laboratory experiment,  is the bond energy in lattice, , Nb is number of broken bonds.  shows influence of substrate lattice and the fact that influence of substrate lattice decreases when cluster size increases, when cluster locates in point (x,y) , here ax and ay are lattice parameter on x and y axes, in our case ax=ay=a, a is lattice parameter of W. Y is model function which is depended on islands sizes, also the dislocation of the lattice. If cluster locates in point (x,y) then ~ , otherwise Y~1, the dislocation can be simulated by Y decreasing in same times. All parameters are non-dimensional, it is traditional for kinetic theory. The  and  are nonlinear functional-coefficients which dependence on clusters sizes.

Gibbs energy includes difference of chemical potential of vapour and liquid phases, interface tensions on surfaces of condensate- vapour, condensate- substrate, substrate- vapour, elastic force of lattice and possibility of relcases of part of connections in lattice, non-equivalence of islands positions on surface. The heterogeneous condensation is considered on substrate and on clusters surfaces.


Fig. 21. The  without cluster and dislocation is presented in this figure.

Fig. 22. This figure presents when the cluster has size equated two lattice parameters and located in place with coordinates (3;2). x and y measured in lattice parameters of W.


Fig. 23 This figure presents  near dislocation which placement can be described by equation y=x



For solve these kinetic equations authors used original computational method of stochastic simulation [9, 12, 1-8]. The main idea of this method is using the fundamental qualities of partial differential equations Fokker-Planck kind, which give us possibility to present physical problem by set Ito-Stratonovich equations with functional- coefficients. SDE are equivalent to kinetic problem formulated with Fokker-Planck formalism. We replace of kinetic equation by these stochastic analogs – stochastic differential stochastic Ito-Stratonovich equations /SDE/. The SDE (analogue of equation for fr(g,t)) looks like following:

where Tk is duration of fluctuating stage, x is stochastic function related with increment of Wiener process, g0 is initial cluster size, gmin and gmax are borders of unstable region of initial size of cluster which calculated from , T is temperature of cluster, , gcr is critical size.

For solve of systems of stochastic equations authors modified Artemiev’s method [11]; it is a second-order accuracy method, with infinite domain of stability. For all i=1,2,...106 trajectories of Wiener stochastic process we can use the following determination of the function , where a1 и a2 are random numbers evenly distributed in region (0,1).

Fig. 24 presents of calculation scheme.

Fig. 25. The distribution function of islands from x and y coordinate at initial time moment is shown on this figure. The distribution function is normalize on 1.

Fig. 26. The distribution function of islands from x and y coordinate at finish moment of time is presented on this figure.


finish moment of time


initial time moment


Fig. 27. This figure presents distribution function of clusters from size at initial and finish moments of time. The radius in lattice parameters of W is shown on abscissa axis.

Fig. 28. The snap of surface at finish moment of time is presented in this figure.

Fig. 29. The ratio of total islands square from total islands square at initial time moment is shown. The time in txy is put off on abscissa axis.

Fig. 30. This figure describes change of logarithm of the ratio of total islands square from total islands square at initial time moment depend on logarithm of time. Time is supposed in sec.

Fig. 31. presents the dependence of rate of stress in layer at finish time to stress on surface at initial time moment from layer depth. Z is layer depth measured in lattice parameter of W.


As we can see from presented pictures, radius of islands distribute uniformly from 5.31 Å to 17.7 Å at initial time moment. From 10-4 sec form of distribution function shows that two most probably radiuses exist. The first size (~ 16 Å) is similar to critical size and corresponds to newly form clusters, second size (~ 61 Å) corresponds to islands which grow including at the expense of fusions during calculation.

The number of islands placed near linear dislocation is more than clusters number far from it approximately at 8 times. So, the thin films formation begins on defects of surface such as dislocations.

Three stage of cover formation during fluctuation stage are discovered. The first stage lasts from 0 to 8×10-7 sec, it is stage of slow development. The second stage continues from 8×10-7 sec to 5×10-5 sec and it is stage of quick growth of thin film. The third stage lasts from 5×10-5 sec to 10-4 sec and it is notable for deceleration of growth velocity. At that, cover square increases at 11 times approximately with respect to cover square at initial moment of time.

The calculations confirm that influence of cover reaches on depth of 5 lattice parameters approximately. At the same time, stress on surface and near surface layers caused by thin film formation does not exceed the stress caused by blisters development. The stress on surface connected with cover growth increase at 21 times during fluctuating stage.




The work is partially supported by Russian Science Support Foundation, RFBR grant 02-01-01004, the program "Nanoparticles and nanotechnology" by Department of Mathematical Sciences RAS 3.5 and scientific school 1388.2003.2.



1.     Bondareva A.L., Zmievskaya G.I.// Mathematical Models of Non-Linear Excitations, Transfer, Dynamics and Control in Condensed Systems and Other Media / Ed. L. Uvarova. ¾ N. Y.: Plenum Publ. Co. 1998. P 241-250.

2.     Bondareva A.L., Zmievskaya G.I. // Proceeding. XXV ICPIG, 2001, Nagoya, Japan, Vol.3, pp. 187-188

3.     Bondareva A.L., Zmievskaya G.I. // Russian Academy of Sciences. Physics (Doklady Physics), 2002, V.66, № 7 pp. 994-997

4.     Bondareva A.L., Zmievskaya G.I.// Russian Academy of Sciences. Physics (Doklady Physics), 2004, Vol. 68, No 3, pp. 336-339

5.     Bondareva A.L., Zmievskaya G.I. // Thermophysics and Aeromechanics, Vol. 10, No. 2, 2003, pp. 255-265

6.     Bondareva A.L., Zmievskaya G.I.// Proc. of ICPIG XXVI  July 15-20 2003 Greifswald, Germany. Volume 4. Editors J.Meichsner, D. Loffhagen, H.-E. Wagner, 2003, pp.119-120

7.     Bondareva A.L., Zmievskaya G.I.// Proc. of 16th International Symposium on Plasma Chemistry. Taormina, Italy June 22-27 2003, Edited by R.d'Agostino,  P. Favia, F. Fracassi, F. Palumbo (full-papers CD), 2003

8.     Bondareva A.L. Stochastic simulation of fluctuating stage of high-temperature blistering, Ph.D. thesis, Moscow, 2002

9.     Zmievskaya G.I.// Plasma Physics Reports, 1997. Vol. 23, No. 4, P. 45.

10. Morozov A.I., Ovchankov P.A., Sigov A.S. in: Coll. of Research Papers, ed. by Yu.S. Sigov, Moscow, Keldysh IAM AS USSR, 1990, P. 183.

11. Artemiev S.S., Averina T.A. Numerical analysis of systems of ordinary and stochastic differential equations.-Utrecht, The Netherlands. 1997. P. 176

12. Zmievskaya G.I. // Rarefied Gas Dynamics, AIAA edited by Shizgal, B.D. and Weaver, D.P. Vol. 159 Washington, DC, 1994, pp. 371-383