Fractal Analysis on Surface Topography of Thin Films: A Review

The topographies of various surfaces have been studied in many fields due to the significant influence that surfaces have on the practical performance of a given sample. A comprehensive evaluation requires the assistance of fractal analysis, which is of significant importance for modern science and technology. Due to the deep insights of fractal theory, fractal analysis on surface topographies has been widely applied and recommended. In this paper, the remarkable uprising in recent decades of fractal analysis on the surfaces of thin films, an essential domain of surface engineering, is reviewed. By summarizing the methods used to calculate fractal dimension and the deposition techniques of thin films, the results and trends of fractal analysis are associated with the microstructure, deposition parameters, etc., and this contributes profoundly to exploring the mechanism of film growth under different conditions. Choosing appropriate methods of surface characterization and calculation methods to study diverse surfaces is the main challenge of current research on thin film surface topography by using fractal theory. Prospective developing trends are proposed based on the data extraction and statistics of the published literature in this field.


Introduction
Surface topography is one of the important indicators in thin film analysis and surface quality evaluation [1]. The quality of surface topography is generally characterized by surface roughness (such as R a and R q ), but many studies [2][3][4][5] have shown that the roughness value is obviously affected by the measurement scale, i.e., scale characteristics of roughness. Sayles et al. [2] found that the standard deviation of height distribution is related to the sample length. When the surface topography is measured under different scales, the calculated roughness values are variant. In general, the roughness value of the film surface will gradually increase with the expansion of the measured scale [5][6][7]. Thus, fractal theory can be used for more accurate analysis results at the microscale than roughness as an evaluation means.
Fractal theory is a nonlinear science, and mainly focuses on complex and random geometric objects that are self-similar/affine and scale-invariant, and it is widely used in the natural sciences and engineering applications [8][9][10]. In the field of surface engineering, Yehoda et al. [11,12] explicitly pointed out in 1985 that thin film surfaces showed fractal characteristics. The self-similarity is regarded as the whole to the local invariant after equal scale transformation in all directions. The self-affine is invariant when scaled by direction independent factors, i.e., the unequal scale transformation from the whole to the local. For films deposited under non-equilibrium conditions, the surface is usually self-affine rather than self-similar. The possible reason is that the newly absorbed atomic layers are correlated with their adjacent, deposited surface during diffusion process [13][14][15][16][17]. In addition, they showed that the scale dependence of surface roughness can be revealed by calculating the fractal dimension (FD). Unlike roughness, FD can be independent on measuring range within a certain scope, thus FD is more objective and reliable.
Scale symmetry is present in many objects and figures. It means that there is still a level of details can be observed even after zooming in or out, which is the same at different scales of observation [18]. FD is the key parameter of fractal theory and basic property of the fractal [19], which evaluates the complexity of a fractal objects. In other words, fractal analysis, compared to traditional surface topography characterization parameters, is not sensitive to resolution [20]. The value reflects the irregularity and degree of fragmentation on the surface. The FD value of a profile is between 1 and 2, while that of a surface is between 2 and 3. The larger the FD is, the greater complexity, irregularity, space-filling ability the topography has. When FD as a parameter is combined with the physical processes, it provides important descriptions and constitutes sufficient measures [21].
At present, the commonly used calculation methods of FD include the box counting (BC) method, the Higuchi method, the power spectrum density (PSD) method, the structure function (SF) method, the autocorrelation function (ACF) method, height-height correlation function, and the traditional roughness (TR) method. According to their calculation principles, their applicable occasions and result accuracies are different. The calculation of the Higuchi/PSD methods is generally based on the profiles of topography, while the TR method is usually based on multiple topography images. The BC/SF/ACF methods can be used for both profile analysis and image analysis, and the latter is used in most cases. However, it is worth noting that the calculation results of FD show certain differences due to different calculation methods, data selection ranges (scaling region), measurement equipment, etc. In this paper, the authors have reviewed the methods of both the fractal analysis on surfaces and the deposition techniques of thin films, summarizing the related results and trends in this field. Meanwhile, the mainstream fractal methods including reality issues are discussed in detail. Future development trends in this field are proposed and discussed.

Simulation Methods for Fractal Features
The surface topography of a thin film usually presents fractal property. However, the exact FD value of an actual film surface cannot be defined in prior, thus the accuracy of calculate FD value by using a certain fractal method is difficult to evaluate. Therefore, the accuracy evaluation of a concerned fractal method has to be carried out with the assistance of the ideal fractal surfaces with known FD values. The common solution is to create artificial topography by using the theoretical fractal functions, then the accuracy of a fractal method can be obtained by deviation of calculated FD value relative to the ideal FD. There are several functions for generating ideal fractal topography, including Weierstrass-Mandelbrot (W-M) function (the most frequently used one), Takagi function, Brownian motion function, etc.

Artificial Profile Generated by W-M Function
The widely used standard fractal profile and surface function is the W-M function, which is continuous but indifferent everywhere [22]. A series of sequence values with fractal characteristics, i.e., the one dimensional profile, can be generated by W-M function. The specific mathematical expression [23,24] is as follows, where D represents the ideal FD of generated sequence, 1 < D < 2; γ is a parameter determining the frequency density, which is greater than 1, and generally taken γ as 1.5 [25]; n represents the number of items accumulated by series; ϕ n represents the random phase, which is to avoid the phase repetition of generated sequence at any position [22,26]. It can be seen from Equation (1) that the W-M function is the superposition of a series of cosine functions with increasing frequency.

Artificial Surface Generated by W-M Function
When W-M function is used to generated two-dimensional ideal fractal surface, the equation is as follows, where z(x, y) is the height of topography at (x, y); L is the measurement scale of surface topography; G is the height scaling coefficient; φ n is a random phase; D i is ideal FD ranging from 2 to 3; M is the number of overlapping components. Different fractal surfaces could be obtained by setting related parameters, such as the measurement scale, FD, height coefficient, and frequency density of the W-M function. The surface generated by mathematical methods has a certain FD, which can be used to measure the accuracy of the calculated results by different methods. Feng [24,27] used W-M function to generate a series of fractal surface with known FD, compared the accuracy of different FD calculation methods, and then optimized the proposed methods based on the comparative results.
for x ∈ [0, 1] where ψ(x) = |2x − 2[x + 1/2]| and ψ n is an n-fold iteration by ψ. Tagaki function is composed of infinitely many piecewise linear continuous functions, and its slope is ±1. The Tagaki function is generally used in the generation of one-dimensional fractal profile [29,30], and the form is simpler than the W-M function.
(2) Fractal Brownian motion function After highly irregular trails of minute particles were found by R. Brown, Wiener proposed a mathematical model to exhibit the random and irregular paths in Brownian motion in 1923. Fractal Brownian motions are stochastic processes of non-stationary and are self-similar, and they are mainly used to describe the irregular topography of mountains, clouds, and landforms in nature [31]. Fractional Brownian motion of index-α(0 < α < 1) can be defined as a Gaussian process Y : [0, ∞] → R.

Surface Characterization Methods
For fractal analysis of thin films, the calculation results of the FD affects by the characterization of film surface, i.e., acquisition of surface data. The usage ratio of characterization methods for fractal calculation of thin film is shown in Figure 1 based on the data statistics of the literature, which can be referred in detail in Appendix A. It can be noticed that atomic force microscopy (AFM) is used most commonly to measure the surface topography of thin films. AFM is based on the closed-loop feedback of weak interaction force signals to obtain the surface topography, so it outputs image height data directly. Except the topographic images with height matrices obtained by AFM or scanning tunneling microscope (STM), the grayscale images obtained by using other instruments such as scanning electron microscope (SEM) have been reported to carry out fractal analysis. However, the FD value of such grayscale images can be different from those of AFM or STM, because the grayscale is not directly relevant to the height of a measured position. In some studies, the grayscale images output by SEM needs to be transformed into topographic data, and then used for subsequent calculations [32,33]. The directly acquired topography data of AFM is more suitable for fractal analysis than the converted data [34]. The calculation results of FD between the two methods are quite different. Therefore, it was necessary to pay attention to the source of topographical data before fractal analysis.

Calculation Methods for FD
With the wide application of fractal geometry, a variety of FD calculation methods have appeared, such as the BC, SF, ACF, PSD, TR method. According to different calculation principles, the above calculation methods can be roughly divided into scanline-based methods and image-based methods. The former are based on the profile, which referred to a standard height profile as (r, z). The latter are based on topographical images directly, which refer to a point plot of height matrix as (x, y, z). The data forms of profile and topographical image are presented in Figure 2. The procedures of calculating the ideal fractal profiles and surfaces with different methods are shown in Figures 3 and 4, respectively. The specific calculation equations and solutions of these methods are demonstrated in Section 2.3.   Liebovitch et al. [35] proposed the BC method, to calculate the FD of dynamic and iterated function systems. The profile point set is completely covered by boxes of size , and N B ( ) is the minimal number of boxes. Hence, its FD (D B ) can be expressed as follows.
In order to calculate FD easily and clearly, it could be written as follows, where C is a constant value. In the double logarithmic coordinate system, D B could be obtained by fitting the slope of the straight line with the least square method, as Figure 3b.
(2) Higuchi method The Higuchi method is a fractal method for time/space signals. For surface data, the FD is calculated by extracting profiles from surface topography. The specific calculation process [36] is: Step 1: Set a time series X(1) − X(N) of length N, and use the delay method to reconstruct the time series into the matrix X m k as: Step 2: Each curve length L m (k) is defined as: Step 3: The total sequence length is similar to the average of the lengths of the k delay-generated sequence curves, Step 4: Obtain the curves lg(L(k)) ∼ lg(1/k) for different k and L(k). If it is a straight line, it means that L(k) ∼ k exists. FD can be obtained by fitting in double logarithmic coordinates as Figure 3c. The Higuchi method has higher accuracy in calculating the FD of time series, but the setting of k max in the realization process is relatively vague. In realistic application, empirical or hypothetical values are often used, which lead to uncertain risks to the accuracy.
(3) Power spectral density method Surface topography can be considered as a series of signals in the spatial domain, and the power spectral density can describe the topography from a frequency perspective. Changing the measurement scale L of the surface topography is equal to changing the cutoff frequency f c at the view of the signal. If a surface topography is fractal, then the shape of the f (λ) spectrum remains the same, even if the observation scale is changed. Thus, the fractal curve S( f ) is as follows: where β = 5 − 2D S , and D S is FD of the profile. The power spectral density method can quantitatively describe the frequency spectrum, and provide abundant information [37]. However, it can be seen in the previous description that the calculation process is based on a linear topography rather than a surface topography. In order to calculate the FD of a whole surface, in the PSD method, it was common to disassemble the surface into a series of profiles, and calculate all their FD values. By adding 1 to the average FD of profiles, FD can be obtained. This calculation process greatly reduces the reliability of the results.

(4) Katz Method
The Katz method was proposed by Katz [38] in 1988 to evaluate the FD of electroencephalography (EEG). It can be regarded as a collection of a series of points (x i , y i ), i = 1, 2, 3, · · · N, where x i represents the gradual increase of time. The following is the calculation formula of FD (D) in the Katz method.
where n is the number of intervals between the midpoint of the waveform point column, i.e., n = N − 1, and represents the maximum plane diameter range of the waveform (i.e., the maximum distance from the starting point of the waveform to any other point), and L refers to the Euclidean geometric length of the waveform. The Katz method is considered to be an algorithm reflecting time-domain amplitude characteristics of signals, which has relevant applications in the fault diagnosis of bearings and gears.

Single-Image-Based Methods
(1) The Cube counting method The essence of cube counting (CC) and the BC method is the same. BC is used for two-dimensional profile and CC is used for surfaces. Some reports also refer to CC method collectively as BC. For the fractal surface topography, a cube (or sphere) with side length δ is used to cover it, and the minimum number of cubes required is N δ . Changing the side length δ of the cube obtains a series of N δ . When δ → 0, there is the relationship: N δ ∝ δ −D C , where D C is the FD. Therefore, the equation of this method is as follows: Limited by the resolution of the measuring instrument, the interval between the data points of the surface topography obtained by actual measurement cannot be infinitely small, resulting in a surface topography with limited resolution, and the ratio between the height and width of the surface topography is large. Therefore, in practice, a cuboid with a scale of L × L × h is usually used to cover the known surface topography. The value of L/h is usually selected from 100 to 1000 (according to the ratio of the height or the pixel width). The minimum number of cuboids required is N(L). The FD is calculated by the slope K of the log N(L) ∼ log(L) line fitted in a double logarithmic coordinate, as Figure 4b. The relationship between D C and the slope K is D C = K.
The BC/CC method is one of the most classic FD calculation methods and widely used in various fields. This method can obtain FD from a single image, but the disadvantage is that the calculation accuracy is relatively low.
(2) Structure function method In surface topography, the structural function of a spatial sequence with fractal features is Equation (15). (15) where · · · represents the averaging operation. The slope K of logS(τ) ∼ log(τ) is fitted by a straight line in a double logarithmic coordinate, as Figure 4d, and its FD is D = (4 − K)/2, where 0 < K < 2. For the surface topography, the FD of each profile must be calculated and the average valueD is calculated, and the value of the FD of the surface topography is D S =D + 1. The structure function method is suitable for self-affine profiles of engineering surfaces and can provide good results. However, as its FD must be transformed from the profiles to a surface, the accuracy is still not very high.

(3) Autocorrelation function method
For a spatial sequence y(x) in the surface topography, its autocorrelation function of the spatial sequence with fractal characteristics is: The structure function has the following relationship with the autocorrelation function: Therefore, the ACF method can be converted into a structural function method for calculation.

(4) Height-height correlation function
Vicsek et al. linked the spatial scaling behavior with time, and the concept of dynamic scaling appeared in 1985 [39,40]. During film growth, the interface width w as Equation (18), and the deposition time t exhibit power-law behavior. And obviously t is proportional to the film thickness l. Therefore, the power law can be written as Equation (19), and β is the growth exponent. The N × N height matrix of AFM can be analyzed by height-height correlation function, i.e., H(r). Since the actual film surface is continuous, the height obtained by AFM is related rather than independent. Along the scan direction, the height correlation can be obtained in Equation (20).
where m is number of pixel points used in calculation, d is the distance between two pixel points. FD is obtained by fitting a power law of the linear region of the double logarithm curve of H(r) and r [41]. As shown in Equation (21), the r value in a certain range can reflect the surface topography, such as self-affine and mounded thin film surfaces [42]. For smaller r, the relationship between H(r) and r is exponential, and FD = 3 − α, where α is roughness index between 0 and 1. For mounded surface, r is larger and H(r) exhibits oscillatory behavior. Thus, the value of r, when H increase to 1 − 1 e , is defined as the lateral correlation length ζ. The variation of the lateral correlation with thickness obeys the power law as Equation (22), where z is the dynamic exponent. If the lateral distance between two points is more than ζ, the correlation of them is little.

Multi-Image-Based Method
Roughness is one of the most common parameters that characterize surface topography, but the roughness value is highly affected by measuring instruments, measurement scales, etc. Previous studies have reviewed that there is a strong correlation between roughness and measurement scale, which is associated with the scaling characteristics of surface roughness. Soumya et al. [43] studied the power function relationship between roughness, measurement scale and FD in a fractal surface, as shown in Equation (23).
R q is the root-mean-squared value of the height measured over the surface; L is the measurement scale, i.e., side length of the scanned area; H is the Hurst exponent. Thus, the FD can be calculated from the slope of a least-square regression fitting line in a log-log coordinate plot of R q -L.
When comparing the BC/ACF/SF/PSD/TR methods, Kulesza et al. [44] found that the first four can calculate the FD from a single image, but the accuracy of the calculation result is low. Traditional roughness (TR) method has a higher accuracy, however, it requires a series of images with different measurement scales, leading to a decrease in calculation efficiency. To deal with this shortcoming, roughness scaling extraction (RSE) was proposed [24] based on the scaling effect and the robust of roughness against low image elements [45].

Roughness Scaling Extraction Method
The RSE method was found to be effective to analyze thin film surfaces [24], the morphological filtering effect of AFM measurement [27] and time-series signals such as EEG [23]. Similar to the TR method, the RSE method is also based on the scaling characteristics of surface roughness shown in Equation (23) and Figure 5; otherwise, the robustness of roughness at low image-resolution is another vital theoretical basis that ensures the validity of the segmentation process. Therefore, the RSE method only requires a single topographical image of the concerned surface. The sub-images at different scales are obtained through segmentation of the original image. Due to the inclination and waviness of the sample surface, the roughness value calculated directly based on the sub-images is significantly higher [5], so the flatten modification would to be carried out for all the sub-images to obtain a more accurate relationship between R q and L.
The flatten modification procedure firstly involves fitting the sub-image matrix and then subtracting the fitting matrix from the original one, in order to remove deviations of the sample surface due to large-scale bending, tilting, etc., so that the high-frequency details of the topography surface is enhanced. The overall procedure of RSE is shown in Figure 6. There are two methods of flattening modification: scanline flattening (RSE-f) for per-row/column data and surface planarization (RSE-p) for image matrix data. The order of flattening modification includes 0, 1, 2, etc., depending on the order of the utilized fitting function. RSE-f1 is used in the following because of its balanced performance of both high efficiency and accuracy with a mean relative error of about 0.64% [24]. The RSE method can not only accurately calculate the FD of a surface, but also acquire a reliabe R q -L curve during the process, which can quickly obtain the roughness of the surface at other scales.

Mechanism and Techniques of Thin Films
Many thin film surfaces, machined surfaces [2,[46][47][48], and other treated surfaces have been shown to have fractal properties. In the field of surface engineering, the surfaces have different degrees of fractal topography, which leads to different performance of the functional film [37,49]. In the field of thin films, different substrate surfaces and preparation methods lead to different fractal micro-structures. In the field of tribology, the fractal characteristics of surfaces affect the state of friction and wear [50,51]. Many theoretical and experimental studies have shown fractal features in various thin film structures. The ratios of the preparation methods in the literature are shown in Figure 7. In previous publications, the applications of fractal analysis obtained more valuable findings than traditional roughness characterization, which was helpful to study more comprehensive surface quality assessments and the study of topography formation and film growth mechanisms. The applications of fractal analysis in various thin films using different preparation methods are introduced in this sections. The films investigated in this review were prepared with physical vapor deposition (PVD), chemical vapor deposition (CVD) and non-vacuum methods, and their surfaces were analysis with fractal methods.

Surface Growth and Factors
In surface growth, clustering and diffusion of atomic or molecule affect surface topography of thin films at the micro scale. Roughness index and dynamics index can be calculated by surface growth model. The surface growth model is expressed by nonlinear partial differential equation, and partial differential of time-variable is expressed by space-variable partial differential derivative. Edwards-Wilkinson equation [52], Kardar-Paris-Zhang equation [53], Lai-Das Sarma equation [54], and Kuramato-Sivashinsky equation [55] can express the process of particle deposition and accumulation into film under chaotic dynamics. However, the order and variable of the expression are different, so the numerical solution of nonlinear partial differential equations in 1 + 1 dimension needs further discussion.
For Edward-Wilkinson equation, it can be expressed as Equation (24) For Kardar-Paris-Zhang equation, it can be expressed as Equation (25).
For Lai-Das Sarma equation, it can be expressed as Equation (26).
The two common statistical models of sedimentary processes are random deposition and ballistic deposition. In the two deposition methods, fractal model can simulate surface height of thin films and reveal growth mechanism to some extent. Ghosh et al. found that the roughness index of the ballistic deposition model calculated by FD fitted with the calculated value of the Kardar-Paris-Zhang continuum growth model [56]. Qi et al. verified the fractal characteristics in the film growth process in Kuramato-Sivashinsky model [57].
The deposition process and final morphology of films are also influenced by shadowing and diffusion effect. The cause for shadowing is that the higher part of the angular surface topography blocks the lower part, resulting in different particle fluxes received at different locations on the surfaces. The film grows in columns, eventually forming mound on surface [58]. In oblique incidence flux and substrate motion, different deposition angles can affect deposition effect of shadowing [59]. When the fractal surface is affected by different degrees of shadowing, FD is also affected. In the shadowing effect model, the shadowing function is a function of incident angle θ, RMS slope S dq and observation length L. The RMS slope is a function of FD, so the shadowing function is related to θ and FD. The rough surface with larger FD has a more significant masking effect [60,61].
To quantify a fractal surface, except FD, there are some advanced fractal parameters can be used as surface entropy, fractal succolarity and fractal lacunarity. Surface entropy is used to describe the height distribution homogeneity [62], which influences the chemical structure and surface tension of surface at micro or nanoscale [63]. The value of surface entropy varies between 0 and 1, with 0 representing maximum heterogeneity. Being close to 0 of surface entropy may indicate topographic irregularities and problems in adhesion [64]. Fractal succolarity relates the intercommunication in a surface [62], which can be seen as the degree of surface percolation. For the film surface, it is assumed that there is fluid entering the upper band of the surface, and the channels composed of similar points can assist fluid to the lower band. Penetration into the inner layer through such a process is the phenomenon of percolation [65,66]. Fractal lacunarity quantifies the distribution of gaps in space as a function of scale [67]. A lower fractal lacunarity indicates that the gaps are uniform, and the surface can be indirectly considered to be homogeneous [63].

Electron Beam Evaporation
Electron beam evaporation (EBE) heats and vaporizes coating material by accelerated electrons. The kinetic energy of the electrons is mostly converted into thermal energy, thereby achieving evaporation. This method has a high evaporation rate and a high purity of the film, and is suitable for fabricating thin films with a high melting point.
Pandey et al. [68] studied BaF 2 thin films prepared on three different substrates of Al, Si, and glass, and they calculated the average roughness, interface width, and surface FD. It was found that the FD values of the film on a glass substrate in the horizontal and vertical directions were close, indicating isotropy. The FD of the BaF 2 film on the Si substrate was higher than the others, indicating that the surface topography was highly complex. The variation of Hurst exponent and crystalline size with thickness on different substrates was systematic studied, which contributed to growth mechanism understanding. Yadav et al. [20] also studied BaF 2 thin films, and analyzed topographical changes over the growth of thin films through roughness R a , interface width, and FD. Among them, the FD and multifractal approach were mainly used as complexity measures to demonstrate the geometrical and physical characters of ions irradiated surface, which were excellent tools for characterizing surfaces.
The general finding in those reports is that fractal structure actually exists in the EBE films. FD can exclude the influence of the scale effect, which is different from roughness. Hence, using fractal analysis is more direct to study the relationship between the surface topography and substrate material, preparation parameters, etc.

Pulsed Laser Deposition
Pulsed laser deposition (PLD) irradiates target with a focused laser pulse, generating a strong laser spray plume. Evaporated atoms in the melt beam are then deposited on the substrate and form a thin film. Through the PLD method, the crystalline state (size distribution and shape of nanocrystals) and the adhesion to the substrate can be controlled by changing the target-to-substrate distance, background pressure, substrate temperature, laser flux and other parameters. As reported, too many parameters could affect the PLD surfaces, which could be effectively quantified and compared by fractal analysis.
Soumya et al. [43] studied ZnS thin films prepared with the PLD method at different annealing temperatures and analyzed the surface topography based on AFM image data. A strong correlation was found between the annealing temperature, surface roughness and FD. Singh et al. [69] studied multiferroic BiFeO 3 films prepared on Si, STO, and sapphire substrates, and they calculated the FD of the topography profiles in rows and columns respectively. It was found that all the samples shown the anti-persistence behaviour by fractal analysis of Hurst exponent, and the anisotropic level on different substrates can be suggested by FD value.

Magnetron Sputtering Deposition
As one of the PVD methods, magnetron sputtering deposition is commonly utilized in the preparation of thin films. The theory is to bombard the target with hundreds of eV ions accelerated in an electric magnetic field. Target atoms obtain energy, and are sputtered into the gas phase. The sputtered atoms re-agglomerate on the substrate and then form a thin film. In the meantime, the secondary electrons generated during sputtering are bound by the electric-magnetic field in the plasma region on the target surface, so a large amount of ions are continuously ionized to bombard the target to improve the deposition rate. In 1999, high power pulsed/impluse magnetron sputtering technology was proposed by Kouznetsov et al. [70], which combined the advantages of high intensity arc technology and magnetron sputtering, attracting attention for its high metal ionization rate. High power impluse magnetron sputtering can improve film quality, reduce surface roughness [71], intensify adhesion of film and increase hardness and strength [72].
Wang et al. [37] applied the fractal theory to analyze the surface topography of a Cu-W thin film during the growth process, and explored the relationship between FD and the resistivity of the thin film. It was found that the surface roughness of the thin film was caused by high-frequency components. As the deposition time increased, the thickness and FD of the film increased. In addition, for the film surface of the uniform phases, the resistivity was positively correlated with FD. Similarly, Talu et al. [73] analyzed the surface topography of Ag-Cu thin films at different deposition times and found that as the deposition time increased from 4 to 24 min, the particle size and FD both increased and then decreased. The resistivity was then found to be negatively correlated with deposition time, which was ascribed to the increment of the contacts between particles. In these reports, FD can be used not only to evaluate surface complexity, but also to study its correlation with film properties.
The thin film growth is affected by the fabrication parameters and is shown in the analysis of FD. Talu et al. [74] studied the effect of annealing temperatures on the surface topography of C-Ni thin films, and analyzed the correlation between the surface texture characteristics and the length scale. Meanwhile, there were two corner frequencies in the profile for calculating the FD, and their values were positively correlated with tip radius and cluster size respectively. Fang et al. [75] analyzed the nano-friction and wear characteristics of ZnO thin films deposited at different sputtering powers. Moreover, the higher sputtering power led to smoother surfaces, lower FDs, a more pronounced preferential orientation, and larger nano-wear rate.
Feng et al. [76] used the RSE method to perform fractal analysis on MgO thin films prepared by the energetic partical self-assisted deposition (EPSAD) method. The FD values of MgO films under different preparation were compared to reveal deposition mechanism during the EPSAD method. Talu et al. [77] focused on zinc oxide film, calculated multiple fractal characteristics of the 3-D texture on the surface based on AFM measurement, and proved that there was a significant correlation between the fractal features and the surface texture. Talu et al. [32] also studied the surface texture of Cu/Co thin films and established a set of friction mathematical models including S q , kurtosis, FD, and corner frequency, from nanoscale to macroscale, so that the friction properties of the films could be controlled and modified. To understand material science such as crystal growth, defects and phase transformations, fractal analysis contribute to explaining the relationship between topography and material performance.
Based on the above reports, FD as the complexity of the film surface is affected by the film preparation parameters. On the other hand, it can be related to the application performance. As an intermediate parameter, the selection of preparation parameter could be optimized according to the FD trend to obtain the desired application performance.

Chemical Vapor Deposition
CVD refers to the method of depositing thin films by a surface chemical reaction with the raw reaction materials being gaseous and at least one of the products being solid. Compared with PVD, CVD applies a vapor phase reaction, by which the surface can form films with a high production efficiency and performance consistency. The preparation process and parameters of CVD and PVD are different in principle, but the fractal properties of thin films fabricated with both methods have been reported.
Talu et al. [78] studied the three-dimensional surface texture of amorphous hydrogenated carbon by the different power of radio frequency-plasma enhanced CVD, and the analyzed parameters includes FD, pseudo-topothesy, anisotropy ratio, skewness, kurtosis, and roughness. Those function allowed for a reasonable prediction of space patterns. The statistical, fractal, and functional surface properties of films contribute to a characterization of micrometer-and nanometer-scale textures. Haniam et al. [79] synthesized a fractal thin film of cobalt oxides by using CVD, and its density and crystallization were affected by various factors, such as laser irradiation time and the gas flow ratio. This meant that the fractal structures could be controlled in CVD process. Fractal structures described by diffusion-limited aggregation model and compared with the FD of the growth pattern by box counting. Li et al. [80] built a 2D diffusion-limited aggregation model, which was based on a typical model of fractal theory. By this model, the CVD growth of various fractal-morphology 2D materials were precisely controlled. Fractal-growth theory applied to CVD-growth process provided favourable predictions for fabrication parameters.

Electroplating Method
Electroplating is the process of plating a thin layer of metals or alloys on the surface of a substrate using electrolysis, which can enhance the corrosion resistance of the metal, increase hardness, and prevent wear.
Nasehnejad et al. [81] calculated the surface FD of Ag films with different thicknesses. The film thickness increased from 250 to 750 nm, and the FD increased from 2.57 to 2.64, indicating that, during the deposition process, as the deposition time increased, the irregular shape of the surface became more prominent. Naseri et al. [82] analyzed TiO 2 thin films prepared at different current densities and found that extremely low or high current densities cause anisotropic structures of the films, but medium current densities would cause uniformly decorated tubes, and the FD was also lower. Gholamreza et al. [83] studied Ag films obtained from different deposition current densities and substrate rotation speeds. The higher the rotation speed, the larger the FD value was, and the FD could be used to describe the change in the overall topography along the growth direction. Fractal analysis provides the lateral development of surface features [42], and it is useful for understanding the growth and scaling of thin films during formation.

Sol-Gel Method
Sol-Gel method is a process of converting sol, a colloidal suspension of solid particles in liquid [84], into thin films by coating and heat treatment [85]. The unique features of the sol-gel process, including the chemical reactions, the ease of doping, and the low synthesis temperature, were conducive to the preparation of industrial coatings.
Ghosh et al. [86] studied the surface microstructures of sol-gel spin coated ZnO thin films in varying precursor quantified by Higuchi method. There was a certain correlation between grain size and the Hurst exponent, and larger average grain sizes were attributed to higher Hurst exponent values. Pandey et al. [87] analyzed three parameters of Al-doped ZnO thin films by AFM images, were thickness, refractive index and FD. The root mean square interface and FD used to quantify the surface topography decreased with increasing thickness. It was found that increasing film thickness, corresponding to densification and a growing number of hillocks, may contribute to the refractive index.

Parameters Related to FD
In traditional surface analysis, roughness is often used to evaluate surface quality at the micro-scale. However, due to the scaling characteristics of surface roughness, it is difficult to compare the roughness at different scales. A surface with many fine structures and a high degree of fragmentation might have a lower surface roughness value, but the actual corresponding quality does not meet the engineering requirements. Therefore, it is difficult to characterize the essential characteristics of micro-engineered surfaces based only on roughness. Previous studies have found that fractal phenomena exist during the nucleation and growth processes of thin films at the microscopic scale. Thus, fractal theory contributes to a comprehensive characterization on surface topography.
As a scale-independent parameter reflecting the surface microstructure properties and surface quality, FD has a strong correlation with preparation process parameters and topographical parameters.
In terms of the preparation process, common process parameters include deposition time, sputtering power, and annealing temperature. Trends reported in some articles were summarized, but corresponding to specific preparation methods and topography, researchers still have to use FD to explore the rules. As the deposition time increases, the irregular shape of the surface becomes more prominent, and the value of FD increases accordingly [37,63]. As the sputtering power increases, the preferential orientation of the grains becomes more pronounced [75,76], and the FD decreases. As the annealing temperature increases, the atomic mobility increases and the FD decreases [27,43,88].
In terms of topographical parameters, the change in FD is consistent with cluster size formation, which is usually interpreted as follows. The formation or aggregation of island-like structures on the surface is often manifested by an increase in cluster size, which causes the rise of surface complexity, and thus the value of FD increases. Moreover, in the process of calculating the FD, the corner of the R q -L curve is usually related to the tip radius and the size of the clusters. In addition, the FD of the profiles in different directions of the topography can be used to quantify the anisotropy of the film.
At present, the relationship between fractal parameters and film surface characteristics and topography cannot be predicted [42]. It is necessary to conduct fractal analysis of specific films and attempt to explore the relationship between FD and performance parameters. Though fractal results, the fabrication parameters can be modified.

Fractal Study of Various Thin Films
In investigations of thin films prepared by different technologies, previous studies have repeatedly found that the microstructure of the topography appearing in the growth process of the films exhibits fractal properties. The FD can be used to measure the complexity of the topography and explore the effect of preparation conditions on topographical quality [89,90]. In recent years, different methods have been used to calculate the FD of various thin films, and the results are summarized in Appendix A.
For the FD value, firstly, it limited by the definition range. The FD for profiles is between 1 and 2, and it is between 2 and 3 for the topography surface. Secondly, the FD value difference of the same type of film under different preparation conditions is almost within 0.5. The FD value directly reflects the irregularity of the surface topography, and the deposition state during film preparation can be characterized by the Hurst exponent. Hurst exponent is also often used to evaluate the randomness of the time series. The relationship between H and D is shown in the Equation (28).
(1) When 0.5 < H < 1, the time series has a long-range correlation [81], i.e., an increasing (decreasing) trend over a certain period of time, and an increasing (decreasing) trend in the next period, and the closer H is to 1, the stronger the correlation is.
(2) When H = 0.5, the time series is irrelevant and is an independent random process; i.e., the current state does not affect the future state. The change at any time step is uncorrelated with the change of previous or later steps.
(3) When 0 < H < 0.5, there is only negative correlation in the time series, showing an anti-persistent state [91], i.e., the time series is increasing (decreasing) in a certain time period, and then decreasing (increasing) in the next time period.
The value of H ranging 0.5 to 1 can be found if there is a long-range correlation (also known as long memory) during the preparation process [92], which can provide detailed information on cluster growth and surface formation regarding continuity.

Application of Fractal Calculating Methods
During the development of fractal theory, the calculation method of FD was constantly being updated, and was increasingly used in the topography of thin films. To investigate the utilization of different fractal methods in recent decades, a schematic diagram of the development in accumulative count of studies using a certain FD method over time is shown in Figure 8a. It can be seen that classic fractal calculation methods of BC and CC are the most commonly used in the field of surface topography analysis. Although the Higuchi/SF/PSD methods had been proposed as early as the 80s~90s, these methods are commonly used for fractal calculations of time series, such as signals. In recent decades, these fractal calculation methods have extended to the surface area and have been widely used in growth mechanism analysis and surface quality evaluations of various types of films. The utilization ratio of these fractal methods is shown in Figure 8b. Representative of traditional methods, BC and CC methods are the most widely used, and various new methods were developed and applied after 2010. The community of thin films is now embracing a significant uprising of fractal analysis on surface topography.
As a core indicator of fractal theory, the FD represents the complexity of fractal geometry. Thus, accurate calculations of FD play an important role in fractal analysis. In order to compare the accuracy of the results of different FD calculation methods, it is necessary to generate surfaces with an ideal FD for verification calculations by mathematical methods. As a general approach, a series of profiles and surfaces with ideal FDs are generated with the W-M function, and FD is calculated using the above methods. In Figure 9, the accuracy of different method is compared. The relative errors of the results calculated with each method significantly different under each ideal FD. For profiles, the calculated results of BC method and KATZ method are significantly different from the ideal FD, even though they are widely used. The results obtained by the Higuchi method are relatively stable when the FD is 1.4∼1.9. For surfaces, in the BC method, when the ideal FD is large, the obtained FD is obviously different from the ideal FD. In the ACF and SF methods, when the ideal FD is low, the calculated FD and the ideal FD are significantly different. In the PSD method, the calculated FD result is relatively accurate only when the ideal FD is larger than 2.4. Therefore, traditional methods cannot provide calculation results of the FD in the full range of the ideal FD (2.0∼3.0) with high accuracy. The RSE method is the most accurate in most cases.

Reality Issue for Fractal Analysis
As early as 1998, Avnir [93] noticed that the calculation of FD usually requires the experiment data with a large magnitude span of scales, which centers around 1.3 orders of magnitude in a histogram of this perspective. Based on the extracted data of relevant literature in this study, the maximum and minimum values of scales are summarized in Figure 10. By using the division operation on these values, the magnitude span of scales for fractal analysis of thin films was obtained, as shown in the distribution histogram of Figure 11. The average order of magnitude span in the past two decades in the field of thin films was 1.86. It can be speculated that the enlarging magnitude span relative to the report by Avnir [93] in 1998 indicated an advancement of fractal analysis, which can be attributed to the improvement in surface measurement technology. Currently, the further enlargement of scan scale of AFM is expected, but probe wear could be accelerated due to excessive data collection; thus, a reasonable selection of sampling length could reduce equipment wear, accelerate FD calculations while ensuring the accuracy [94]. There must be further improvements in surface measurements in the next decade.  However, the limitation of a further enlarging of the span might be not technical in the future. The difficulty lying ahead could be the intrinsic characteristics of thin films. As shown in Figure 10, the common maximum scale is around 10,000 nm, while the minimum is around 1 nm. Generally, such a maximum scale corresponds to the typical size of clusters of thin films. If a larger scale is to be investigated, the surface topography is dominated by other factors such as waviness, which are usually dependent on the fabrication process other than cluster growth. The minimum scale is similar to the typical size of stabilized nuclei. If a smaller scale is to be investigated, the individual existence of atoms or ions would be considered, and even the sub-atom structures might be taken into consideration, and these are both quite different from the atomic diffusion process during the nucleation and growth of thin films.
The results and discussion above could be related to the reality issue for fractal analysis, i.e., the scaling region. Fractal analysis can theoretically deal with self-similarity or selfaffine features with boundless scales, but the real features in nature are all within a certain range of scales. In 1989, Yokoya et al. [95], when they studied the geographical fractal phenomenon, put forward the argument that fractal feature exists only in the scale range with upper and lower bounds, not a global scope, and the simulation function of fractal Brownian motion was used to derive the specific size of range. The method of cutting the scale-free region based on a genetic algorithm proposed by Fei et al. [96] was applied to the fractal analysis of ground surface. The summation of the average variance can be applied for searching the optimal position of the initial and terminal points. Zuo et al. [97] found that the linearity of the scaling region increased with the sampling length, and the minimum sampling length was greater for a smooth profile than that for a rough profile. The advantage of scaling region interception to improve accuracy was also found in a fractal analysis of EEG signals [98].
For real features in nature, measurement data within the scaling region rather than the global region can invariability exhibit fractal properties with the scale with a high efficiency. When calculating the FD, the linear part in the double logarithmic coordinates can be considered as the scaling region, as shown in Figure 5. Therefore, the identification of the scaling region should play a crucial role in fractal analysis.

Prospective Developing Trends
FD is a measure of surface irregularities and discontinuities. Fractal analysis can be associated with cluster size [73][74][75]86], deposition state [37,63,76], topography [69,88], material composition estimations [99] and other parameters to obtain more objective results than the traditional roughness characterization, which is favorable for further interpreting the growth mechanism of thin film surfaces [68,79,83,100]. The specific value and variation trend of FD can play an important role in the comparative investigations of thin films under different conditions. The optimization of calculating FD has also become a principal part of the development of fractal theory.
As the understanding of fractal theory gradually deepens, the calculation of FDs evolves on the basis of traditional methods, and both the calculation efficiency and accuracy are going to be improved. The comparison of mean relative error in this paper assists researchers to utilize the fractal methods with high accuracy.

Multifractal Analysis for Thin Films
In order to obtain more detailed fractal descriptions, it is necessary to increase the parameters that describe different fractal sets, so the multifractal theory has been developed and applied widely in last decade. A multifractal is a staggered superposition of a large number of single fractals with different dimensions. It is generally characterized by the multifractal spectrum f (a), revealing the complexity and singularity of a class of topography, where a determines the strength of the singularity property, and f (a) describes the density of the distribution. As an evolution based on single fractals, the multifractal theory is developing fast [101][102][103][104][105] and becoming a powerful tool for surface analysis. It would be promising to open up new research fields for the application of fractal theory and exciting research achievements. Multifractal detrended fluctuation analysis can determine higher-dimensional fractal and multifractal scaling behavior in time series, which can better filter out the trend components of evolution, retain the fluctuation components, and assist researchers to extract effective information [106,107]. Moreover, with the emerging of more indicators for surface analysis (including FD, multifractal parameters and traditional parameters), it would be convenient to establish a practical scheme of machine learning, such as an artificial neural network model [108], which has demonstrated high accuracy and efficiency in fractal analysis.

Fractal for Surface Mechanism and Properties
Fractal characters are found in natural topography [8], crystals [109], plants [110], retinal blood vessels [111], etc. For the field of surface engineering, research on the growth mechanism of fractal structures and its impact on the practical applications will be further explored systematically.
In terms of scientific theoretical research, the micro model of the fractal growth mechanism has been established to explain the surface characteristics of thin films [80]. However, whether the theoretical model of fractal growth is consistent with the actual film deposition process has not been systematically studied. The details of film nucleation and cluster growth need to be characterized with more powerful instruments for a comparison with the established models. Such work could further explore the principle of film deposition and reveal the formation mechanism of natural substances from a fractal perspective.
In terms of practical application, many of above-mentioned previous studies have linked fractal characteristics with the application properties of specific film surfaces, such as gas sensing characteristics [49], refractive properties [87], wettability [102], capacitance behavior [112], and resistivity performance [37]. Since the artificial surface of fractal structure can characterize rough surface profile [113], which is random but is also multi-scale and self-similar [25], it is also used in modeling research on the effect of surface topography on performance [114,115]. For studies related to rough surface properties, such as contact process [116], surface liquid film evaporation [117], and post-treatment surface structure [118], fractal analysis can be used to model surfaces. Then the FD as surface evaluation parameter, combined with other physical performance parameters. As a parameter for evaluating the surface complexity, FD is a functional carrier that supports researchers to quantitatively associate surface topography with both fabrication parameters and application performance. In the future, such a functional carrier will be used in more researches on scientific questions and engineering applications.

Summary
This article presents an overview of main methods of fractal analysis, ranging from calculation methods for FD to simulation methods of fractal characteristics, and a variety of thin film fabrication techniques. In fractal analysis, the simulation of arificial fractals for evalating the accuarcy of fractal methods, the source of surface data, and the calculation method of FD, are introduced in detail. For calculating FD, methods include scanline-based methods (including BC, Higuchi, PSD, and Katz), single-image-based methods (including CC, SF, ACF and height-height correlation function), and the RSE method based on the TR method. By comparing their usage and accuracies, it provides a reference for researchers to choose appropriate fractal calculation methods. Physical vapor deposition methods, such as electron beam evaporation, pulsed laser deposition, and magnetron sputtering deposition, chemical vapor deposition, and non-vacuum methods, such as electroplating and sol-gel are described. The relationship between film deposition parameters and FD was analyzed, along with the current developments in fractal calculation methods. FD as a topography parameter can feedback the process parameters and application properties in the comprehensive characterization. In nature, the scale range in which the fractal structure exists is not global, so it is a reality issue to determine the scale range selected for fractal calculation in thin films. The scaling region, which would affect the efficiency and accuracy of calculations, was discussed by investigating the literature over the past 2 decades. Compared to the result of 1998, the average order of magnitude span expended to 1.86, which is predicted to be further improved. Through the analysis of literature data on fractal calculations of the film surface, future development trends in this field are proposed including multifractal analysis, surface modeling and practical application. In addition to offering a broad introduction to the common methods of film deposition and fractal analysis, this review can provide an extensive set of references and statistical data that can be used by researchers in the fields of both thin films and fractal analysis.

Data Availability Statement:
The data presented in this study are available from the corresponding author with reasonable request. Moreover, the codes of analysis and simulation on fractal profiles and surfaces have been established, and the collaboration investigation by using them is welcome by the authors.

Acknowledgments:
The authors would like to thank Shaozhu Xiao, Xiangsong Zhang, Binbin Liu, Shaocong Wang, Junlong Huang, Jun Li, Yousheng Xia and Zengxian Ma of Tsinghua University for their contributions to fractal studies in our lab.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: