*3.1. Convection–Di*ff*usion 2D Model for the Stochastic Coating Deposition Process*

The creation of a wear resistant layer of the coating combines two parallel processes: convection and diffusion. Convection causes the rectilinear motion of particles from a cathode to the surface on which a coating is to be deposited, while diffusion gives rise to the consideration of random factors influencing particle trajectories, as a result of which the density of the substance being deposited becomes inhomogeneous, and the process represents the diffusion limited by the clustering or, alternatively, the clustering limited by the diffusion. A 2D model of the coating deposition process was considered, according to which the area where particles were deposited was plane (Figure 1). This assumption can be quite logically extended to a 3D model, because the flow of the particles moving through the fixed plane parallel to the 0*XZ* plane is quasi-isotropic. The typical linear size of the deposition area is about *L*. The particles are emitted from the line of *y* = *b* and are deposited on the line of *y* = *a*.

**Figure 1.** 2D model of the cross-section of the area perpendicular to the plane in which the particles emitted from the cathode in the direction of the substrate are being deposited.

*Coatings* **2020**, *10*, 927

In case if no diffusion occurs, a particle leaving a point with the coordinate (*x*1, *b*) enters a point with the coordinate (*x*1, *a*), while moving along a straight line connecting these points. In case of diffusion, a particle is moving from the point (*x*1, *b*) to the point (*x*1, *a*) along a random trajectory, and the probability of its getting into point (*x*1, *a*) is not equal to 1 (Figure 2). Let it be assumed that at *t* = 0, the particle with the coordinate (*x*1, *b*) begins to move along the straight line connecting it with the point with the coordinate (*x*1, *a*). The one-dimensional process with *y* = *y*(*t*) for constant *x* = *x*<sup>1</sup> is first considered.

**Figure 2.** Diagram reflecting random particle trajectories (random walks) moving from the point (*x*1, *b*) to the point (*x*1, *a*).

Let at time *t* = *t*<sup>0</sup> the particle has coordinates of *y* = *y*0(*a* < *y*<sup>0</sup> < *b*), *x* = *x*<sup>1</sup> (0 < *x*<sup>1</sup> < *L*). The probability density *f(y)* and likelihood *P(y)* of this event are as follows:

$$f(y) = \delta(y - y\_0), \quad \int\_a^b f(y) dy = P(y = y\_0) \tag{1}$$

When *t* > *t*0, the coordinate *y* = ξ(*t*) is a random time function, then its probability density can be written in the following form [43]:

$$f(y) = <\delta(y) > \tag{2}$$

where the brackets <> denote the averaging over realizations.

Let us assume that ξ(*t*) satisfies a Langevin equation of the Brownian-motion type.

$$\begin{aligned} \frac{d\underline{\xi}}{dt} &= \nu(\underline{\xi}, t) + f(\underline{\xi}, t), \; \underline{\xi}(0) = y\_0 = \underline{\xi}\_0 \\ < f(\underline{\xi}, t)> = 0, \; B(\underline{\xi}, t, \underline{\xi}', t') &= < f(\underline{\xi}, t) f(\underline{\xi}', t')> = 2\delta(t - t') F(\underline{\xi}, \underline{\xi}', t) \end{aligned} \tag{3}$$

where *B*(ξ, *t*, ξ , *t* ) is a correlation function.

By (3), the differentiation of (2) with respect to *t* yields an equation for the probability density function of the particle's position:

$$\begin{aligned} \frac{\partial f\_t}{\partial t} + \frac{\partial}{\partial y} [ [v(y, t) + A(y, t)] f\_t(y) - \frac{\partial^2}{\partial y^2} [F(x, y, t) f\_t(y)] = 0 \\ A(y, t) = \left. \frac{\partial}{\partial y} F(y, y', t) \right|\_{y' = y} \end{aligned} \tag{4}$$

The Fokker-Planck-Kolmogorov differential Equation (4) should be solved with the initial condition (1). The physical meaning of the Equation (4) is in description of the stochastic processes for the motion of the coating's deposited particle in the gas phase from the cathode to the substrate. The coefficients in the equation at *ft*(*y*) are the particle transport and diffusion parameters, respectively. According to the Formulas (3) and (4), the equation terms with parameters *A* and *F* are specified by fluctuations of the field *ft*(*y*, *t*), affecting *v*(ξ, *t*), which is a deterministic function describing the convective particle transport. If *ft*(*y*, *t*) = *f*(*y*), i.e., is constant in time, then *A* and *F* are also constant in time, which corresponds to the constant diffusion coefficient *F* and convection *A* = 0. Cases when *f*(*y*, *t*) is non-stationary require a model specification.

The real process of the layer-by-layer formation of a coating is of stochastic nature due to the fact that the layer of particles which started their motion from a source plane does not remain planar, while the particles interact with each other and cannot simultaneously form microclusters during the deposition, with correlations formed between them. Such microconglomerates, being deposited on the substrate, will be considered as nuclei of clusters in the coating layer, which introduce randomness into the deposition process and make the coating structure heterogeneous, fractal. Depending on the mentioned random factors, the deposition time is a random variable for different particles.

Figure 3 exhibits an initial view of the surface of the Ti–TiN–(Ti,Cr,Al)N coating with a pronounced cluster structure.

**Figure 3.** Cluster structure of the surface of the Ti–TiN–(Ti,Cr,Al)N coating after deposition.

The model of a random process reaching the boundaries of its variation area makes it possible to estimate the probability of this event and the average time required for that [43]. Let *P*(*y*0, *t*) denote the probability of the event that the random process will reach the boundary *a or b* with the initial condition of *P*(*y*, *o*) = 0 and the boundary conditions of *P*(*a*, *t*) = *P*(*b*, *t*) = 1, then *P*(*y*0, *t*) satisfies the equation [43]:

$$\begin{aligned} \frac{\partial P(y\_0, t)}{\partial t} &= K\_1(y\_0) \frac{\partial P}{\partial y\_0} + \frac{1}{2} K\_2(y\_0) \frac{\partial^2 P}{\partial y\_0^2} \\ K\_1(y\_0) &= v(y\_0) + \frac{AF'}{4} K\_2(y\_0) = \frac{AF^2}{2}, \; F' = \frac{dF}{dy\_0} \end{aligned} \tag{5}$$

The drift coefficient *K*<sup>1</sup> and the diffusion coefficient *K*<sup>2</sup> are to be specified for every model. As is known [43], the Equation (5) can be solved only for some particular models. When the initial and boundary conditions are formulated, the probability to reach the boundary increases. Figure 4 exhibits in graphic form the probabilities *P*—to reach the boundary and *P*—not to reach the boundary when the normalization condition is satisfied.

$$P + \overline{P} = 1\tag{6}$$

**Figure 4.** Graphic form depicting the probabilities *P*—to reach the boundary and *P*—not to reach the boundary.

At Figure 4, *P* + *P* = 1 is the normalization condition:

$$P\_i = \frac{S\_i}{S}, \ \overline{P\_i} = 1 - \frac{S\_i}{S} = \frac{\overline{S\_i}}{S}, i = 1, 2, \dots \tag{7}$$

*Si* is the area of the rectangle *P(ti),*

*Si* is the area of the rectangle *P(ti)*.

Thus, *P(ti)* describes a change in the rectangle, which area is numerically equal to the probability of reaching the boundary (of the substrate), while *P(ti)* describes a change in the rectangle, which area is numerically equal to the probability of not reaching the boundary.

The time it takes for a particle to reach the substrate is a random variable over a set of realizations. Then to calculate the average time to reach < *t* >, a formula for mathematical expectation can be used as follows:

As far as *P*(*y*0, *t*) = 1 − *P*(*y*0, *t*), then

$$\begin{split} t < t > &= T = \int\_{0}^{\infty} t \quad \frac{\partial P(y\_0 t)}{\partial t} dt = -\int\_{0}^{\infty} t \frac{\partial \overline{P}}{\partial t} dt = \left\{ \frac{\text{integenu by parts}}{t = K, \frac{\partial \overline{P}}{\partial t}, dt = \partial \upsilon} \right\} \\ &= t \overline{P}|\_{0}^{\infty} + \int\_{0}^{\infty} \overline{P}(y\_0, t) dt = \int\_{0}^{\infty} \overline{P}(y\_0, t) dt, \text{ as far as } t \overline{P}(y\_0, t)|\_{t = \infty} \\ &= 0, \ t \overline{P}(y\_0, t)|\_{t = 0} = 0 \end{split} \tag{8}$$

Given (5), an equation for *T*(*y*0, *t*) can be derived from (7), as follows:

$$ = T = \int\_0^\infty t \frac{\partial P(y\_0, t)}{\partial t} dt = K\_1 \int\_0^\infty t \frac{\partial P}{\partial y\_0} dt + K\_2 \int\_0^\infty t \frac{\partial^2 P}{\partial y\_0^2} dt\tag{9}$$

Integration by parts:

$$\int\_0^\infty t \frac{\partial P}{\partial t} dt = tP|\_0^1 - \int\_0^\infty P(y\_0, t) \, \partial t = -1\tag{10}$$

$$\int\_{0}^{\infty} t \frac{\partial P}{\partial y\_{0}} dt = \frac{d}{dy\_{0}} \int\_{0}^{\infty} tP dt = \frac{dT}{dy\_{0}}\tag{11}$$

*Coatings* **2020**, *10*, 927

$$\int\_0^\infty t \frac{\partial^2 P}{\partial y \, ^2} dt = \frac{d^2}{dy^2} \int\_0^\infty t P dt = \frac{d^2 T}{dy^2} \tag{12}$$

Through adding, we get:

$$K\_2(y\_0)\frac{d^2T}{dy\_0} + K\_1(y\_0)\frac{dT}{dy\_0} + 1 = 0\tag{13}$$

The boundary conditions for (13) are:

$$T = 0 \text{ at } y\_0 = a \text{ or } y\_0 = b \tag{14}$$

The form of solution for (13) depends on the form of the coefficients *K*1(*y*0), *K*2(*y*0), i.e., the form of the equation of motion. If the drift and diffusion coefficients are constant, then the source of particles coincides with the *y*-axis, i.e., *y*<sup>0</sup> = 0, and the real and imaginary substrates are located symmetrically: *a* = −*b*, then

$$T = \frac{b}{K\_1} t h \left(\frac{K\_1}{K\_2} b\right) \tag{15}$$

Function graph (15) is depicted in Figure 5.

**Figure 5.** Dimensionless relationship for *T* and *b* at various values of *K*1/*K*2. 1: *K*1/*K*<sup>2</sup> = 10, 2: *K*1/*K*<sup>2</sup> = 1, 3: *K*1/*K*<sup>2</sup> = 0.1

Figure 5 depicts the dimensionless relationship for *T* and *b* at various values of *K*1/*K*2.

*b* is the distance from the cathode to the substrate surface (in conventional units).

Figure 5 proves that the larger the drift coefficient *K*<sup>1</sup> is, the shorter time to deposition is, while the larger the diffusion coefficient *K*<sup>2</sup> is, the longer the time to deposition is. A model of coating synthesis during the layer-by-layer deposition of particles on a substrate was considered. Under the ideal conditions, a layer of particles emitted by a plane source is deposited on a substrate as a layer, and then the layer thickness *h*(*t*) is a deterministic function. Under the real conditions, the coating thickness is a random function, as mentioned above, and it can be represented as follows:

$$H(t) = h(t) + v(t),<\upsilon> = 0,<\upsilon(t)v(t\_1)> = \frac{N\_0}{2}\delta(t - t\_1) \tag{16}$$

where ν(*t*) is a random function of white-noise type.

After cooling down, a layer transforms, and the layer thickness at the time *t* is a certain random process, represented in the following form:

$$H\_{\rm KOH}(t) = h\_{\rm KOH}(t) + \upsilon\_{\rm KOH}(t) \tag{17}$$

The coating synthesis process should end when *H*KOH(*t*) = *Q* is reached. The probability that the *H*KOH(*t*) function would reach the given thickness *Q* during the time *t* was considered. The quality of the coating is determined by the fact that its thickness in the final state should be a given constant value. The probability *P*(*y*0, *t*) of reaching the boundary *Q* is found by the formula:

$$P(y\_0, t) = 1 - \overline{P}(y\_0, t) = 1 - \int\_{-\infty}^{\underline{Q}} \overline{f}(y, y\_0, t) dy \tag{18}$$

where *f*(*y*0, *y*, *t*) satisfies the Equation [43]

$$\frac{\partial \overline{f}}{\partial t} = \alpha \frac{\partial}{\partial y} [y \overline{f}] + \frac{\alpha^2 N\_0}{4} \frac{\partial^2 \overline{f}}{\partial y^2} \tag{19}$$

with boundary conditions:

$$f(y\_0, Q, t) = f(y\_0, a, t) = \overline{f}(y\_0, b, t) = 0\tag{20}$$

The initial condition is as follows:

$$\overline{f}(y\_0, y, t\_0) = \delta(y - y\_0) \tag{21}$$

If *h*(*t*) = *Q*, i.e., a layer of the required thickness is deposited to the substrate, then the probability *P*(*y*0, *t*) has the following form:

$$\begin{array}{l} P(y\_0, t) = 2 \left[ 1 - \Phi \left( \frac{Q - y\_0 e^{-at}}{\sigma(t)} \right) \right] \\ \sigma^2(t) = \frac{N\_0 \alpha}{4} \left( 1 - e^{-2at} \right) \end{array} \tag{22}$$

We considered Φ as the probability integral (standard Gaussian distribution function) [44].

A model of fluctuations in the direction of the deposited particle motion was considered. As a result of such fluctuation, a particle can be deposited not in the intended place, and that can lead to non-uniform density of the coating. In the 0*xy* plane, the position of a particle is characterized by the values of *x*(*t*), *y*(*t*). Instead of the parameter *t*, a parameter *s* is introduced, defined as follows:

$$s = \int Vdt\tag{23}$$

where *V* is velocity of the particle motion.

At each point, particles deviate from the rectilinear trajectories of the motion along straight lines parallel to the *y*-axis (Figure 6).

**Figure 6.** Deviation from the rectilinear trajectories of the motion along straight lines parallel to the *y*-axis.

The probability density *f*(θ|*s*, 0) satisfies the following equation [45]:

sin θ <sup>∂</sup> *<sup>f</sup>* <sup>∂</sup><sup>θ</sup> <sup>=</sup> *<sup>B</sup>* 2 ∂ ∂θ sin θ <sup>∂</sup> *<sup>f</sup>* ∂θ *B* = lim Δ*t*→0 ε2/Δ*t*, ε<sup>2</sup> = *r*2< (Δϕ) <sup>2</sup> >, < Δ*r*.*r*Δϕ > *x* = *r* cos ϕ, *y* = *r* sin ϕ (24)

The value of *B* is related to the correlation function of the particle deposition rate

$$
\psi\_V(\rho) =  - \tag{25}
$$

through the formula as follows:

$$B = -\frac{-2}{\int\_0^\infty \frac{\partial}{\partial \rho} (\rho^2 \frac{\partial \Psi\_x}{\partial \rho}) \frac{d\rho}{\rho}}\tag{26}$$

The solution to the Equation (26) is being sought with the initial condition of:

$$f(\theta|0,0) = \delta(\theta)/2\pi\sin\theta\tag{27}$$

and looks as follows:

$$f(\theta|\mathbf{s} - \mathbf{s}\_1, \mathbf{0}) = \frac{1}{4\pi} \sum\_{n=0}^{\infty} (2n+1) P\_n(\cos \theta) e \frac{-n(n+1)}{2} B(\mathbf{s} - \mathbf{s}\_0) \tag{28}$$

where *Pn*(cos θ) are Legendre polynomials [45].

At *s* → ∞:

$$\lim\_{s \to \infty} f(\theta | s - s\_0, 0) = \omega\_1(\theta) = \frac{1}{4\pi} \tag{29}$$

It is convenient to describe the condition for a change in the particle velocity vector through cos θ, then cos θ = *P*1(θ) and the average value of the direction cosine of the velocity vector is calculated by the following formula:

$$<\cos\theta> = \int\_0^{\pi} \cos\theta V \sin\theta d\theta \tag{30}$$

The following differential equation < cos θ > can be derived from (24):

$$\frac{d<\cos\theta>}{ds} = -B < \cos\theta>\tag{31}$$

The solution of Equation (30) with the boundary condition of < cos θ >= 1 at *s* = 0 is of the form:

$$<\cos\theta> = e^{-Bs} \tag{32}$$

At *s* → ∞, < cos θ >= 0, and in accordance with (31), that means at a sufficiently large distance from a source to a substrate, all orientations of the particle motion are equally probable. For short distances:

$$<\cos\theta> = 1 - \frac{\theta^2}{2} = 1 - Bs\tag{33}$$

which follows that:

$$\text{H}^2 = 2\text{Bs} \tag{34}$$

The mean square of the particle deviation ρ from the initial (rectilinear) direction along the *y*-axis is calculated by the following formula [46]:

$$
\langle \rho \rangle = \frac{4s}{3B} - \frac{2}{B^2} \left( 1 - e^{-Bs} \right) + \frac{2}{2B^2} \left( 1 - e^{-3Bs} \right) \tag{35}
$$

It follows from (35) that when *Bs* 1, then:

$$
\sqrt{\langle \rho^2 \rangle} = \sqrt{\frac{B}{6}} s^{3/2} \tag{36}
$$

The mean square distance in a straight line to the particle being deposited is determined as follows:

$$
\langle r^2 \rangle = \frac{2s}{B} - \frac{2}{B^2} \left( 1 - e^{-Bs} \right) \tag{37}
$$

The Formula (37) at *Bs* 1 is reduced to the formula for the average square of displacement of a heavy molecule in light gas, obtained by Smoluchowski [47–50]:

$$
\langle r^2 \rangle = s^2 \Big( 1 - \frac{Bs}{3} \Big) \tag{38}
$$

The fractional degree in the Formula (35) indicates that the deposition of coating particles is described by the fractal dynamics.

An ideal model of coating deposition represents an emission of particles from a certain conventionally planar surface of a cathode with the same velocity (temperature) along a normal to a surface (plane). The particles move along the parallel trajectories and are being layer-by-layer deposited on the surface due to the complex planetary motion of the substrate surface (Figure 7a). The process is described by the convection equation.

**Figure 7.** Diagram for (**a**) the formation of the relief of the coating being deposited; (**b**) the predicted cracking for the multilayer coating.

Under the real conditions, while emitted, particles demonstrate velocity fluctuations in value and direction, and during the particle motion, the process becomes even more irregular. The cloud has density fluctuations, and diffusion of particles occurs, as a result of which the density and hence the stiffness and strength of the layer being deposited fluctuate over the surface area. Furthermore, the very surface of the substrate is in the process of the complex planetary motion (rotation). As a result, at the initial stages of the coating destruction, damages and then cracks arise in the places of the deposited coating, in which the density is lower [28–30,49]. In terms of mathematics, this is described by a model of Brownian motion of a particle (Equation (3)), and the probability density satisfies the Equation (4).

The analysis of the resulting model indicates that the very coating deposition process implies the formation of a fractal structure, which further contributes to the formation of stress concentration areas, which, in turn, stimulate the process of cracking. At the initial stages of the coating destruction, damages and then cracks arise in the places of the deposited coating, in which the density is lower. During the formation of the multilayer coating structure, mismatched fractal structures of the layers overlap each other, and thus the clustering effect is largely leveled out. In turn, the above reduces the likelihood of cracking in the coating structure and makes it possible to predict an increase in the performance properties of the coating [27–31,51–53] (see Figure 7b).
