Mapping between two models of etching process

T.Patsahan1, A.Taleb2, J.Stafiej3, J.-P.Badiali2 1 Institute for Condensed Matter Physics,National Academy of Sciences of Ukraine, 1 Svientsitskii Str., 79011 Lviv, Ukraine 2 Laboratoire d’Electrochimie et Chimie Analytique, ENSCP et Université P. et M. Curie, UMR 7575, 4 Place Jussieu, 75005 Paris, France 3 Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland


Introduction
There is observed an increasing interest not only to the surface growth but also to the etching models.The interest to the growth models is stimulated largely by the molecular beam epitaxy technique [1].Etching has received much less attention and mainly in relation to corrosion problems because of their economical significance [2].A renewed interest in etching comes from the fact that it can be considered as a method of obtaining sophisticated, nanostructured materials by means of a sacrificial surface disaggregation of the bulk material or by anodic dissolution of oxydes [3].The dynamics of the simple growth and etching models has already been studied using the scaling concepts [4].The same type of reasoning has been used in the study of crossover phenomena in the surface dynamics.Thus, a lot of insight in the surface dynamics comes from the scaling concepts [5][6][7][8].
In parallel to this kind of approach, mainly due to physicists, it is interesting to analyze the surface properties from the point of view of chemical-physicists that have developed an enormous literature in the domain of corrosion.The problematic issue is then: to what extent the classical models based on chemical kinetics, i.e., on a mean field approach, can be used in order to analyze the simulation results in a quantitative way.We work on this line and we have shown that in the case of very simple models the interpretation of basic quantities as the mean position of the corrosion front requires a sophisticated analysis [9].Here we consider two etching models.The first model, model A, assumes that the etching probability of a given site depends on the number of nearest neighbours in the solution.The physical idea behind this model is that one effect of the aggressive species present in the solution is to avoid the formation of tips on the surface, therefore we mimic a solution having a strong smoothing character.We compare the results of model A with those of a model B in which the etching property is the same on each site.From a statistical mechanics point of view we have to deal with two different kinds of models: one is non-local while the second one is purely local.It is then reasonable to inquire whether these two models are qualitatively different or not.In order to answer this question we will try to find if model B is appropriate in order to predict the results for model A, i.e., we are looking for some simple relation between to models to fit their results.As we shall see we may consider that the model B can be a mean field approximation of model A.
In the following section we describe the two models more in detail.We define the matching condition for the simpler model.We describe our simulation and compare the results of two models for the matching condition and conclude the paper in the last section.

Two models
In both models we assume that the etched material disappears in a surface reaction.Working at a mesoscopic level we avoid a lot of details associated with atomic level and introduce a simple lattice representation of the reaction (for more details and justification see [9,10]).Thus, we consider that the lattice sites can be occupied by 3 species.The first species is the bulk material of the etched block M .The second species is the etching solution S. By construction of the models the sites M and S cannot be nearest neighbours in terms of the lattice connectivity.They are separated by a layer of surface reactive sites R. The surface reactive sites have at least one nearest neighbour occupied by the etching solution and they are immediately exposed to etching components of the solution.Thus they can undergo the dissolution reaction: R → S which amounts to converting the etched site R to an S site.Now the nearest neighbours of the etched R site have an S site as the nearest neighbour (nn).If these sites happen to be of type M they are immediately converted to the R sites according to: Both models share the above scheme.However, in model A we assume that this probability is proportional to the number n S of S sites being nearest neighbours of the site tested for etching.We write it as follows: where C n is the coordination number for the lattice i. e. the number of nearest neighbours for each lattice site.Thus P ms is the highest possible probability of the R site destruction which happens if the site forms a disconnected island completely surrounded by solution sites and P ms /C n is the lowest possible probability of destruction which happens when the R site has only one S site as its nearest neighbour [10].
In model B the probability of reaction 1 is a constant that we denote P dis .In particular it does not depend on the geometrical situation of the etched site.We write it as follows: Our purpose is to find out whether model A with a linear dependence of P R→S on n S can be simplified to the case of model B where P R→S is constant.To this end we propose a relation converting the parameter P ms of model A to the parameter P dis of model B. Actually, P dis can be presented as an averaged etching probability observed during simulation of model A for the given P ms : where f (n S ) -is a probability to find one of the reactive sites connected to n S solution sites.The distributions f (n S ) are extracted from the simulations of model A for the different parameters of P ms = 0.05, 0.20 and 0.80 (figure 1).Using these distributions the corresponding values of P dis are calculated (table 1).The probability distributions f (n S ) for model B for the adjusted parameters of P dis can been seen in figure 1 as well.From the situations sketched in figure 1 we easily learn that the singular topmost sites (corners) will be more readily etched than the bottommost valley sites leading to a roughness reduction compared to model B.

Simulation details
Using conception of cellular automata technique a two-dimensional lattice with 4 or von Neumann connectivity [11] is simulated.According to above mentioned we determine three types of lattice sites: reactive (surface) sites (R), material (M) and solution sites (S).The simulation box is of the size of N x × N y and the periodical boundary conditions are applied in the x-direction.At every simulation step N s surface sites can be etched with the probability P R→S defined according to the model considered.Due to the site destructions the island formations can be observed.The numbers and the sizes of the islands formed during simulations are collected to build histograms of size distribution.Each simulation run takes from 200000 to 700000 steps depending on the characteristic to be calculated and/or on the parameter of P ms (P dis ).The obtained results are averaged through the number of simulation runs.Due to the time consumption of calculations this number depends on the maximum number of simulation steps in each run.For instance, when the number of simulation steps is 700000 only 50 runs are used, while in the case of the number of simulation steps equal to 200000 we use 300 runs to average our results.
To estimate the absolute values of the obtained characteristics, the averaging through the time in the saturation regime is made, which allows us to achieve a satisfactory accuracy of the most of our results.To collect the islands appearing during simulation and to extract the data about them as well as to remove them from the simulation box, the Hoshen-Kopelman algorithm is applied [12].This algorithm is quite consumptive, but rather general, thus it is efficient for an arbitrary system where intensive islands production is considered.Moreover, this is very convenient in calculating the island sizes which is needed in our study.

Results
To compare the model A and model B two types of roughness are calculated from simulations.The first one is presented as the ratio of the number of surface sites N s to the number of surface sites of the plain surface: Another type of roughness describes the width of the etched front and it is expressed as the root of mean-square deviation of the front position h i from its average position h(t):  In figure 2 the roughness r(t) is presented for the models A and B. For both models this characteristic saturates rapidly after 2000 simulation steps.As it was expected for the model B the higher values of roughness r are observed.The more detailed data for the roughness r are collected in table 2. Also the ratio of the roughness obtained for model B to the roughness obtained from model A is calculated for the parameters considered.It is noticed that this ratio is constant.This result is obtained with a quite high accuracy which is less than 0.4%.An increase of the parameter P ms and the corresponding parameter P dis leads to a decrease of the roughness r, although in the case of model B this decrease is faster.P ms r A r B r B /r A 0.05 1.447 ± 0.003 1.691 ± 0.004 1.169 ± 0.004 0.20 1.403 ± 0.003 1.636 ± 0.004 1.166 ± 0.003 0.80 1.237 ± 0.002 1.443 ± 0.003 1.166 ± 0.003 A width of the front σ cannot be calculated with the same accuracy as for r due to the big oscillations of this characteristic around its average value (figure 3).Also a remarkable difference in time needed for σ saturation is observed for different parameters of P ms and P dis .The lower is the value of P ms (P dis ), the larger is the time to get a saturation for σ.Therefore, for P ms = 0.05 (P dis = 0.0189) the maximum number of steps of simulations are taken equal to 700000 while for P ms = 0.80 (P dis = 0.2964) it is set to 200000.In table 3 one can find the average values of σ with the corresponding ratios between them.Taking into account the accuracies for the calculated ratios which are around 5% it can be concluded that the ratio σ B /σ A is constant and greater than one as well as for r B /r A .The obtained data indicate that the front width for model B is slightly larger than in the case of model A. Also similar to the roughness r, the front width σ decreases when P ms increases.Except roughnesses, a change of the mean position of the etched front with time can be obtained from simulations using the following expression: At large t when a stationary regime is observed h(t) depends linearly on t.A slope of such a dependence describes the dynamics of the etching process and represents a rate of material etching.In table 4 one can see the values of the slopes for the models A and B. The ratio h B /h A between them is not constant, but if one plots these data (figure 4), a strict linear dependence on P ms will be observed.A standard deviation of the linear fit in this case is less than 0.002 and the linearization parameters are derived with an accuracy less than 0.3%.It is worth noting that the slope for the model B is larger than for the model A. It means that a yield per one simulation step for the model B is higher than for the model A and the etching process is faster.Also an increase of the parameter P ms (P dis ) leads to the acceleration of this process for the both models.
Finally, the normalized distributions of island sizes formed during simulations are presented in figure 5.The only maximum of these characteristics appear at the size equal to one.For the model B this maximum is lower than for the model A, therefore more islands with larger sizes can be produced during simulation of the model B. This result could be predicted from the roughness.The higher are the roughnesses r and σ the higher is the probability of the large islands production.Since the roughness decreases with P ms increasing, the number of the islands with sizes larger than one decreases.Thus one can see a rise of the maximum at the distributions considered.The points presented in figure 5 are approximated by functions in the form of y = a/x b .The power b for the given parameters of models A and B vary within the range of b = 1.9 − 3.0.

Conclusions
In conclusion we can say that the two models describe the etching process of a solid material and show qualitatively the same results.From the results one can see that the second model is sufficient to describe the surface roughness and qualitatively the non-locality of the first model plays a minor role in the roughness formation.Also, it is worth noting that for the both models a saturation regime for σ is observed at long times.

Figure 1 .
Figure 1.The histograms of the probability function f (nS) for the models A (clean bars) and B (shaded bars) and for the different parameters of Pms and P dis .

Figure 2 .
Figure 2. Variation of the roughness r versus the number of simulation steps for the models A and B.

Figure 3 .
Figure 3.The evolution of the root of mean square deviation from the average position, σ, for the models A and B.

Figure 4 .
Figure 4.The ratio of the slopes of the mean position of the corrosion front for the models A and B.

Figure 5 .
Figure 5. Normalized distribution functions of island sizes formed during simulations for the model A (closed squares and solid lines) and the model B (open squares and dashed lines).

Table 1 .
The correspondence between some values of parameters of Pms and P dis .

Table 2 .
The roughness r for the models A and B.

Table 3 .
The root of the mean square deviation from the average position σ for the models A and B.

Table 4 .
The slope of the mean position of the corrosion front for the models A and B.